EXPERIMENT – 01
Aim: Write a program to draw direction fields for the system of differential equations with constant
coefficients.
Concept: A direction field (or slope field) is a graphical representation used to visualize the solutions of a
system of differential equations without needing to solve them analytically.
𝑑𝑥
= 𝑎𝑥 + 𝑏𝑦
𝑑𝑡
𝑑𝑦
= 𝑐𝑥 + 𝑑𝑦
𝑑𝑡
The variables x and y represent the state of the system and the plane containing all possible (x, y)
coordinate pairs is called the phase plane. At any given point (x, y) in this plane the equation tell us the
rate of change (dx/dt, dy/dt).
Procedure:
A = [2 1;
-1 1];
[x, y] = meshgrid(-5:0.5:5);
u=A(1,1)*x+A(1,2)*y;
v=A(2,1)*x+A(2,2)*y;
figure;
quiver(x,y,u,v);
axis tight;
xlabel("x")
ylabel("y")
title('Direction field for linear system')
eigenvalues=eig(A);
disp('Coefficient Matrix A: ')
disp(A);
disp('Eigenvalues of A: ')
disp(eigenvalues);
Output:
EXPERIMENT – 02
Aim: Write a program to find the steady state solution draw its phase plane plot for the system of ODE and
plot the solution of the system.
Concept: For the system
𝑑𝑥
= 𝑎𝑥 + 𝑏𝑦
𝑑𝑡
𝑑𝑦
= 𝑐𝑥 + 𝑑𝑦
𝑑𝑡
Steady state solution is obtained by setting
𝑑𝑥 𝑑𝑦
=0 =0
𝑑𝑡 𝑑𝑡
This gives equilibrium points of the system. Phase plane plots show trajectories around equilibria. Time
domain plots show system evolution.
Procedure:
a=1;
b=2;
c=-3;
d=-1;
steady_state = [0 0];
disp('Steady state solutions (x*,y*): ');
disp(steady_state);
[x,y]=meshgrid(-5:0.5:5 , -5:0.5:5);
dx=a*x+b*y;
dy=c*x+d*y;
figure;
quiver(x,y,dx,dy,'LineWidth',1);
hold on;
plot(0,0,'ro','MarkerSize',8,'LineWidth',2);
xlabel('x')
ylabel('y')
title('Phase Plane plot of the system');
axis equal;
grid on;
legend('vector Field', 'steady State');
tspan = [0 10];
x0 = [2 -1];
odefun = @(t,x) [
a*x(1)+b*x(2);
c*x(1)+d*x(2);
];
[t,x]=ode45(odefun , tspan, x0);
figure;
plot(t,x(:,1),'b',t,x(:,2),'r','LineWidth',2);
xlabel('Time');
ylabel('State Variables ');
title('Solution of the ODE System ');
legend('x(t)','y(t)');
grid on;
Output:
EXPERIMENT – 03
Aim: Write a program to find the critical points, Jacobian matrix, eigenvalues and eigenvectors of a system of
differential equation with constant coefficients.
Concept: Critical points must satisfy
𝑑𝑥 𝑑𝑦
=0 =0
𝑑𝑡 𝑑𝑡
The Jacobian matrix is
𝑑𝑥̇ 𝑑𝑥̇
𝑑𝑥 𝑑𝑦
𝐽=
𝑑𝑦̇ 𝑑𝑦̇
[ 𝑑𝑥 𝑑𝑦 ]
Eigenvalues determine stability. Eigenvectors indicate direction of trajectories.
Procedure:
syms x y
f=x+2*y+1;
g=3*x+4*y;
S=solve([f==0 , g==0],[x,y]);
critical_point = double([S.x;S.y]);
A=jacobian([f;g],[x,y]);
A=double(A);
[V, D] = eig(A);
EigenValues=diag(D);
EigenVectors=V;
disp("-----Results-----");
disp("Critical Points: ");
disp(critical_point);
disp('Jacobian Matrix (Coefficient Matrix): ');
disp(A);
disp('EigenValues: ');
disp(EigenValues);
disp('EigenVectors (Column correspond to EigenValues): ');
disp(EigenVectors);
Output:
EXPERIMENT – 04
Aim: Develop a MATLAB function to check the stability of a linear system using the Routh – Hurwitz Criteria.
Concept: The characteristic equation is
𝒂𝒏 𝒔𝒏 + 𝒂𝒏−𝟏 𝒔𝒏−𝟏 + ⋯ + 𝒂𝟎 = 𝟎
Routh-Hurwitz criterion determines stability without solving roots. A routh array is constructed from
coefficients. If all the elements of the first column are positive, the system is stable.
Procedure:
coeff=[1,1,2,24];
n=length(coeff);
cols=ceil(n/2);
RouthTable=zeros(n,cols);
RouthTable(1,:)=[coeff(1:2:n)];
RouthTable(2,:)=[coeff(2:2:n)];
for i=3:n
for j=1:cols -1
numerator = RouthTable(i-1,1)*RouthTable(i-2 ,j+1)-......
RouthTable(i-2,1)*RouthTable(i-1,j+1);
denominator = RouthTable(i-1,1);
if denominator ==0
denominator = 1e-10;
end
RouthTable(i,j)=numerator/denominator;
end
end
first_col=RouthTable(:,1);
sign_changes = 0;
for k =1:length(first_col)-1
if sign(first_col(k)) ~= sign(first_col(k+1))
sign_changes = sign_changes+1;
end
end
fprintf('------------------------------------------\n');
fprintf('Routh Table: \n');
disp(RouthTable);
if sign_changes == 0 && all(first_col ~= 0)
fprintf('Result: System Is Stable \n');
else
fprintf('Result: System Is Unstable \n');
fprintf('Sign changes: %d \n',sign_changes);
end
fprintf('-------------------------------------------\n')
Output:
EXPERIMENT – 05
Aim: Write a program to check the stability of a linear system using the Lyapunov function.
Concept: A non-linear system is
𝑑𝑥
= 𝑓(𝑥, 𝑦)
𝑑𝑡
𝑑𝑦
= 𝑔(𝑥, 𝑦)
𝑑𝑡
Linearization is done near equilibrium [Link] Jacobian Matrix approximates the system locally. Eigenvalues
of the Jacobian determine local stability.
Procedure:
A = [-2 1; -3 -4];
Q=eye(size(A));
P = lyap(A' , Q);
disp('matrix P Obtained from lyapunov equation: ');
disp(P);
eigenvalues_P = eig(P);
if all(eigenvalues_P>0)
disp('p is Positive difinite');
disp('System is Asymptotically Stable(Lyapunov Stable)');
else
disp('P is NOT Positive definite');
disp('Cannot conclude Stability');
end
disp('EigenValues of P are: ');
disp(eigenvalues_P);
Output:
EXPERIMENT – 06
Aim: Simulate the Lotka-Volterra prey predator and analyze the population dynamics over time.
Concept: The Lotka-Volterra prey predator model is described by the system of non-linear differential
equations.
𝑑𝑥
= 𝑎𝑥 − 𝑏𝑥𝑦
𝑑𝑡
𝑑𝑦
= 𝑐𝑥𝑦 − 𝑑𝑦
𝑑𝑡
Where: x(t) → Prey Population
Y(t)→ Predator Population
a → Prey Growth rate
b → Predation rate
c → Predator growth due to prey
d → Predator death rate
The term ax represents the natural growth of the prey population, while the term bxy represents the reduction
in prey due to predation.
The term cxy represents the increase in predator population due to prey availability, while the term dy
represents the natural decay of the predator population.
The system exhibits oscillatory behavior around an equilibrium point, representing cyclic population dynamics.
Hopf bifurcation occurs when eigenvalues cross the imaginary axis:
lambda= α(r) ± jω
As parameter r varies equilibrium loses stability. A limit cycle is created. The system exhibits oscillatory
behaviour.
Procedure:
a =1.1;
b=0.4;
c=0.1;
d=0.4;
params = [a;b;c;d];
y0 = [10; 5];
tspan = 0:0.1:50;
[t, y] = ode45 (@(t,y) LVeqn(t,y,params), tspan, y0);
figure;
plot(t, y(:,1), 'LineWidth', 1);
hold on;
plot(t, y(:,2), 'LineWidth', 1);
hold off;
legend ('Prey', 'Predator');
set (gca, 'LineWidth', 1, 'FontSize', 10);
title('Prey-Predator Model');
xlabel ('Time');
ylabel('Population');
figure;
plot (y(:,1), y(:,2), 'LineWidth', 1);
xlabel('Prey');
ylabel('Predator');
title('Prey vs Predator');
set (gca, 'LineWidth', 1, 'FontSize', 10);
u = 0:1:20;
v = 0:1:10;
[u,v] = meshgrid (u,v);
V1 = a*u-b*u.*v;
V2 = c*u.*v-d*v;
VF = sqrt (V1.^2 + V2.^2);
hold on;
quiver (u, v, V1./VF, V2./VF, 0.5);
axis([0 20 0 10]);
hold off;
function df = LVeqn(~, var, params)
a = params(1);
b = params(2);
c = params(3);
d = params(4);
x = var(1);
y = var(2);
df = zeros(2,1);
df(1) = a*x - b*x*y;
df(2) = c*x*y - d*y;
end
Output:
EXPERIMENT – 07
Aim: Implement the SIR model for an epidemic and analyze how varying parameters (infection rate , recovery
rate) affect the spread of disease.
Concept: The SIR model is
𝑑𝑆
= −𝛽𝑆𝐼
𝑑𝑡
𝑑𝐼
= 𝛽𝑆𝐼 − 𝛾𝐼
𝑑𝑡
𝑑𝑅
= 𝛾𝐼
𝑑𝑡
Here, 𝛽 is infection rate and 𝛾 is recovery rate. The model explains diseases spread dynamics.
Procedure:
beta = 0.3;
gamma = 0.1;
S0 = 0.99;
I0 = 0.01;
R0 = 0;
x0 = [S0 I0 R0];
tspan = [0 160];
sir = @(t,x)[
-beta*x(1)*x(2);
beta*x(1)*x(2)-gamma*x(2);
gamma*x(2)
];
[t,x] = ode45(sir, tspan ,x0);
figure;
plot(t,x(:,1),'b',t,x(:,2),'r',t,x(:,3),'g','LineWidth',2);
xlabel('Time');
ylabel('Population Fraction');
title('SIR Model Simulation');
legend ('Susceptible', 'Infected', 'Recovered');
grid on;
beta_values = [0.2 0.4 0.6];
gamma = 0.1;
figure;
hold on;
for beta = beta_values
sir = @(t, x)[
-beta*x(1)*x(2);
beta*x(1)*x(2) - gamma*x(2);
gamma*x(2)];
[t, x] = ode45 (sir, tspan, x0);
plot(t, x(:,2), 'LineWidth', 2);
end
xlabel('Time');
ylabel('Infected Population');
title('Effect of Infection Rate (\beta)');
legend('\beta = 0.2', '\beta = 0.4', '\beta = 0.6');
grid on;
beta = 0.4;
gamma_values = [0.05 0.1 0.2];
figure;
hold on;
for gamma = gamma_values
sir = @(t, x) [
-beta*x(1)*x(2);
beta*x(1)*x(2) - gamma*x(2);
gamma*x(2)
];
[t, x] = ode45 (sir, tspan, x0);
plot(t, x(:,2), 'LineWidth', 2);
end
xlabel('Time');
ylabel('Infected Population');
title('Effect of Recovery Rate (\gamma)');
legend ('\gamma = 0.05', '\gamma = 0.1', '\gamma 0.2');
grid on;
Output:
EXPERIMENT – 08
Aim: Simulate and analyze the saddle node bifurcation in a simple dynamical system.
Concept: A standard system showing saddle node bifurcation is
𝑑𝑥
= 𝑟 − 𝑥2
𝑑𝑡
r is the bifurcation parameter
for r < 0: No equilibrium points
for r > 0: Two equilibria (one stable and one unstable)
for r = 0: One equilibrium (Bifurcation point)
Procedure:
r = linspace(-2 ,2 ,400);
x1 = sqrt(r);
x2 = -sqrt(r);
x1(imag(x1) ~= 0) = NaN;
x2(imag(x2) ~= 0) = NaN;
figure;
plot(r,x1,'b','LineWidth',2);
hold on;
plot(r,x2,'r--','LineWidth',2);
plot(0,0,'ko','MarkerSize',8,'LineWidth',2);
xlabel('Bifurcation Parameter r');
ylabel('Equilibrium points x');
title('Saddle-node Bifurcation Diagram');
legend('Stable Equilibrium ','Unstable Equilibrium','Bifurcation point');
grid on;
tspan = [0 10];
x0 = 0.5;
r_values = [-1, 0, 1];
figure ;
hold on;
for r = r_values
f = @(t,x) r-x.^2;
[t,x] = ode45(f, tspan,x0);
plot(t,x,'LineWidth',2);
end
xlabel('Time');
ylabel('State x(t)');
title('Time response for different r values');
legend('r=-1','r=0','r=1');
grid on;
Output:
EXPERIMENT – 09
Aim: Simulate a transcritical bifurcation in a simple system.
Concept: A standard system showing transcritical bifurcation is
𝑑𝑥
= 𝑟𝑥 − 𝑥 2
𝑑𝑡
r is the bifurcation parameter.
Equilibrium points are: x = 0, x = r
At r = 0, stability of equilibria is exchanged.
Procedure:
r = linspace( -2 ,2 , 400);
x1 = zeros(size(r));
x2 = r;
figure ;
plot(r , x1,'r--','LineWidth',2);
hold on;
plot(r,x2,'b','LineWidth',2);
plot(0,0 ,'ko','MarkerSize',8,'LineWidth',2);
xlabel('Bifurcation parameter r');
ylabel('Equilibrium Points x');
title('Transcritical Bifurcation Diagram');
legend('x = 0 (Unstable to Stable)' , 'x = r (Stable to Unstable)',...
'Bifurcation Point');
grid on;
tspan = [0 10];
x0 = 0.5;
r_values = [-1, 0.5 , 1];
figure;
hold on;
for r = r_values
f = @(t,x) r*x - x.^2;
[t,x] = ode45(f,tspan,x0);
plot(t,x,'LineWidth',2);
end
xlabel('Time');
ylabel('State x(t)');
title('Time Response Showing Transcritical Bifurcation');
legend('r = -1','r = 0', 'r = 1');
grid on;
Output: