SECTION 4.3 - THE METHOD OF UNDETERMINED COEFFICIENTS

Contents

Here I solve the inhomogeneous fifth order ODE

$$y^v  \,\mbox{--}\,3y^{iv} + 3y''' \,\mbox{--}\, 3y'' + 2y' = x^5 + 6e^{7x}
\,\mbox{--}\, 2e^{\,\mbox{--}\,x}\cos(3x) + 9 \sin(7x).$$

Solving by hand (letting MATLAB do each calculation)

First I have MATLAB carry out the same steps I would do if I were solving it by hand. Then I see what happens if I ask MATLAB to solve the original equation.

Here is a MATLAB function (defined by a function-M file) which will apply the left side of the equation to a function f.

type diffop
function [L] = diffop(f,x)

% DIFFOP applies the differential operator L which takes y to y^v - 3y^{iv} + 3y'''- 3y'' + 2y' to the
% function f(x).  It requires two arguments, a function and an independent variable.

L = diff(f,x,5) - 3*diff(f,x,4) + 3*diff(f,x,3) - 3*diff(f,x,2) + 2*diff(f,x);

The first step is to solve the corresponding homogeneous equation, by plugging

$$y = e^{rx}$$

into the left side.

syms r x
L = diffop(exp(r*x),x)
 
L =
 
2*r*exp(r*x) - 3*r^2*exp(r*x) + 3*r^3*exp(r*x) - 3*r^4*exp(r*x) + r^5*exp(r*x)
 

I divide the left side by exp(rx). I have to give the simplify command to get MATLAB to carry out the division.

charpoly = simplify(L/exp(r*x))
 
charpoly =
 
r*(r^2 + 1)*(r - 1)*(r - 2)
 

Now I have to find the roots of the characteristic polynomial. solve assumes the right side of the equation is 0.

roots = solve(charpoly)
 
roots =
 
  0
  1
  2
  i
 -i
 

I know that 1, exp(x), exp(2x), cos(x) and sin(x) form a fundamental set of solutions of the homogeneous equation.

The right side of the equation has 4 terms, each of a different form, so I break the inhomogeneous problem into 4 problems. For the first, the right side is x^5. So, since 1 solves the homogeneous equation, I look for a solution of the form y_1 equal to x times a polynomial of degree 5. Here is a loop which builds it.

y_1 = 0;
for k = 0:5
vars{k+1} = ['a' num2str(k)];
u = vars{k+1}* x^(6-k);
y_1 = u + y_1;
end
y_1
 
y_1 =
 
a0*x^6 + a1*x^5 + a2*x^4 + a3*x^3 + a4*x^2 + a5*x
 

I plug it into the left side and subtract x^5. The expression is very lengthy. It helps to use collect to collect all the terms involving the same power of x (and factor the power out of the terms) and to use pretty so it doesn't run off the screen.

Y = collect(diffop(y_1,x) - x^5,x); pretty(Y)
               5                    4                            3
  (12 a0 - 1) x  + (10 a1 - 90 a0) x  + (360 a0 - 60 a1 + 8 a2) x  +
  
                                        2
     (180 a1 - 1080 a0 - 36 a2 + 6 a3) x  + (720 a0 - 360 a1 + 72 a2 - 18 a3 + 4 a4) x +
  
     120 a1 - 72 a2 + 18 a3 - 6 a4 + 2 a5

I want to solve the system of equations that I get from setting the coefficents equal to 0. I can either do this by copying and pasting the coefficients into the solve command or using a for loop to calculate the coefficients and set them equal to 0. Here I use a loop to do it. The value of the coefficient of x^j is the jth derivative of Y evaluated at 0. The corresponding equation is indexed by j+1.

for k = 1:6
    eqn = subs(diff(Y,x,k-1),x,0);
    eqns{k} = [char(eqn) '=0'];
end

Now I solve the resulting system.

range = 1:6;
sol = solve(eqns{range},vars{range});
for k = range
    vals{k} = sol.(vars{k});
end

I substitute these into y_1.

Y_1 = subs(y_1,vars,vals)
 
Y_1 =
 
x^6/12 + (3*x^5)/4 + (15*x^4)/8 + (15*x^3)/4 + (285*x^2)/8 + (765*x)/8
 

I verify that this is a solution. To get MATLAB to do the algebra after applying the left side to Y_1, I have to use simplify.

simplify(diffop(Y_1,x))
 
ans =
 
x^5
 

The right side of the second equation is 6 exp(7x), which doesn't solve the homogeneous equation so I look for a solution of the form

syms a

y_2 = a*exp(7*x)
 
y_2 =
 
a*exp(7*x)
 

I plug it into the left side of the equation and subtract 6 exp(7x). This time I want to collect the terms involving exp(7x).

Y = collect(diffop(y_2,x)-6*exp(7*x),exp(7*x))
 
Y =
 
(10500*a - 6)*exp(7*x)
 

I need the coefficient of exp(7x).

C = coeffs(Y,exp(7*x))
 
C =
 
10500*a - 6
 

I solve this equation.

[a] = solve(C)
 
a =
 
1/1750
 

I evaluate y_2 using this.

Y_2 = eval(y_2)
 
Y_2 =
 
exp(7*x)/1750
 

I verify that this is a solution.

diffop(Y_2,x)
 
