undet-coeff.mws

>    restart;

Undetermined coefficients

About this worksheet

Author and Date

Matthias Kawski
Arizona State University

kawski@asu.edu
 http://math.asu.edu/~kawki

Original version: July 2005. MAPLE release 8

>   

Content, Purpose and Use

This worksheet develops a hands-on approach to the method of undetermined
coefficients to solve inhomogeneous costant coefficient linear differential equations
whose right hand sides are combinations of polynomial, sinusoids, and exponential
functions.

Except for some fancier commands to use MAPLE to extract coefficients of "like

terms" and for some display purposes an attempt is made to keep the MAPLE
commands simple and accessible to relative novices.

The main idea is invite to some "trial-and-error" approach, lots of experimentation.
of course, the main objective in the background is to develop the thinking in terms
of differential operators (as opposed to differential equations), encouraging the

thinking in terms of inputs and outputs. This should in particular strengthen the
understanding of the role of linearity in this method amd setting.

Updates and log of modifications.

None yet.

>   

From a naive trial and error approach to linear algebra

Introduction

>   

Define a differential operator

Define a linear differential operator -- here with constant coefficients.
Note that for this to be useful, one needs to substitute an expression dependingt on t for the
formal parameter u.

>    L:=u->diff(u,t,t)+4*diff(u,t)+3*u;

L := proc (u) options operator, arrow; diff(u,t,t)+4*diff(u,t)+3*u end proc

Note that MAPLE does not accept the formal parameter IN THE DEFINITION  of this
procedure to be called u(t) -- however, when using the procedure, it should be called using
either an explicit formula depending on t, or using y(t).

>    badL:=u(t)->diff(u(t),t,t)+4*diff(u(t),t)+3*u(t); 

Error, symbol or symbol::type expected in parameter list

A few brief tries

>    L(y(t));

diff(y(t),`$`(t,2))+4*diff(y(t),t)+3*y(t)

>    L(exp(t));

8*exp(t)

>    L(exp(r*t));

r^2*exp(r*t)+4*r*exp(r*t)+3*exp(r*t)

>   

The problem: solve an inhomogeneous differential equation

The problem is to find the general solution of the inhomgeneous linear differential equation

>    IDE:=L(y(t))=g(t):
IDE;

diff(y(t),`$`(t,2))+4*diff(y(t),t)+3*y(t) = g(t)

As we saw earlier, the general solution of this inhomogeneous DE is the sum of the
general solution of associated the homogeneous DE and any one solution of the

inhomogeneous DE.

>    HDE:=L(y(t))=0:
HDE;

diff(y(t),`$`(t,2))+4*diff(y(t),t)+3*y(t) = 0

>    homsol:=dsolve(HDE,y(t)):
homsol;

y(t) = _C1*exp(-3*t)+_C2*exp(-t)

The problem thus is reduced to finding just one solution of the inhomogeneosu DE.
The method of undetermined coefficients works particulalry well with constant coefficient
linear differentail equations whose right side is a linear combination of (products of)
exponentials, sines and cosines, and polynomials.

>   

The first example: linear combination of exponentials

Let us consider a simple example of a right hand side:

>    exa1:=L(y(t))=3*exp(-2*t)-5*exp(-4*t):
exa1;

diff(y(t),`$`(t,2))+4*diff(y(t),t)+3*y(t) = 3*exp(-2*t)-5*exp(-4*t)

Start w/ trial and error, trying different functions which L may map to terms similar
to those on the right hand side:

>    L(exp(-2*t));

-exp(-2*t)

>    L(-3*exp(-2*t));

3*exp(-2*t)

>    L(exp(-4*t));

3*exp(-4*t)

>    L(-5/3*exp(-4*t));

-5*exp(-4*t)

>    got_it:=-3*exp(-2*t)-5/3*exp(-4*t);
L(got_it);

got_it := -3*exp(-2*t)-5/3*exp(-4*t)

3*exp(-2*t)-5*exp(-4*t)

Simple trial and error, and taking advantage of the linearity does the job.

A little more formal approach is to start with "UNDETERMINED COEFFICIENTS"
and equate the result to the right hand side of the differential equation

>    L(A*exp(-2*t)+B*exp(-4*t))=rhs(exa1);

