%% SECTION 4.3 - THE METHOD OF UNDETERMINED COEFFICIENTS %% % 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 %% % 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) %% % 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)) %% % Now I have to find the roots of the characteristic polynomial. *solve* % assumes the right side of the equation is 0. roots = solve(charpoly) %% % 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 %% % 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) %% % 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) %% % 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)) %% % 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) %% % 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)) %% % I need the coefficient of exp(7x). C = coeffs(Y,exp(7*x)) %% % I solve this equation. [a] = solve(C) %% % I evaluate y_2 using this. Y_2 = eval(y_2) %% % I verify that this is a solution. diffop(Y_2,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) %% % 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)) %% % 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)) %% % The coefficient of cos(3x) is the second term in E1 eq1 = E1(2) %% % Now I find the coefficient of sin(3x). E2 = coeffs(Y,sin(3*x)); eq2 = E2(2) %% % 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] %% % I evaluate y_3 using this. Y_3 = eval(y_3) %% % I verify that this solves the equation. diffop(Y_3,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) %% % I plug this into the left side of the equation and subtract 9 sin(7x). Y = diffop(y_4,x) - 9*sin(7*x) %% % Again I use *coeffs* to get the equations. E1 = coeffs(Y,cos(7*x)); eq1 = E1(2) %% E2 = coeffs(Y,sin(7*x)); eq2 = E2(2) %% % I solve the equations. [a,b] = solve(eq1,eq2) %% % The resulting solution of the equation is Y_4 = eval(y_4) %% % I verify that it solves the equation. diffop(Y_4,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) %% 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)) %% % This is a general solution of the equation. I compare it with our % particular solution. simplify(Y - y_p) %% % I get a solution of the homogeneous equation, which is what I should % get.