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

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

thanks for any advice in advance

% 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

Answer 1

[Stephan]

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

[raha ahmadi]

Dear Stephen

I just wanted to express how much I appreciate your help. Thank you

my code is messy sorry about that

your code works. thnks very much. I need to solve with Runge Kutta algorithem. Is there any way to solve this problem with ode45?

Best regards

[Stephan]

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);

Answer 2

[raha ahmadi]

Thank you so much for your great help

with Best wishes