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.