ans =
 
6*exp(7*x)
 

The third term on the right is - 2exp(-x)cos(3x) which doesn't solve the homoegeneous equation, so I look for a solution of the form

clear a
syms a b
y_3 = a*exp(-x)*cos(3*x)+b*exp(-x)*sin(3*x)
 
y_3 =
 
(a*cos(3*x))/exp(x) + (b*sin(3*x))/exp(x)
 

I plug this into the left side of the equation and subtract - 2exp(-x)cos(3x) When I do this, every term will contain a factor of exp(-x), so I divide the result by that.

Y = simplify((diffop(y_3,x) + 2*exp(-x)*cos(3*x))/exp(-x))
 
Y =
 
2*cos(3*x) + 330*a*sin(3*x) - 300*b*sin(3*x) - 300*a*cos(3*x) - 330*b*cos(3*x)
 

I could get the equations by evaluating Y and its derivative at 0. Instead I will use coeffs command. I notice that Y is a linear "polynomial" in cos(3x) and sin(3x), so I want the coefficients of cos(3x) and of sin(3x).

E1 = coeffs(Y,cos(3*x))
 
E1 =
 
[ 330*a*sin(3*x) - 300*b*sin(3*x), 2 - 330*b - 300*a]
 

The coefficient of cos(3x) is the second term in E1

eq1 = E1(2)
 
eq1 =
 
2 - 330*b - 300*a
 

Now I find the coefficient of sin(3x).

E2 = coeffs(Y,sin(3*x)); eq2 = E2(2)
 
eq2 =
 
330*a - 300*b
 

I solve the resulting equations using solve. The solve command assumes that the right side of each equation is 0.

[a,b] = solve(eq1,eq2); [a,b]
 
ans =
 
[ 2/663, 11/3315]
 

I evaluate y_3 using this.

Y_3 = eval(y_3)
 
Y_3 =
 
(2*cos(3*x))/(663*exp(x)) + (11*sin(3*x))/(3315*exp(x))
 

I verify that this solves the equation.

diffop(Y_3,x)
 
ans =
 
-(2*cos(3*x))/exp(x)
 

The last term on the right is 9 sin(7x). This doesn't solve the homogeneous equation so I look for a solution of the form

clear a b
syms a b
y_4 = a*cos(7*x) + b*sin(7*x)
 
y_4 =
 
b*sin(7*x) + a*cos(7*x)
 

I plug this into the left side of the equation and subtract 9 sin(7x).

Y = diffop(y_4,x) - 9*sin(7*x)
 
Y =
 
15792*b*cos(7*x) - 15792*a*sin(7*x) - 7056*b*sin(7*x) - 7056*a*cos(7*x) - 9*sin(7*x)
 

Again I use coeffs to get the equations.

E1 = coeffs(Y,cos(7*x)); eq1 = E1(2)
 
eq1 =
 
15792*b - 7056*a
 
E2 = coeffs(Y,sin(7*x)); eq2 = E2(2)
 
eq2 =
 
- 15792*a - 7056*b - 9
 

I solve the equations.

[a,b] = solve(eq1,eq2)
 
a =
 
-141/296800
 
 
b =
 
-9/42400
 

The resulting solution of the equation is

Y_4 = eval(y_4)
 
Y_4 =
 
- (141*cos(7*x))/296800 - (9*sin(7*x))/42400
 

I verify that it solves the equation.

diffop(Y_4,x)
 
ans =
 
9*sin(7*x)
 

The resulting particular solution of the original equation is obtained by adding these 4 solutions.

y_p = Y_1 + Y_2 + Y_3 + Y_4; pretty(y_p)
                                                                   2       3       4
  765 x   141 cos(7 x)   exp(7 x)   9 sin(7 x)   2 cos(3 x)   285 x    15 x    15 x
  ----- - ------------ + -------- - ---------- + ---------- + ------ + ----- + ----- +
    8        296800        1750       42400      663 exp(x)     8        4       8
  
        5    6
     3 x    x    11 sin(3 x)
     ---- + -- + -----------
      4     12   3315 exp(x)

Solving with Matlab

What happens if you let MATLAB solve the equation?

eqn =  'D5y - 3*D4y + 3*D3y - 3*D2y + 2*Dy= x^5 + 6*exp(7*x) - 2*exp(-x)*cos(3*x) + 9*sin(7*x)';
Y = dsolve(eqn, 'x');
pretty(simple(Y))
  765 x   C18   141 cos(7 x)   exp(7 x)   9 sin(7 x)
  ----- - --- - ------------ + -------- - ---------- + C19 cos(x) + C22 exp(x) +
    8      2       296800        1750       42400
  
                                    2       3       4      5    6
                  2 cos(3 x)   285 x    15 x    15 x    3 x    x
     C20 sin(x) + ---------- + ------ + ----- + ----- + ---- + -- + C21 exp(2 x) +
                  663 exp(x)     8        4       8      4     12
  
     11 sin(3 x)
     ----------- + 765/16
     3315 exp(x)

This is a general solution of the equation. I compare it with our particular solution.

simplify(Y - y_p)
 
ans =
 
C19*cos(x) - C18/2 + C22*exp(x) + C20*sin(x) + C21*exp(2*x) + 765/16
 

I get a solution of the homogeneous equation, which is what I should get.