matlab 无法使用ODE45求解器求解ODE

h7appiyu  于 2022-11-15  发布在  Matlab
关注(0)|答案(1)|浏览(324)

这是我的原始代码。

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

下面的代码是我用来编写代码解算器的语法。

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

当我运行上面的文件没有足够的输入参数时，我得到以下错误。
Biogearsmcdaniel(第109行)dydt(1)=(SP/kPT)y(1)(1-y(1))-(thetaPy(1)/1+kPTBy(1))-(kPTMT*y(2)y(1))-kPTNTy(3)*y(1)；出错
在Untitled148>@(t，y)Biogearsmcdaniel(第3行)[t，y]=ode15s(@(t，y)Biogearsmcdaniel，tspan，yo)中出错；
Odearguments出错(第92行)f0=ode(t0，y0，args{：})；%ODE15I将args{1}设置为yp0。
Ode15s(第153行)odearguments(odeIsFuncHandle、odeTreatAsMFile、solver_name、ode、tspan、y0、Options、varargin)出错；
未完成148中的错误(第3行)[t，y]=ode15s(@(t，y)Biogearsmcdaniel，tspan，yo)；
请帮帮我。

matlab

1条答案

按热度按时间
c0vxltue

c0vxltue1#

这是John BG jgb2012@sky.com
我修补了您的启动脚本，使其运行时没有错误。
但是，不确定所获得的结果是否就是您要寻找的结果

启动并运行脚本所需的步骤：
1.-您没有正确输入ode15，在告诉函数名称时，必须定义该函数的输入。
2.-我已将每个函数放在不同的文件中
3.-我已经全局设置了大额参数，在主支持函数之外可能不需要，但你可以从这里拿
4.-以下是不会崩溃的主要支持函数的代码

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.-这是打电话的人

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.-不需要在for循环内重复保持，只需重叠每个结果图
是的，所有这些都是重复的，太多行，但在脚本工作后立即发布是一种好做法；过度工作的脚本经常会回到最初解决的错误。
希望能有所帮助
邮箱：John BG jgb2012@sky.com

赞(0）分享  回复(0）举报  2022-11-15
我来回答

相关问题

    查看更多

    热门标签

    更多
    JavaquerypythonNode开发语言requestUtil数据库Table后端算法LoggerMessageElementParser

    最新问答

    更多
    • 蜀ICP备13028337号-1 大数据知识库 https://www.saoniuhuo.com © All rights reserved
    • 本站内容来源互联网,如果侵犯您的权益请联系我们删除, 联系方式：448109455@qq.com