Method of Variation of Parameters to Solve 2nd Order Differential Equations in MATLAB

MATLAB can be used to solve numerically second and higher-order ordinary differential equations. In this article, we will see the method of variation of parameters to Solve 2^nd Order Differential Equations in MATLAB.

Step 1: Let the given 2nd Order Differential Equation in terms of ‘x’ is:

Step 2: Then, we reduce it to its Auxiliary Equation(AE) form: _{r}{^2}+ pr +Q = 0

Step 3: Then, we find the Determinant of the above AE by the Relation:

Step 4: If the Determinant found above is Positive (2 Distinct Real roots r₁ & r₂), then the Complementary Function(CF) will be:

Step 5: If the Determinant found above is Zero (1 Unique Real Root ,r₁=r₂=r), then the Complementary Function(CF) will be: 

Step 6: If the Determinant found above is Negative (Complex Roots, r = α ± iβ), then the Complementary Function(CF) will be:

Step 7: In all the Above 3-Cases, the Coefficient of C₁ is termed ‘y₁‘, and the Coefficient of C₂ is termed ‘y₂‘_.

Step 8: Then we find the Wronskian(W) by the Relation: 

Step 9: After finding ‘W’, we find the Particular Integral (PI) by the Relation: 

Step 10: Finally the General Solution(GS) of the 2^nd Order Differential Equation is found by the Relation:  

Example:

Matlab

% MATLAB code for Method of Variation of Parameters
% to Solve 2nd Order Differential
% Equations in MATLAB
clear all
clc
% To Declare them as Variables
syms r c1 c2 x
disp("Method of Variation of Parameters to Solve
2nd Order Differential Equations in MATLAB | GeeksforGeeks")

E=input("Enter the coefficients of the 2nd Order Differential equation");
X=input("Enter the R.H.S of the 2nd order Differential equation");

% Coefficients of the 2nd Order Differential Equations
AE=a*r^2+b*r+c;
a=E(1); b=E(2); c=E(3);
S=solve(AE);

% Roots of Auxiliary Equation (AE)
r1=S(1); r2=S(2);

% Determinant of Auxiliary Equation (AE)
D=b^2-4*a*c;
if D>0
    y1=exp(r1*x);
    y2=exp(r2*x);
    % Complementary Function
    cf=c1*y1+c2*y2
elseif D==0
    y1=exp(r1*x);
    y2=x*exp(r2*x);
    % Complementary Function
    cf=c1*y1+c2*y2
else
    alpha=real(r1);
    beta=imag(r2);
    y1=exp(alpha*x)*cos(beta*x);
    y2=exp(alpha*x)*sin(beta*x);
     % Complementary Function
    cf=c1*y1+c2*y2
end
W=simplify(y1*diff(y2,x)-y2*diff(y1,x));

% Particular Integral
PI=simplify((-y1)*int(y2*X/W) + (y2)*int(y1*X/W))

% General Solution
GS=simplify(cf+PI)

Output:

Input: y'' -6y' + 25 = e2x + sinx + x

Input: y'' -2y' + 3 = x3 + cosx