-A*exp(-2*t)+3*B*exp(-4*t) = 3*exp(-2*t)-5*exp(-4*t)

In this case it is immediate to read of the coefficients which yield a solution. The general
solution then is

>    sol1:=y(t)=rhs(homsol)+got_it:
sol1;

y(t) = _C1*exp(-3*t)+_C2*exp(-t)-3*exp(-2*t)-5/3*exp(-4*t)

check it!

>    L(rhs(sol1))-rhs(exa1);;

0

>   

The second example: polynomial right hand side

>    exa2:=L(y(t))=15*t^4-7*t^2:
exa2;

diff(y(t),`$`(t,2))+4*diff(y(t),t)+3*y(t) = 15*t^4-7*t^2

Start w/ trial and error, trying different functions which L may map to terms similar
to those on the right hand side:

>    L(t^4);

12*t^2+16*t^3+3*t^4

>    L(5*t^4);

60*t^2+80*t^3+15*t^4

>    L(t^3);

6*t+12*t^2+3*t^3

>    L(5*t^4-80/3*t^3);

-260*t^2-160*t+15*t^4

>    L(t^2);

2+8*t+3*t^2

>    L(5*t^4-80/3*t^3+253/3*t^2);

-7*t^2+1544/3*t+506/3+15*t^4

>    L(t);

4+3*t

>    L(5*t^4-80/3*t^3+253/3*t^2-1544/9*t);

-7*t^2-4658/9+15*t^4

>    L(1);

3

>    got_it:=5*t^4-80/3*t^3+253/3*t^2-1544/9*t+4658/27;
L(got_it);

got_it := 5*t^4-80/3*t^3+253/3*t^2-1544/9*t+4658/27

15*t^4-7*t^2

>   

This step by step procedure should be enjoyable ... but it easy to again make this into a one-step
linear algebra problem.... (but its solution basically involves exactly the same amount of work as

above recursive procedure).

>    try_t4:=A*t^4+B*t^3+C*t^2+D*t+E:
L(try_t4)=rhs(exa2);

12*A*t^2+6*B*t+2*C+16*A*t^3+12*B*t^2+8*C*t+4*D+3*A*t^4+3*B*t^3+3*C*t^2+3*D*t+3*E = 15*t^4-7*t^2

For polynomials it is easy to make MAPLE colelct terms and equate coefficients to zero ... but many

novice users may prefer to read of the equati0on by hand, and re-eneter the systems to be solved
manually.  
The following tries to use MAPLE to extract these equations:

>    poleq:=collect(L(A*t^4+B*t^3+C*t^2+D*t+E)-rhs(exa2),t)=0:
poleq;

(3*A-15)*t^4+(16*A+3*B)*t^3+(12*A+12*B+3*C+7)*t^2+(6*B+8*C+3*D)*t+2*C+3*E+4*D = 0

>    algsys:={seq(coeff(lhs(poleq),t^k)=0,k=1..4),subs(t=0,poleq)}:
 for display purposes only:

>    array(map(q->[q],convert(algsys,list)));

matrix([[6*B+8*C+3*D = 0], [12*A+12*B+3*C+7 = 0], [16*A+3*B = 0], [3*A-15 = 0], [2*C+3*E+4*D = 0]])

>    algsol:=solve(algsys,{A,B,C,D,E}):
algsol;

{A = 5, B = -80/3, D = -1544/9, C = 253/3, E = 4658/27}

>    yp(t):=subs(algsol,try_t4);

yp(t) := 5*t^4-80/3*t^3+253/3*t^2-1544/9*t+4658/27

In this case it is immediate to read of the coefficients which yield a solution. The general
solution then is

>    sol2:=y(t)=rhs(homsol)+yp(t):
sol2;

y(t) = _C1*exp(-3*t)+_C2*exp(-t)+5*t^4-80/3*t^3+253/3*t^2-1544/9*t+4658/27

check it!

>    L(rhs(sol2))-rhs(exa2);;

0

>   

The third example: sinusoids

Very typical for applications are (much longer) expressions of the following form

>    f9:=1*sin(1*t)+1/3*sin(3*t)+1/5*sin(5*t)+1/7*sin(7*t)+1/9*sin(9*t):
exa3:=L(y(t))=f9:
exa3;

