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 Δyn+1+yn = 2, It can also be rewritten as: yn+2-yn+1+yn = 2 (Since, Δyn = yn+1-yn)
Order of a Difference Equation:
The Order of a Difference Equation can be found by the Relation:
Example: Consider the Equation 3yn+2-yn+1+5yn = -12, The Order of it is:
Hence, The Order of the above Equation is ‘2’. So, It can be called a 2nd Order Difference Equation specific.
Homogeneous: If the R.H.S of the above equation is zero, then it can be called a 2nd Order Homogeneous Difference Equation in specific.
Steps to Solve a 2nd 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 r1 & r2), then the yn will be:
yn = C1r1n + C2r2n
Step 5:Else if the Determinant found above is Zero (2 Equal Real Root,r1=r2=r), then the yn will be:
yn = (C1n + C2)r1n
Step 6: Else if the Determinant found above is Negative (Complex Roots, r = α ± iβ), then the yn will be:
yn = Pn(C1cos(nθ) + C2sin(nθ)), where P = √(α2+β2) and θ = tan-1(β/α)
If any Initial Conditions are given, like y0 = m and y1 = 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 C1 & C2 found in the equation ‘yn‘.
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 | w3wiki" ) % 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 | w3wiki" ) % 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: