Using Mathematica to Teach Linear Differential Operators and the Method of Undetermined Coefficients

Itai Seggev
Knox College / Wolfram Research

Joint Mathematics Meetings

January 9, 2011

New Orleans, Louisiana

Problem Statement

Given a linear differential operator P(D) and a polynomial-exponential forcing function f, we wish to find

the general solution to the homogeneous equation P(D)[y]=0

one solution to the driven equation P(D)[y]=f by means of an ansatz

We can then solve any IVP for P(D), but

we need to enter P(D) into Mathematica

Ingredients

Solve[]

CoefficientList[]

ComplexExpand[], Re[], Im[] (if dealing with complex and or trigonometric equations)

Some magic code:

Clear[DPlusS]
DPlusS /: Power[DPlusS[t_, s_:0], n_] := Function[f,Nest[(D[#, t] + s #)&, f, n]]
DPlusS[t_, s_:0]^n_[f___] ^:= ((DPlusS[t, s])^n)[f]
DPlusS[t_,s_:0][f_]:=D[f,t]+s f
Unprotect[Plus, Times];
Plus  /: Plus[ (a_:1) o1_Function , cruft__][f_] := a  o1[f] + Plus[cruft][f]
Plus  /: Plus[ (a_:1)DPlusS[t_, s___] , cruft__][f_] := a  DPlusS[t, s][f] + Plus[cruft][f]
Times /: Times[a_, DPlusS[t_,s___] ][f_] := a DPlusS[t,s][f]
Times /: Times[a_, o1_Function ][f_] := a o1[f]
Protect[Plus, Times];

Playing with DPlusS[t,s] = D+s

The basic operator

DPlusS[t, s][y[t]]

Monomials

Positive powers:

DPlusS[t, 0]^2[y[t]]
DPlusS[t]^2[y[t]]

Positive powers with shift:

DPlusS[t, s]^2[y[t]]//Expand

The identity operator

DPlusS[t, s]^0[y[t]]

A polynomial

Clear[PofD]
PofD[t_,s_ :0] := L DPlusS[t,s]^2 + R DPlusS[t, s] + (((DPlusS[t,s])^0)/C)
PofD[t,0][y[t]]

PofD[t,I ω][y[t]]//Simplify

Step 1: Finding the basic solution set

Simple roots

Assuming a constant coefficient operator, we know the solution is given by the characteristic roots

PofD[t_,s_ :0] := ((DPlusS[t,s])^2)+5DPlusS[t,s]+6((DPlusS[t,s])^0)

The characteristic polynomial is obtain with D->0, i.e., computing P(D)[1]

PofD[t,s][1]
Solve[% == 0, s]
Exp[s t] /. %

So, we have our basic solution set.

Repeated roots

PofD[t_,s_ :0] := ((DPlusS[t,s])^2)+4DPlusS[t,s]+4((DPlusS[t,s])^0)
PofD[t,s][1]
Solve[% == 0, s]
Exp[s t] /. %

Our set is linearly dependent, we need to replace one vector with t (e^(-2t))

Step 2: Find the particular solution (real case)

Clear[PofD]
PofD[t_,s_ :0] := ((DPlusS[t,s])^2)+5DPlusS[t,s]+6((DPlusS[t,s])^0)

A couple of forcing functions

f1[t_] := (t^2) +3
f2[t_] := Exp[7t](t+2)

(Subscript[f, 1])

We know that the particular solution has the same form as the driving polynomial in the case of a simple polynomial driving term.

y1[t_] := a2 (t^2) + a1 t + a0

We must substitute in and solve for the coefficients.  CoefficientList[] to the rescue!

CoefficientList[f1[t],t]
CoefficientList[PofD[t][y1[t]],t]
Solve[%==%%, {a0,a1,a2}]
y1[t] /. First[%]

(Subscript[f, 2])

When there is an exponential term, it is conventional to separate the solution into its polynomial and exponential parts.

Clear[y2,h2]
h2[t_] := a1 t + a0
y2[t_] := h2[t] Exp[7t]

We then equate the polynomial part of the driving terms with the shifted operator acting on the polynomial ansatz:

CoefficientList[f2[t]/Exp[7t],t]
CoefficientList[PofD[t,7][h2[t]],t]
Solve[%==%%, {a0,a1}]
y2[t] /. First[%]

Alternate solution

While it is conventional to do the separation, it is not strictly required.   CoefficientList[] will consider the exponential part of the "coefficient", allowing its cancellation.

CoefficientList[f2[t],t]
CoefficientList[PofD[t][y2[t]],t]
Solve[%==%%, {a0,a1}]
y2[t] /. First[%]

Step 2: Find the particular solution (trig/complex cases)

Sinusoidal driving term

Clear[PofD,f]
PofD[t_,s_ :0] := ((DPlusS[t,s])^2)+5DPlusS[t,s]+6((DPlusS[t,s])^0)
f[t_] := Sin[7t](t+2)

Since it we have a sinusoid, we need to complexify the problem, solving with a driving term (Subscript[f, c])(t)=(t+2)(e^(7 i t)).

Clear[yc,hc]
hc[t_] := a1 t + a0
fc[t_] :=  Exp[7I t](t+2)
yc[t_] := hc[t] Exp[7 I t]

Finding the complex solution is as before.

CoefficientList[fc[t]/Exp[7 I t],t]
CoefficientList[PofD[t, 7I][hc[t]],t]
Solve[%==%%, {a0,a1}]
yc[t_] = yc[t] /. First[%]

We can then find the solution by taking the imaginary part of our complex solution.

y[t_] = ComplexExpand[Im[yc[t]]]//Simplify

Notice that the ComplexExpand[] is essential, since the value above is not valid if t is complex.

Im[yc[t]]//Simplify

For free, we also get the solution for a cosine driving term.

ComplexExpand[Re[yc[t]]]//Simplify

Alternate solution

As before, we are not obligated to separate out the polynomial if we don't want to.

yc[t_] := hc[t] Exp[7 I t]
CoefficientList[fc[t],t]
CoefficientList[PofD[t][yc[t]],t]
Solve[%==%%, {a0,a1}]

Bonus slide: additional cases

Degeneracy

A degeneracy occurs when the characteristic polynomial has a repeated root, and/or the shift equals one of the roots.

Clear[PofD, f]
PofD[t_,s_ :0] := ((DPlusS[t,s])^2)+4DPlusS[t,s]+4((DPlusS[t,s])^0)
f[t_] := (t^2+2t+3)Exp[-2t]

Clear[y,h]
h[t_] := a2 t^2+a1 t + a0
y[t_] := h[t] Exp[-2t]

Now, what happens if we proceed with our normal ansatz without paying attention to the degeneracy?

CoefficientList[f[t]/Exp[-2t],t]
CoefficientList[PofD[t,-2][h[t]],t]
Solve[%==%%, {a0,a1,a2}]

We get no solution: something has gone wrong.  But, there is more information.  The  two lowest terms have dropped out entirely: we have a double degeneracy (repeated root equal to shift)

Hence, we use a polynomial of degree two higher (but possessing the same number of terms).

h[t_] := a4 t^4+a3 t^3 + a2 t^2
CoefficientList[f[t]/Exp[-2t],t]
CoefficientList[PofD[t,-2][h[t]],t]
Solve[%==%%, {a2,a3,a4}]
y[t] /. First[%]

Higher order operators

Clear[PofD, f]
PofD[t_,s_ :0] := ((DPlusS[t,s])^4)+ 4((DPlusS[t,s])^3)+3((DPlusS[t,s])^2)-4DPlusS[t,s]-4((DPlusS[t,s])^0)
f[t_] := (t^2+2t+3)Exp[-t]

Characteristic roots are -2, -2, 1, and -1.

PofD[t,s][1]//Simplify

We have single degenerac: shift equal to a simple root.  So we need a polynomial of one degree higher.  Rest of process is unchanged.

Clear[y,h]
h[t_] := a3 t^3 + a2 t^2+ a1 t
y[t_] := h[t] Exp[-t]
CoefficientList[f[t]/Exp[-t],t]
CoefficientList[PofD[t,-1][h[t]],t]
Solve[%==%%, {a1,a2,a3}]
y[t] /. First[%]

RLC circuit (complete example from lecture notes)

An RLC circuit obeys

L (Overscript[q, ..]) + R (OverDot[q]) + (q/C)==ℰ(t)

If we wanted to express the relationship in terms of current, we would need to differentiate the previous equation to get:

L (Overscript[I, ..]) + R (OverDot[I]) + (I/C)==(OverDot[ℰ])(t).

Getting equation (2) is helpful because expressing charge in terms of current involves a definite integral: q(t)=q((Subscript[t, 0]))+(∫_(Subscript[t, 0])^t)I(s)ds.  If we substitute this equation into  equation (1), we'd get a so-called integro-differential equation.   Differentiating both sides with respect to t gives us a plain old differential equation.

For convenience, we put a bunch of assumptions into the global variable $Assumptions which will shorten many Simplify[] commands.

$Assumptions={a∈Reals, R >0, ω >0, L>0, C>0, t∈Reals}

We define our differential operator

Clear[PofD]
PofD[t_,s_ :0] := L DPlusS[t,s]^2 + R DPlusS[t, s] + ((DPlusS[t,s])^0)/C

We want to solve equation (2) assuming (OverDot[ℰ])(t)=a cos(ωt).  First, we must find the characteristic roots:

Solve[PofD[t,s][1]==0,s]//Simplify

As long as R is non-zero (i.e., there actually is a resistor in the circuit), the characteristic roots have non-zero  real part and hence cannot equal  s=iω.  This eliminates degenerate cases of undetermined coefficients.  We use Mathematica to find the solution, but as in example 3 this is unnecessary.   Since the polynomial multiplying the cosine is a constant, the polynomial h(t) in the driven solution is

h(t)=(a/(P(i ω)))

We ask Mathematica to verify:

Clear[h,y,z]
h[t_] := a0
Solve[CoefficientList[PofD[t,I ω][h[t]],t]=={a}, a0]
z[t_] = h[t] (E^(I ω t ))/. %[[1]]

Since we're focused on a cosine driving term, we need to take the real part of z.  The following would be the steps corresponding to finding the real part by hand.

We make the denominator real by multiplying by its complex conjugate:

((1/C)+I R ω-L (ω^2))((1/C)-I R ω-L (ω^2))//Simplify

We do the same thing to the numerator and use Euler's identity to convert to sines and cosines.

ExpToTrig[a (E^(I t ω)) ((1/C)-I R ω-L (ω^2))]//Expand

We then take combine the numerator and denominator.  Since the denominator is now real, this simply means take the real part of the numerator and keep the denominator unchanged.

yprovisional= FullSimplify[((Re[%])/(%%))]

Notice that we could have gotten to (Subscript[y, provisional]) in one step using ComplexExpand[]:

yprovisional=ComplexExpand[Re[z[t]]]

Also, here is the imaginary part. It looks similar, but the roles of sine and cosine have been reversed:

ComplexExpand[Im[z[t]]]

Let θ=arctan ((1/C)-L (ω^2), -R ω).

φ=ArcTan [(1/C)-L (ω^2), -R ω]

The two-variable arctangent gives an angle between -π and π, taking into account the quadrant of the point.

Then the expression in the Evaluate[] in the cell below is equivalent to (Subscript[y, provisional]). Notice that they have the same denominator, and the equality of the numerators follows from the angle addition formulas and then basic facts cos(arctan(x,y))=x, sin(arctan(x,y))=y.  Just to be sure, we ask Mathematica to verify:

FullSimplify[Evaluate[(( a Cos[ω t+φ])/(((((1/C)-L (ω^2)))^2)+ (((R ω))^2)))//TrigExpand]==yprovisional,TransformationFunctions->{Automatic,(#/.{Cos[ArcTan[d_,n_]]->d, Sin[ArcTan[d_,n_]] ->n})&}]

So, the general solution is one of the following   The only remaining issue is what is the size of   4L versus C (R^2), which affects the homogeneous solution.  When 4L is smaller, we have two decaying exponentials.  When 4L is larger, we have two oscillating, decaying exponentials.  When they are equal, we have an exponential and a polynomial exponential.   This is the electrical equivalent of a critically damped harmonic oscillator, in which the system returns to equilibrium as quickly as possible.  Clearly it is a very interesting case from an engineering standpoint.  So here are the three possible solutions:

yOverDamped[t_] := ((a Cos[ω t+φ])/(((((1/C)-L (ω^2)))^2)+ (((R ω))^2))) + c1 (E^(((-R+(Sqrt[-((4 L)/C)+(R^2)]))/(2 L))t))+ c2 (E^(((-R-(Sqrt[-((4 L)/C)+(R^2)]))/(2 L))t))
yCriticallyDamped[t_] :=  ((a Cos[ω t+φ])/(((((1/C)-L (ω^2)))^2)+ (((R ω))^2))) + c1 (E^((((-C R)/(2 L)))t))+c2 t (E^((((-C R)/(2 L)))t))
yUnderDamped[t_] :=  ((a Cos[ω t+φ])/(((((1/C)-L (ω^2)))^2)+ (((R ω))^2))) + c1 (E^((((-C R)/(2 L)))t))Cos[((Sqrt[-((4 L)/C)+(R^2)])/(2 L))t]+ c2 (E^((((-C R)/(2 L)))t))Sin[((Sqrt[-((4 L)/C)+(R^2)])/(2 L))t]

All in three cases, however, the homogeneous solution decays to zero, leaving just the driven solution asymptotically.   This piece has its largest amplitude when the denominator is minimized, which for fixed R occurs at (1/C)=L (ω^2), or ω=(Sqrt[(1/(LC))]).  Recall that (Sqrt[(1/(LC))]) is what we found for the natural frequency of an LC circuit, so the amplitude is maximized at resonance.  Indeed, if there is no resistance, the amplitude will grow without bound, as we found on problem 3.4.21.

Being finished with our assumptions, we reset them by assigning True to $Assumptions

$Assumptions=True

Is that a parameter which I see varying before me?

Who said we need constant coefficients?

Clear[PofD,y]
PofD[t_] := ((DPlusS[t])^2)-(2/t) DPlusS[t]+(2/(t^2))((DPlusS[t])^0)

Basic solution set

We can construct a basic solution set in the standard manner

DSolve[{PofD[t][y[t]]==0, y[1]==1, y'[1]==0}, y[t],t]
DSolve[{PofD[t][y[t]]==0, y[1]==0, y'[1]==1}, y[t],t]
Wronskian[y[t] /. Join[%,%%], t]

By adding and subtracting the first to/from the second, we would get {t,(t^2)} is a basic solution set

Clear[y1, y2]
y1[t_] := t
y2[t_] := t^2
PofD[t][y1[t]]
PofD[t][y2[t]]
Wronskian[{y1[t],y2[t]},t]

Variation of Parameters

Consider a logarithmic forcing function

f[t_] := ((Log[t])^2)

We can construct the

yGen[t_] = Simplify[Integrate[((-y2[s]f[s])/(Wronskian[{y1[s],y2[s]},s])), {s,1,t}]y1[t] + Integrate[((y1[s]f[s])/(Wronskian[{y1[s],y2[s]},s])), {s,1,t}] y2[t] + c1 y1[t] + c2 y2[t], Assumptions->t>0]//Expand

PofD[t][yGen[t]]//Simplify
yGen[1]
yGen'[1]

Plot[Evaluate[yGen[t] /. {c1->0,c2->0}], {t, 0, 2}]

Use and Reception

Code given as black box to students

advanced CS students can analyze/appreciate it

Students assigned these problems, including higher order equations

Students had little difficulty mastering the use these operators and showed they could self correct if started with wrong polynomial

Emphasize operator/linear algebra view

There is a single theory of linear equations!