diff(y(t),`$`(t,2))+4*diff(y(t),t)+3*y(t) = sin(t)+1/3*sin(3*t)+1/5*sin(5*t)+1/7*sin(7*t)+1/9*sin(9*t)

These arise typically as Fourier expansions of naturally occurring "signals" -- we may get a little
more of this later. Just to get some idea look at the graph of this "9-th order Fourier approximation"

of a rectangular wave (just like an on-off digital signal, think digital TV or CD, or DVD, ...).

>    plot(f9,t=-5..20,-3..3,numpoints=500);

[Maple Plot]

You can imagine what will happen if we take 100 more terms.... (overlaid in blue)

>    plot([sum(1/(2*k-1)*sin((2*k-1)*t),k=1..50),f9],t=-5..20,
     -3..3,color=[blue,red],thickness=[2,1],numpoints=500);

[Maple Plot]

After all the previous work, we expect that a promising Ansatz should contain at leats a linear
combination of the terms of right hand side -- but for this class, in this exploration let us just look
at the first five terms....

>    L(b[1]*sin(1*t)+b[3]*sin(3*t)+b[5]*sin(5*t)+b[7]*sin(7*t)+b[9]*sin(9*t));

2*b[1]*sin(t)-6*b[3]*sin(3*t)-22*b[5]*sin(5*t)-46*b[7]*sin(7*t)-78*b[9]*sin(9*t)+4*b[1]*cos(t)+12*b[3]*cos(3*t)+20*b[5]*cos(5*t)+28*b[7]*cos(7*t)+36*b[9]*cos(9*t)

Not a real big surpise, the output also includes many cosine terms. Thus let us improve our
Ansatz to:

>    try3:=a[1]*cos(1*t)+a[3]*cos(3*t)+a[5]*cos(5*t)+a[7]*cos(7*t)+a[9]*cos(9*t)+
      b[1]*sin(1*t)+b[3]*sin(3*t)+b[5]*sin(5*t)+b[7]*sin(7*t)+b[9]*sin(9*t);

try3 := a[1]*cos(t)+a[3]*cos(3*t)+a[5]*cos(5*t)+a[7]*cos(7*t)+a[9]*cos(9*t)+b[1]*sin(t)+b[3]*sin(3*t)+b[5]*sin(5*t)+b[7]*sin(7*t)+b[9]*sin(9*t)

Evaluating L on this input and equating it to the given right hand side yields

>    L(try3)=f9;

12*b[3]*cos(3*t)-22*a[5]*cos(5*t)+28*b[7]*cos(7*t)+2*a[1]*cos(t)-36*a[9]*sin(9*t)+4*b[1]*cos(t)+20*b[5]*cos(5*t)-6*a[3]*cos(3*t)-12*a[3]*sin(3*t)-28*a[7]*sin(7*t)-46*a[7]*cos(7*t)-78*a[9]*cos(9*t)-78*b...
12*b[3]*cos(3*t)-22*a[5]*cos(5*t)+28*b[7]*cos(7*t)+2*a[1]*cos(t)-36*a[9]*sin(9*t)+4*b[1]*cos(t)+20*b[5]*cos(5*t)-6*a[3]*cos(3*t)-12*a[3]*sin(3*t)-28*a[7]*sin(7*t)-46*a[7]*cos(7*t)-78*a[9]*cos(9*t)-78*b...
12*b[3]*cos(3*t)-22*a[5]*cos(5*t)+28*b[7]*cos(7*t)+2*a[1]*cos(t)-36*a[9]*sin(9*t)+4*b[1]*cos(t)+20*b[5]*cos(5*t)-6*a[3]*cos(3*t)-12*a[3]*sin(3*t)-28*a[7]*sin(7*t)-46*a[7]*cos(7*t)-78*a[9]*cos(9*t)-78*b...

or

>    trig_eq:=L(try3)-f9=0:
trig_eq;

12*b[3]*cos(3*t)-22*a[5]*cos(5*t)+28*b[7]*cos(7*t)+2*a[1]*cos(t)-36*a[9]*sin(9*t)+4*b[1]*cos(t)+20*b[5]*cos(5*t)-6*a[3]*cos(3*t)-12*a[3]*sin(3*t)-28*a[7]*sin(7*t)-46*a[7]*cos(7*t)-78*a[9]*cos(9*t)-78*b...
12*b[3]*cos(3*t)-22*a[5]*cos(5*t)+28*b[7]*cos(7*t)+2*a[1]*cos(t)-36*a[9]*sin(9*t)+4*b[1]*cos(t)+20*b[5]*cos(5*t)-6*a[3]*cos(3*t)-12*a[3]*sin(3*t)-28*a[7]*sin(7*t)-46*a[7]*cos(7*t)-78*a[9]*cos(9*t)-78*b...
12*b[3]*cos(3*t)-22*a[5]*cos(5*t)+28*b[7]*cos(7*t)+2*a[1]*cos(t)-36*a[9]*sin(9*t)+4*b[1]*cos(t)+20*b[5]*cos(5*t)-6*a[3]*cos(3*t)-12*a[3]*sin(3*t)-28*a[7]*sin(7*t)-46*a[7]*cos(7*t)-78*a[9]*cos(9*t)-78*b...

One may think that it is easy to just collect "like-terms" and equate coefficients to zero -- but it seems
to require more advanced thinking methods to do this nicely.
The novice user may prefer to do this step BY HAND, and then MANUALLY re-enter the system
of simultaneous linear equations to be solved.
But the author of this worksheet feels compelled to find some way to get these equations directly...
using brute-force computation of Fourier coefficients:

>    algsys:={
  seq(simplify(1/Pi*int(lhs(trig_eq)*sin((2*k-1)*t),t=0..2*Pi))=0,k=1..5),
  seq(simplify(1/Pi*int(lhs(trig_eq)*cos((2*k-1)*t),t=0..2*Pi))=0,k=1..5)};

