Solving 2nd Order Homogeneous Difference Equations in MATLAB

Difference Equations are Mathematical Equations involving the differences between usually successive values of a function (‘y’) of a discrete variable. Here, a discrete variable means that it is defined only for values that differ by some finite amount, generally a constant (Usually ‘1’). A difference equation can also be defined as a relation between the difference of an unknown function at one or more general values how’s the argument. For example equation Δy_n+1+y_n = 2, It can also be rewritten as: y_n+2-y_n+1+y_n = 2 (Since, Δy_n = y_n+1-y_n)

Order of a Difference Equation:

The Order of a Difference Equation can be found by the Relation: 

Example: Consider the Equation 3y_n+2-y_n+1+5y_n = -12, The Order of it is:

Hence, The Order of the above Equation is ‘2’. So, It can be called a 2^nd Order Difference Equation specific.

Homogeneous: If the R.H.S of the above equation is zero, then it can be called a 2^nd Order Homogeneous Difference Equation in specific. 

Steps to Solve a 2^nd Order  Homogeneous Difference Equation:

Step 1: Let the given 2nd Order Difference Equation is:

ayn+2+byn+1+cyn = 0

Step 2:Then, we reduce the above 2nd Order Difference Equation to its Auxiliary Equation(AE) form:

ar2+br+c = 0

Step 3:Then, we find the Determinant of the above Auxiliary Equation(AE) by the Relation:

Det = (b2 − 4ac)

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

yn = C1r1n + C2r2n

Step 5:Else if the Determinant found above is Zero (2 Equal Real Root,r₁=r₂=r), then the y_n will be:

yn = (C1n + C2)r1n

Step 6: Else if the Determinant found above is Negative (Complex Roots, r = α ± iβ), then the y_n will be: 

yn = Pn(C1cos(nθ) + C2sin(nθ)),
where P = √(α2+β2) and θ = tan-1(β/α)

If any Initial Conditions are given, like y₀ = m and y₁ = n, Then we substitute these two Equations in the above Equation and then we find the C1 & C2 by solving the equation we formed earlier by substitution. Then we resubstitute the C₁ & C₂ found in the equation ‘y_n‘.

MATLAB Functions used in the Below Code are:

• disp(‘txt’): This Method displays the message ‘txt’ to the User.
• input(‘txt’): This Method displays the ‘txt’ and waits for the user to input a value and press the Return key.
• solve(eq): This Method solves the ‘eq’ for the variable present in it.
• abs(z): This method returns the complex modulus of z.
• angle(z): This method returns the phase angle in the interval [-π,π] of z.
• subs(y, old, new): This method returns a copy of y, replacing all occurrences of old with new.
• simplify(eq): This Method performs algebraic simplification of ‘eq’.

Example 1:

Matlab

% MATLAB CODE (WITHOUT INITIAL CONDITIONS):
% To clear all variables from the current workspace
clear all

% To clear all the text from the Command Window, resulting in a clear screen
clc
disp("Solving 2nd Order Homogeneous Difference Equations in MATLAB | GeeksforGeeks")

% To Declare them as Variables
syms r c1 c2 n
F=input('Input the coefficients [a,b,c]: ');

% Coefficients of the 2nd Order Difference Equation
a=F(1);b=F(2);c=F(3);

% Auxiliary Equation Formed
eq=a*(r^2)+b*(r)+c;
S=solve(eq);

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

% Determinant of the Auxiliary Equation Formed
D=b^2-4*a*c;
if D>0
    y1=(r1)^n;
    y2=(r2)^n;
    yn=c1*y1+c2*y2;
elseif D==0
    y1=(r1)^n;
    y2=n*((r2)^n);
    yn=c1*y1+c2*y2;
else
    P=abs(r1);
    t=angle(r1);
    y1=((P)^n)*cos(n*t);
    y2=((P)^n)*sin(n*t);
    yn=c1*y1+c2*y2;
end
disp ('The solution of the difference equation yn=')
disp(yn)

Output:

Example 2:

Matlab

% MATLAB CODE (WITH INITIAL CONDITIONS):
% To clear all variables from the current workspace
clear all
clc
disp("Solving 2nd Order Homogeneous Difference Equations in MATLAB | GeeksforGeeks")

% To Declare them as Variables
syms r c1 c2 n
F=input('Input the coefficients [a,b,c]: ');

% Coefficients of the 2nd Order Difference Equation
a=F(1);b=F(2);c=F(3);

% Auxiliary Equation Formed
eq=a*(r^2)+b*(r)+c;
S=solve(eq);

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

% Determinant of the Auxiliary Equation Formed
D=b^2-4*a*c;
if D>0
    y1=(r1)^n;
    y2=(r2)^n;
    yn=c1*y1+c2*y2;
elseif D==0
    y1=(r1)^n;
    y2=n*((r2)^n);
    yn=c1*y1+c2*y2;
else
    P=abs(r1);
    t=angle(r1);
    y1=((P)^n)*cos(n*t);
    y2=((P)^n)*sin(n*t);
    yn=c1*y1+c2*y2;
end
IC=input('Enter the initial conditions in the form [y0,y1]:');

% Initial conditions
y0=IC(1);y1=IC(2);
eq1=subs(yn,n,0)-y0;
eq2=subs(yn,n,1)-y1;

% Finding C1 & C2
[c1,c2]=solve(eq1,eq2);
yn=simplify(subs(yn));
disp ('The solution of the difference equation yn=')
disp(yn)

Output: