# SECTION 4.3 - THE METHOD OF UNDETERMINED COEFFICIENTS

## Contents

Here I solve the inhomogeneous fifth order ODE

## 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

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.