algsys := {-22*a[5]+20*b[5] = 0, 2*a[1]+4*b[1] = 0, 2*b[1]-1-4*a[1] = 0, -46*b[7]-28*a[7]-1/7 = 0, 36*b[9]-78*a[9] = 0, -20*a[5]-1/5-22*b[5] = 0, 28*b[7]-46*a[7] = 0, 12*b[3]-6*a[3] = 0, -12*a[3]-1/3-6...
algsys := {-22*a[5]+20*b[5] = 0, 2*a[1]+4*b[1] = 0, 2*b[1]-1-4*a[1] = 0, -46*b[7]-28*a[7]-1/7 = 0, 36*b[9]-78*a[9] = 0, -20*a[5]-1/5-22*b[5] = 0, 28*b[7]-46*a[7] = 0, 12*b[3]-6*a[3] = 0, -12*a[3]-1/3-6...

For display purposes only -- we get a system of 10 linear equations in 10 unknowns, which really
is decoupled into five systems of two simultaneous equations in two unknowns each.
Linearity at play again!!!! we could have looked at one sinusoid at a time, and then added all the
solutions together.

>    array(map(q->[q],convert(algsys,list)));

matrix([[-22*a[5]+20*b[5] = 0], [2*a[1]+4*b[1] = 0], [2*b[1]-1-4*a[1] = 0], [-46*b[7]-28*a[7]-1/7 = 0], [36*b[9]-78*a[9] = 0], [-20*a[5]-1/5-22*b[5] = 0], [28*b[7]-46*a[7] = 0], [12*b[3]-6*a[3] = 0], [...

>    myvars:={seq(a[2*k-1],k=1..5),seq(b[2*k-1],k=1..5)};

myvars := {a[9], a[5], a[7], a[3], b[7], b[9], b[5], a[1], b[1], b[3]}

>    algsol:=solve(algsys,myvars):
op(algsol);

b[9] = -13/11070, b[3] = -1/90, a[3] = -1/45, b[1] = 1/10, a[1] = -1/5, a[7] = -1/725, b[7] = -23/10150, a[5] = -1/221, b[5] = -11/2210, a[9] = -1/1845

This yields the particular solution

>    yp:=subs(algsol,try3);

yp := -1/5*cos(t)-1/45*cos(3*t)-1/221*cos(5*t)-1/725*cos(7*t)-1/1845*cos(9*t)+1/10*sin(t)-1/90*sin(3*t)-11/2210*sin(5*t)-23/10150*sin(7*t)-13/11070*sin(9*t)

check it:

>    L(yp);

sin(t)+1/3*sin(3*t)+1/5*sin(5*t)+1/7*sin(7*t)+1/9*sin(9*t)

>    sol3:=y(t)=rhs(homsol)+yp:
sol3;

y(t) = _C1*exp(-3*t)+_C2*exp(-t)-1/5*cos(t)-1/45*cos(3*t)-1/221*cos(5*t)-1/725*cos(7*t)-1/1845*cos(9*t)+1/10*sin(t)-1/90*sin(3*t)-11/2210*sin(5*t)-23/10150*sin(7*t)-13/11070*sin(9*t)

check it!

>    L(rhs(sol3))-rhs(exa3);

0

>   

The fourth example: terms in the kernel of L

It becomes a little more involved when the right hand side of the inhongeneous differential
equaiton contains terms that coincide with solutions of the associated homgeneous differential
equation:

>    g4:=3*exp(-3*t)-5*t^2*exp(-t):
exa4:=L(y(t))=g4:
exa4;

diff(y(t),`$`(t,2))+4*diff(y(t),t)+3*y(t) = 3*exp(-3*t)-5*t^2*exp(-t)

Start w/ trial and error, trying different functions which L may map to terms similar
to those on the right hand side -- the first experiments:

>    L(exp(-3*t));

0

>    L(-t^2*exp(-t));

-2*exp(-t)-4*t*exp(-t)

>    L(t*exp(-t));

2*exp(-t)

>    L(exp(-t));

0

By taking linear combinations of terms of the above form it appears impossible to get match the
right hand side:

>    simplify(L((A*t^2+B*t+C)*exp(-3*t)+D*exp(-t)))=g4;

-2*exp(-3*t)*(-A+2*A*t+B) = 3*exp(-3*t)-5*t^2*exp(-t)

It should not come as a big surprise that a successful Ansatz also includes higher order
powers of t multiplying these functions. A complete theoretical explanation is found under

then name of "annihilators" (-- the idea is that these right hand sides are themselves solutions
of linear constant coefficent hmogeneous differentail equations) -- for our purposes it suffices

to take an open-minded trial-and-error approach:

>    try4:=(A*t^3+B*t^2+C*t)*exp(-t)+(D*t)*exp(-3*t);

try4 := (A*t^3+B*t^2+C*t)*exp(-t)+D*t*exp(-3*t)

which leads to the equation

>    eq4:=L(try4)-g4=0:
eq4;

(6*A*t+2*B)*exp(-t)+2*(3*A*t^2+2*B*t+C)*exp(-t)-2*D*exp(-3*t)-3*exp(-3*t)+5*t^2*exp(-t) = 0

Again, many may prefer to manually extract the coefficients that are to be equated to zero -- here

a brute force effort to use MAPLE to extraxt these.

>    eq0:=coeff(collect(lhs(eq4),exp(-t)),exp(-t))=0:
eq0;
eq1:=coeff(collect(lhs(eq4),exp(-3*t)),exp(-3*t))=0:
eq1;

6*A*t+2*B+6*A*t^2+4*B*t+2*C+5*t^2 = 0

-2*D-3 = 0

>    pre4:=collect(eq0,t):
pre4;

(6*A+5)*t^2+(4*B+6*A)*t+2*C+2*B = 0

>    eqs:=seq(coeff(lhs(pre4),t,k)=0,k=1..2),subs(t=0,pre4):
eqs;

4*B+6*A = 0, 6*A+5 = 0, 2*C+2*B = 0

For display purposes only -- we get a system of 4 linear equations in 4 unknowns

>    array(map(q->[q],[eq1,eqs]));

matrix([[-2*D-3 = 0], [4*B+6*A = 0], [6*A+5 = 0], [2*C+2*B = 0]])

>    algsol4:=solve({eq1,eqs},{A,B,C,D}):
algsol4;;

{D = -3/2, C = -5/4, A = -5/6, B = 5/4}

This yields the particular solution

>    yp:=subs(algsol4,try4);

yp := (-5/6*t^3+5/4*t^2-5/4*t)*exp(-t)-3/2*t*exp(-3*t)

check it:

>    simplify(L(yp));

3*exp(-3*t)-5*t^2*exp(-t)

>    sol4:=y(t)=rhs(homsol)+yp:
sol4;

y(t) = _C1*exp(-3*t)+_C2*exp(-t)+(-5/6*t^3+5/4*t^2-5/4*t)*exp(-t)-3/2*t*exp(-3*t)

check it!

>    simplify(L(rhs(sol4))-rhs(exa4));

0

>   

>   

Exercise 1 Products of exponentials, sinusoids, and polynomials

>    ex1:=L(y(t))=3*t^3*exp(-2*t)*sin(t)-5*t^2*exp(-t):
ex1;

diff(y(t),`$`(t,2))+4*diff(y(t),t)+3*y(t) = 3*t^3*exp(-2*t)*sin(t)-5*t^2*exp(-t)

Getting started:

>    L(t^3*exp(-2*t)*sin(t));

6*t*exp(-2*t)*sin(t)+6*t^2*exp(-2*t)*cos(t)-2*t^3*exp(-2*t)*sin(t)

What other terms do you think you should include in your Ansatz....?

>   

Exercise 2: Higher order differential equations

A crazy effort to create a linear differential operator whose characteristic polynomial
as designated roots ... a more advanced way would use inverse Laplace transforms. But
in any case the following seems to wok -- cook-up your own polynomial, and auto-

matically get a matching differential operator

>    cookup:=expand((r^2+2*r+2)^2*(r+1)^2);

cookup := r^6+6*r^5+17*r^4+28*r^3+28*r^2+16*r+4

>    solve(cookup=0,r);

-1, -1, -1+I, -1-I, -1+I, -1-I

>    exps:=map(q->diff(q,r)/q*r,convert(cookup,list));
cofs:=map(q->subs(r=1,q),convert(cookup,list));

exps := [6, 5, 4, 3, 2, 1, 0]

cofs := [1, 6, 17, 28, 28, 16, 4]

>    L:=u->convert(zip((a,b)->a*'diff'(u,seq(t,i=1..b)),cofs[1..-2],exps[1..-2]),`+`)+cofs[-1]*u:
L(y(t));

diff(y(t),`$`(t,6))+6*diff(y(t),`$`(t,5))+17*diff(y(t),`$`(t,4))+28*diff(y(t),`$`(t,3))+28*diff(y(t),`$`(t,2))+16*diff(y(t),t)+4*y(t)

>    g2:=(3+5*t)*exp(-t)*sin(t)+t^2*exp(-t):
ex2:=L(y(t))=g2:
ex2;

diff(y(t),`$`(t,6))+6*diff(y(t),`$`(t,5))+17*diff(y(t),`$`(t,4))+28*diff(y(t),`$`(t,3))+28*diff(y(t),`$`(t,2))+16*diff(y(t),t)+4*y(t) = (3+5*t)*exp(-t)*sin(t)+t^2*exp(-t)

Can you get the right Ansatz for the method of undetermined coefficients to succeed?

>    eval(L((A+B*t)*exp(-t)*sin(t)))=g2;

0 = (3+5*t)*exp(-t)*sin(t)+t^2*exp(-t)

What other terms do you think you should include in your Ansatz....?

>   

>   

>   

If all fails, you know that dsolve will give the answer in one step:

>    dsolve(ex2,y(t));

y(t) = (9/4*t*exp(-t)+15/8*t^2*exp(-t)-645/64*exp(-t))*cos(t)+(3/8*t^2*exp(-t)+5/24*exp(-t)*t^3-9/2*exp(-t)-55/8*t*exp(-t))*sin(t)+771/128*exp(-t)+1/12*exp(-t)*t^4-2*t^2*exp(-t)+_C1*exp(-t)+_C2*t*exp(-...
y(t) = (9/4*t*exp(-t)+15/8*t^2*exp(-t)-645/64*exp(-t))*cos(t)+(3/8*t^2*exp(-t)+5/24*exp(-t)*t^3-9/2*exp(-t)-55/8*t*exp(-t))*sin(t)+771/128*exp(-t)+1/12*exp(-t)*t^4-2*t^2*exp(-t)+_C1*exp(-t)+_C2*t*exp(-...

>   

>   

>