ode45 solver code for solving a system of three coupled equations does not work

1 view (last 30 days)
raha ahmadi on 16 Jun 2020
Hi ,
I have a system of three ordinary equations and want to solve them numerically. I wrote them with two methods but they have not any result .Hoe can I have output and plot the solution of these equations? sorry but I dont know the difference between two methods
I really appreciate if anyone can help
% syms y(t)
% [V] = odeToVectorField(diff(y, 2) == (1 - y^2)*diff(y) - y)
% M = matlabFunction(V,'vars', {'t','Y'})
% sol = ode45(M,[0 20],[2 0]);
% fplot(@(x)deval(sol,x,1), [0, 20])
close all
clear all
clc
T0=300;
mili=1e-3;
P=90*mili;
Pc=20*mili;
R1=24;%k/w
R2=7;
R3=9.4;
micro=1e-6;
lA=300*micro;
wA=1.2*micro;
tA=0.15*micro;%tickness
vA=lA*wA*tA;%volume
aA=lA*wA;%area
%pc structure
lP=250*micro;
wP=1.5*micro;
tP=0.25*micro;%tickness
vP=lP*wP*tP;
aP=lP*wP;%area
%Bragg section
lB=300*micro;
leff=135*micro;
wB=1.5*micro;
tB=0.25*micro;%tickness
vB=lB*wB*tB;
aB=lB*wB;
nm=1e-9;
RHOv=4.825e3;
V=[vA,vP,vB,vA,vP,vB];
Cv=3.124e2;
Cc=2.46e-7;
Cs=5.32e-5;
Cj=Cv*RHOv.*V;
C=mean(Cj);
m11=-1/(R2*Cc);
m12=-m11;
m13=1/(Cc*R1);
m21=1/(R2*Cs);
m22=-(1/R2+1/R3)*(1/Cs);
m23=0;
m31=6/(R2*Cc);
m32=-m31;
m33=-1/(R1*C)-6/(Cc*R1);
U1=Pc/Cc;
U2=T0/(R3*Cs);
U3=P/C-(6*Pc)/Cc;
syms y(t)
dyd(1)= m11*y(1)+ m12*y(2)+ m13*y(3)+U1;
dyd(2)= m21*y(1)+ m22*y(2)+ m23*y(3)+ U2;
dyd(3)= m31*y(1)+ m32*y(2)+ m33*y(3)+ U3;
dydt=matlabFunction(dyd,'vars',{'t','y'})
sol=ode45(dydt,[0 1],[300 300 0]);
fplot(@(t)deval(sol,t,1),[0 1])
%% metho2
function dydt = odefcn(t,y,U1,U2,U3)
T0=300
mili=1e-3;
P=90*mili;
Pc=20*mili;
R1=24;%k/w
R2=7;
R3=9.4;
micro=1e-6;
lA=300*micro;
wA=1.2*micro;
tA=0.15*micro;%tickness
vA=lA*wA*tA;%volume
aA=lA*wA;%area
%pc structure
lP=250*micro;
wP=1.5*micro;
tP=0.25*micro;%tickness
vP=lP*wP*tP;
aP=lP*wP;%area
%Bragg section
lB=300*micro;
leff=135*micro;
wB=1.5*micro;
tB=0.25*micro;%tickness
vB=lB*wB*tB;
aB=lB*wB;
nm=1e-9;
RHOv=4.825e3;
V=[vA,vP,vB,vA,vP,vB];
Cv=3.124e2;
Cc=2.46e-7;
Cs=5.32e-5;
Cj=Cv*RHOv.*V;
C=mean(Cj);
m11=-1/(R2*Cc);
m12=-m11;
m13=1/(Cc*R1);
m21=1/(R2*Cs);
m22=-(1/R2+1/R3)*(1/Cs);
m23=0;
m31=6/(R2*Cc);
m32=-m31;
m33=-1/(R1*C)-6/(Cc*R1);
U1=Pc/Cc;
U2=T0/(R3*Cs);
U3=P/C-(6*Pc)/Cc;
dydt=zeros(3,1);
dydt(1)= m11*y(1)+ m12*y(2)+ m13*y(3)+U1;
dydt(2)= m21*y(1)+ m22*y(2)+ m23*y(3)+ U2;
dydt(3)= m31*y(1)+ m32*y(2)+ m33*y(3)+ U3;
tspan=[0 1];
y0=[300 300 0];
[t,y] = ode45(@(t,y)odefcn(t,y,U1,U2,U3),tspan,y0);
figure
plot(t,y(:,1))
figure
plot(t,y(:,2))
figure
plot(t,y(:,3))
end

Stephan on 16 Jun 2020
Your ode appears to be stiff - therefore i recommend to use ode15s instead of ode45. Also the behaviour of your system can be seen much better with t=[0 0.01] instead of t=[0 1]. Note that there are several constants that are unused in your code. You see that these values are underlined in orange by Matlab inside the editor window.
However - here's a working code:
[t,y] = solveODE;
subplot(3,1,1)
plot(t,y(:,1))
subplot(3,1,2)
plot(t,y(:,2))
subplot(3,1,3)
plot(t,y(:,3))
function [t,y] = solveODE
%% Constants
T0=300;
mili=1e-3;
P=90*mili;
Pc=20*mili;
R1=24;%k/w
R2=7;
R3=9.4;
micro=1e-6;
lA=300*micro;
wA=1.2*micro;
tA=0.15*micro;%tickness
vA=lA*wA*tA;%volume
aA=lA*wA;%area
%pc structure
lP=250*micro;
wP=1.5*micro;
tP=0.25*micro;%tickness
vP=lP*wP*tP;
aP=lP*wP;%area
%Bragg section
lB=300*micro;
leff=135*micro;
wB=1.5*micro;
tB=0.25*micro;%tickness
vB=lB*wB*tB;
aB=lB*wB;
nm=1e-9;
RHOv=4.825e3;
V=[vA,vP,vB,vA,vP,vB];
Cv=3.124e2;
Cc=2.46e-7;
Cs=5.32e-5;
Cj=Cv*RHOv.*V;
C=mean(Cj);
m11=-1/(R2*Cc);
m12=-m11;
m13=1/(Cc*R1);
m21=1/(R2*Cs);
m22=-(1/R2+1/R3)*(1/Cs);
m23=0;
m31=6/(R2*Cc);
m32=-m31;
m33=-1/(R1*C)-6/(Cc*R1);
U1=Pc/Cc;
U2=T0/(R3*Cs);
U3=P/C-(6*Pc)/Cc;
%% solve the system
tspan=[0 0.01];
y0=[300 300 0];
[t,y] = ode15s(@odefcn,tspan,y0);
%% ODE function
function dydt = odefcn(~,y)
dydt=zeros(3,1);
dydt(1)= m11*y(1)+ m12*y(2)+ m13*y(3)+U1;
dydt(2)= m21*y(1)+ m22*y(2)+ m23*y(3)+ U2;
dydt(3)= m31*y(1)+ m32*y(2)+ m33*y(3)+ U3;
end
end
Stephan on 17 Jun 2020
You can also use ode45 - i also tried at the first attempt - there was no error message arising, but when i noticed that it takes a very long time to calculate, i guessed that the problem may be stiff and changed the solver to ode15s.
So ode45 will also work, but will be time consuming and inefficient. However all you have to do is change this line:
[t,y] = ode15s(@odefcn,tspan,y0);
to
[t,y] = ode45(@odefcn,tspan,y0);

raha ahmadi on 17 Jun 2020
Thank you so much for your great help
with Best wishes

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!