CODEWITHSUNDEEP
matlab
Ayush Ranjan
Ayush Ranjan

Reputation: 1

Not able to solve ODE using ODE45 solver

This is my original code.

function dydt = Biogearsmcdaniel(t,y)
SP= 0.60;
kPT=3.7*10^4;
thetaP =1.35*10^-4;
kPTB= 3.1;
kPTMT=6.3*10^-3;
kPTNT=6.1*10^-4;
SMT=2.6*10^-2;
Mv=0.3;
kMT=6.43*10^-5;
kMTB=36.0;
SNT=7.0*10^-7;
Nv =1*10^8;
kNTB=36.0;
kNTMT=0.16;
kNT=6.1*10^-2;
SB =4.6*10^-2;
kBPT=26.0;
kBR= 0.14;
kBNT=4.0*10^-8;
kPBNA=5.8;
xPN= 0.5;
kPS =6.9*10^3;
xPS=1.3*10^4;
kMP=1.01;
xMP =-37.5;
kMD=5.0*10^-2;
xMD = 0.75;
xMTNF=0.4;
kM6=0.1;
xM6 =1.0;
xM10 =0.297;
kMR=0.05;
SM= 1.0;
kMA=0.2;
kNP= 33.75;
xNP=56.25;
kND= 0.05;
xND =0.4;
kNTNF= 0.2;
xENOSP =1.015; 
k10MA=0.1;
xNTNF= 2.0; 
kENOS= 4.0;
k10TNF= 1.485;
kN6 =1.5;
kNO3= 0.46;
x10TNF=0.05;
xN6 =1.0;
kNOMN= 2.0; 
k106= 5.1*10^-2;
xN10= 0.2; 
kTNFN =2.97;
x106 =8.0*10^-2;
kNR=0.05;
kTNFM=0.1;
k10R= 0.1;
SN =1.0; 
xTNF10= 7.9*10^-2;
x1012=1.0*10^-2;
kNA=0.5;
xTNF6= 5.9*10^-2;
k10= 0.35;
kINOSN =1.5;
kTNF= 1.4;
S10= 1.0*10^-2;
kINOSM =0.1;
k6N= 0.2;
k12M= 0.303;
kINOSEC= 0.1;
k6M= 3.03;
x1210= 0.2525;
xINOSTNF =0.05;
k6TNF= 1.0;
k12= 5.0*10^-2;
kINOSd= 0.05;
x6TNF= 0.1; 
kD= 0.15;
kINOS6= 2.0;
k6NO= 2.97;
kD6= 0.125;
xINOS6 =0.1;
x6NO= 0.4; 
xD6= 0.85;
xINOS10= 0.1;
x610= 0.1782;
xDNO= 0.5;
xINOSNO= 0.3;
x66= 0.5;
kINOS= 0.101;
k6= 0.7;
kENOSEC= 0.05;
xENOSTNF= 0.4;
k10N=0.1;
S6= 1.0*10^-3;

dydt=zeros(19,1);%DAEs 

function f=HU1(x,n,h)

    f=x^h/1+(x/n)^h;

end

function g=HU2(x,n,h)

    g=x^h/(x^h+n^h);

end

function r=HUD(x,n,h)

    r=1/1+(x/n)^h;

end

R=1;

