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

1 view (last 30 days)
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

Accepted Answer

Stephan
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
  2 Comments
raha ahmadi
raha ahmadi on 17 Jun 2020
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
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);

Sign in to comment.

More Answers (1)

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

Tags

Products


Release

R2018b

Community Treasure Hunt

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

Start Hunting!