dydt(1)=(SP/kPT)*y(1)*(1-y(1))-(thetaP*y(1)/1+kPTB*y(1))-(kPTMT*y(2)*y(1))-kPTNT*y(3)*y(1);
dydt(2)=SMT*y(3)*Mv/1+kMTB*y(4)-kMT*y(2);
dydt(3)=SNT*R*Nv/((1+kNTB*y(4))*(1+kNTMT*y(2)))-kNT*y(3);
dydt(4)=(SB/1+kBPT)*y(4)*(1-y(4))-kBR*R*y(4)-kBNT*y(3)*y(4);
dydt(5)=SP*y(5)+thetaP*y(1)/1+kPTB*y(1)-kPS*y(5)/xPS+y(5)-kPBNA*y(9)*HU2(y(5),xPN,2);
dydt(6)=(-(kMP*HU2(y(5),xMP,2)+kMD*HU2(1-y(18),xMD,4)*(HU2(y(14),xMTNF,2)+kM6*HU2(y(15),xM6,2)))*HUD(y(16),xM10,2)*y(6)-kMR*(y(6)-SM));
dydt(7)=(kMP*HU2(y(5),xMP,2)+kMD*HU2(1-y(18),xMD,4))*(HU2(y(14),xMTNF,2)+kM6*HU2(y(15),xM6,2))*HUD(y(16),xM10,2)*y(6)-kMA*y(7);
dydt(8)=-((kNP*HU2(y(5),xNP,1))+kND*HU1(1-y(18),xND,2)+kNTNF*HU1(y(14),xNTNF,1)+kN6*HU1(y(15),xN6,2))*HUD(y(16),xN10,2)*y(8)-kNR*(y(8)-SN);
dydt(9)=(kNP*HU2(y(5),xNP,1))+kND*HU1(1-
y(18),xND,2)+kNTNF*HU1(y(14),xNTNF,1)+kN6*HU1(y(15),xN6,2)*HUD(y(16),xN10,2)*y(8)-kNA*y(9);
dydt(10=(kINOSN*y(9)+kINOSM*y(7)+kINOSEC*HU1(y(14),xINOSTNF,2)+kINOS6*HU1(y(15),xINOS6,2))*HUD(y(16),xINOS10,2)*HUD(y(19),xINOSNO,4)-kINOSd*y(10);
dydt(11)=kINOS*(y(10)-y(11));
dydt(12)=kENOSEC*HUD(y(14),xENOSTNF,1)*HUD(y(5),xENOSP,1)-kENOS*y(12);
dydt(13)=kNO3*(y(19)-y(13));
dydt(14)=(kTNFN*y(9)+kTNFM*y(7))*HUD(y(16),xTNF10,2)*HUD(y(15),xTNF6,3)-kTNF*y(14);
dydt(15)=(k6N*y(9)+y(7))*(k6M+k6TNF*HU2(y(14),x6TNF,2)+k6NO*HU2(y(19),x6NO,2))*HUD(y(16),x610,2)*HUD(y(16),x66,1)-k6*(y(15)-S6);
dydt(16)=(k10N*y(9)+y(7))*(k10MA+k10TNF*HU2(y(14),x10TNF,4)+k106*HU2(y(15),x106,4))*((1-k10R)*HUD(y(17),x1012,4)+k10R)-k10*(y(16)-S10);
dydt(17)=k12M*y(7)*HUD(y(16),x1210,2)-k12*y(17);
dydt(18)=kD*(1-y(18))*(y(18)-0.05)-(y(18)-0.05)*kD6*HU2(y(15),xD6,6)*(1/xDNO^2+y(19)^2);
dydt(19)=y(11)*(1+kNOMN*(y(7)+y(9)))+y(12);

end

The below code is the syntax used by me to write for the solver of the code.

tspan = 0:120;

yo=[1 1000 1000 0 0 1000 0 1000 0 0 0 0 0 1.1 1.02 1.05 1.2 0 0];

[t,y]=ode15s(@(t,y)Biogearsmcdaniel,tspan,yo);

for i=1:19

    figure(i) 

    plot(t,y(:,i))

    hold on 

end

I am getting the following error when I run the above file Not enough input arguments.

Error in Biogearsmcdaniel (line 109) dydt(1)=(SP/kPT)y(1)(1-y(1))-(thetaPy(1)/1+kPTBy(1))-(kPTMT*y(2)y(1))-kPTNTy(3)*y(1);

Error in untitled148>@(t,y)Biogearsmcdaniel (line 3) [t,y]=ode15s(@(t,y)Biogearsmcdaniel,tspan,yo);

Error in odearguments (line 92) f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode15s (line 153) odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in untitled148 (line 3) [t,y]=ode15s(@(t,y)Biogearsmcdaniel,tspan,yo);

Please help me out.

Upvotes: 0

Views: 73

Answers (1)

John Bofarull Guix
John Bofarull Guix

Reputation: 820

this is John BG [email protected]

I patched your start script to run without errors.

However, not sure if the obtained result is what you are looking for

enter image description here

steps done to get your script up and running:

1.- you were not feeding ode15s correctly, on telling the name of the function, one has to define the inputs of such function.

2.- I have put each function in a different file

3.- I have made global the large amount parameters, it may not be necessary outside the main support function, but you can take if from here

4.- here's the code of the main support function that doesn't crash

function dydt = Biogearsmcdaniel(t,y)
global SP
global kPT
global thetaP
global kPTB
global kPTMT
global kPTNT
global SMT
global Mv
global kMT
global kMTB
global SNT
global Nv
global kNTB
global kNTMT
global kNT
global SB
global kBPT
global kBR
global kBNT
global kPBNA
global xPN
global kPS
global xPS
global kMP
global xMP
global kMD
global xMD
global xMTNF
global kM6
global xM6
global xM10
global kMR
global SM
global kMA
global kNP
global xNP
global kND
global xND
global kNTNF
global xENOSP
global k10MA
global xNTNF
global kENOS
global k10TNF
global kN6global
global kNO3
global x10TNF
global xN6global
global kNOMN
global k106
global xN10
global kTNFNglobal
global x106global 
global kNR
global kTNFM
global k10R
global SN
global xTNF10
global x1012
global kNA
global xTNF6
global k10
global kINOSN
global kTNF
global S10
global kINOSM
global k6N
global k12M
global kINOSEC
global k6M
global x1210
global xINOSTNF
global k6TNF
global k12
global kINOSd
global x6TNF
global kD
global kINOS6
global k6NO
global kD6
global xINOS6
global x6NO
global xD6
global xINOS10
global x610
global xDNO
global xINOSNO
global x66
global kINOS
global k6
global kENOSEC
global xENOSTNF
global k10N
global S6

SP= 0.60;
kPT=3.7*10^4;
thetaP =1.35*10^-4;
kPTB= 3.1;
kPTMT=6.3*10^-3;
kPTNT=6.1*10^-4;
SMT=2.6*10^-2;
Mv=0.3;
kMT=6.43*10^-5;
kMTB=36.0;
SNT=7.0*10^-7;
Nv =1*10^8;
kNTB=36.0;
kNTMT=0.16;
kNT=6.1*10^-2;
SB =4.6*10^-2;
kBPT=26.0;
kBR= 0.14;
kBNT=4.0*10^-8;
kPBNA=5.8;
xPN= 0.5;
kPS =6.9*10^3;
xPS=1.3*10^4;
kMP=1.01;
xMP =-37.5;
kMD=5.0*10^-2;
xMD = 0.75;
xMTNF=0.4;
kM6=0.1;
xM6 =1.0;
xM10 =0.297;
kMR=0.05;
SM= 1.0;
kMA=0.2;
kNP= 33.75;
xNP=56.25;
kND= 0.05;
xND =0.4;
kNTNF= 0.2;
xENOSP =1.015; 
k10MA=0.1;
xNTNF= 2.0; 
kENOS= 4.0;
k10TNF= 1.485;
kN6 =1.5;
kNO3= 0.46;
x10TNF=0.05;
xN6 =1.0;
kNOMN= 2.0; 
k106= 5.1*10^-2;
xN10= 0.2; 
kTNFN =2.97;
x106 =8.0*10^-2;
kNR=0.05;
kTNFM=0.1;
k10R= 0.1;
SN =1.0; 
xTNF10= 7.9*10^-2;
x1012=1.0*10^-2;
kNA=0.5;
xTNF6= 5.9*10^-2;
k10= 0.35;
kINOSN =1.5;
kTNF= 1.4;
S10= 1.0*10^-2;
kINOSM =0.1;
k6N= 0.2;
k12M= 0.303;
kINOSEC= 0.1;
k6M= 3.03;
x1210= 0.2525;
xINOSTNF =0.05;
k6TNF= 1.0;
k12= 5.0*10^-2;
kINOSd= 0.05;
x6TNF= 0.1; 
kD= 0.15;
kINOS6= 2.0;
k6NO= 2.97;
kD6= 0.125;
xINOS6 =0.1;
x6NO= 0.4; 
xD6= 0.85;
xINOS10= 0.1;
x610= 0.1782;
xDNO= 0.5;
xINOSNO= 0.3;
x66= 0.5;
kINOS= 0.101;
k6= 0.7;
kENOSEC= 0.05;
xENOSTNF= 0.4;
k10N=0.1;
S6= 1.0*10^-3;


dydt=zeros(19,1);%DAEs 

% function f=HU1(x,n,h)
% these functions in separate files
%     f=x^h/1+(x/n)^h;
% 
% end

% function g=HU2(x,n,h)
% 
%     g=x^h/(x^h+n^h);
% 
% end

% function r=HUD(x,n,h)
% 
%     r=1/1+(x/n)^h;
% 
% end

R=1;

dydt(1)=(SP/kPT)*y(1)*(1-y(1))-(thetaP*y(1)/1+kPTB*y(1))-(kPTMT*y(2)*y(1))-kPTNT*y(3)*y(1);
dydt(2)=SMT*y(3)*Mv/1+kMTB*y(4)-kMT*y(2);
dydt(3)=SNT*R*Nv/((1+kNTB*y(4))*(1+kNTMT*y(2)))-kNT*y(3);
dydt(4)=(SB/1+kBPT)*y(4)*(1-y(4))-kBR*R*y(4)-kBNT*y(3)*y(4);
dydt(5)=SP*y(5)+thetaP*y(1)/1+kPTB*y(1)-kPS*y(5)/xPS+y(5)-kPBNA*y(9)*HU2(y(5),xPN,2);
dydt(6)=(-(kMP*HU2(y(5),xMP,2)+kMD*HU2(1-y(18),xMD,4)*(HU2(y(14),xMTNF,2)+kM6*HU2(y(15),xM6,2)))*HUD(y(16),xM10,2)*y(6)-kMR*(y(6)-SM));
dydt(7)=(kMP*HU2(y(5),xMP,2)+kMD*HU2(1-y(18),xMD,4))*(HU2(y(14),xMTNF,2)+kM6*HU2(y(15),xM6,2))*HUD(y(16),xM10,2)*y(6)-kMA*y(7);
dydt(8)=-((kNP*HU2(y(5),xNP,1))+kND*HU1(1-y(18),xND,2)+kNTNF*HU1(y(14),xNTNF,1)+kN6*HU1(y(15),xN6,2))*HUD(y(16),xN10,2)*y(8)-kNR*(y(8)-SN);
dydt(9)=(kNP*HU2(y(5),xNP,1))+kND*HU1(1-y(18),xND,2)+kNTNF*HU1(y(14),xNTNF,1)+kN6*HU1(y(15),xN6,2)*HUD(y(16),xN10,2)*y(8)-kNA*y(9);
dydt(10)=(kINOSN*y(9)+kINOSM*y(7)+kINOSEC*HU1(y(14),xINOSTNF,2)+kINOS6*HU1(y(15),xINOS6,2))*HUD(y(16),xINOS10,2)*HUD(y(19),xINOSNO,4)-kINOSd*y(10);
dydt(11)=kINOS*(y(10)-y(11));
dydt(12)=kENOSEC*HUD(y(14),xENOSTNF,1)*HUD(y(5),xENOSP,1)-kENOS*y(12);
dydt(13)=kNO3*(y(19)-y(13));
dydt(14)=(kTNFN*y(9)+kTNFM*y(7))*HUD(y(16),xTNF10,2)*HUD(y(15),xTNF6,3)-kTNF*y(14);
dydt(15)=(k6N*y(9)+y(7))*(k6M+k6TNF*HU2(y(14),x6TNF,2)+k6NO*HU2(y(19),x6NO,2))*HUD(y(16),x610,2)*HUD(y(16),x66,1)-k6*(y(15)-S6);
dydt(16)=(k10N*y(9)+y(7))*(k10MA+k10TNF*HU2(y(14),x10TNF,4)+k106*HU2(y(15),x106,4))*((1-k10R)*HUD(y(17),x1012,4)+k10R)-k10*(y(16)-S10);
dydt(17)=k12M*y(7)*HUD(y(16),x1210,2)-k12*y(17);
dydt(18)=kD*(1-y(18))*(y(18)-0.05)-(y(18)-0.05)*kD6*HU2(y(15),xD6,6)*(1/xDNO^2+y(19)^2);
dydt(19)=y(11)*(1+kNOMN*(y(7)+y(9)))+y(12);

end

5.- and here is the caller

clear all;clc;close all

global SP
global kPT
global thetaP
global kPTB
global kPTMT
global kPTNT
global SMT
global Mv
global kMT
global kMTB
global SNT
global Nv
global kNTB
global kNTMT
global kNT
global SB
global kBPT
global kBR
global kBNT
global kPBNA
global xPN
global kPS
global xPS
global kMP
global xMP
global kMD
global xMD
global xMTNF
global kM6
global xM6
global xM10
global kMR
global SM
global kMA
global kNP
global xNP
global kND
global xND
global kNTNF
global xENOSP
global k10MA
global xNTNF
global kENOS
global k10TNF
global kN6global
global kNO3
global x10TNF
global xN6global
global kNOMN
global k106
global xN10
global kTNFNglobal
global x106global 
global kNR
global kTNFM
global k10R
global SN
global xTNF10
global x1012
global kNA
global xTNF6
global k10
global kINOSN
global kTNF
global S10
global kINOSM
global k6N
global k12M
global kINOSEC
global k6M
global x1210
global xINOSTNF
global k6TNF
global k12
global kINOSd
global x6TNF
global kD
global kINOS6
global k6NO
global kD6
global xINOS6
global x6NO
global xD6
global xINOS10
global x610
global xDNO
global xINOSNO
global x66
global kINOS
global k6
global kENOSEC
global xENOSTNF
global k10N
global S6

SP= 0.60;
kPT=3.7*10^4;
thetaP =1.35*10^-4;
kPTB= 3.1;
kPTMT=6.3*10^-3;
kPTNT=6.1*10^-4;
SMT=2.6*10^-2;
Mv=0.3;
kMT=6.43*10^-5;
kMTB=36.0;
SNT=7.0*10^-7;
Nv =1*10^8;
kNTB=36.0;
kNTMT=0.16;
kNT=6.1*10^-2;
SB =4.6*10^-2;
kBPT=26.0;
kBR= 0.14;
kBNT=4.0*10^-8;
kPBNA=5.8;
xPN= 0.5;
kPS =6.9*10^3;
xPS=1.3*10^4;
kMP=1.01;
xMP =-37.5;
kMD=5.0*10^-2;
xMD = 0.75;
xMTNF=0.4;
kM6=0.1;
xM6 =1.0;
xM10 =0.297;
kMR=0.05;
SM= 1.0;
kMA=0.2;
kNP= 33.75;
xNP=56.25;
kND= 0.05;
xND =0.4;
kNTNF= 0.2;
xENOSP =1.015; 
k10MA=0.1;
xNTNF= 2.0; 
kENOS= 4.0;
k10TNF= 1.485;
kN6 =1.5;
kNO3= 0.46;
x10TNF=0.05;
xN6 =1.0;
kNOMN= 2.0; 
k106= 5.1*10^-2;
xN10= 0.2; 
kTNFN =2.97;
x106 =8.0*10^-2;
kNR=0.05;
kTNFM=0.1;
k10R= 0.1;
SN =1.0; 
xTNF10= 7.9*10^-2;
x1012=1.0*10^-2;
kNA=0.5;
xTNF6= 5.9*10^-2;
k10= 0.35;
kINOSN =1.5;
kTNF= 1.4;
S10= 1.0*10^-2;
kINOSM =0.1;
k6N= 0.2;
k12M= 0.303;
kINOSEC= 0.1;
k6M= 3.03;
x1210= 0.2525;
xINOSTNF =0.05;
k6TNF= 1.0;
k12= 5.0*10^-2;
kINOSd= 0.05;
x6TNF= 0.1; 
kD= 0.15;
kINOS6= 2.0;
k6NO= 2.97;
kD6= 0.125;
xINOS6 =0.1;
x6NO= 0.4; 
xD6= 0.85;
xINOS10= 0.1;
x610= 0.1782;
xDNO= 0.5;
xINOSNO= 0.3;
x66= 0.5;
kINOS= 0.101;
k6= 0.7;
kENOSEC= 0.05;
xENOSTNF= 0.4;
k10N=0.1;
S6= 1.0*10^-3;

tspan = [0:120];

yo=[1 1000 1000 0 0 1000 0 1000 0 0 0 0 0 1.1 1.02 1.05 1.2 0 0];


figure(1)
hold all

for i=1:19

%     figure(i) 

[t,y]=ode15s(@(t,y) Biogearsmcdaniel(t,yo),tspan,yo);

    plot(t,y(:,i))

    hold on 

end

grid on

6.- there's no need to repeat the hold on inside the for loop, just overlap each resulting plot

yes, all repeated, too many lines, but it's good practice to post as soon as the script is working; overworking scripts often go back to errors initially solved.

Hope it helps

John BG [email protected]

Upvotes: 0

Related Questions