hello, I have three complex differential equations (dTb/dt , dTw/dt , dTg/dt) and I solved it using runge kutta method , matlab run successfully but all result on commend window is "NoN",I'm sure from my equations.

3 次查看（过去 30 天）
mohammad abu abbas
评论： mohammad abu abbas 2018-4-6
is rung kutta method solve very complex equation such equation in this file
  3 个评论
显示 1更早的评论
mohammad abu abbas
mohammad abu abbas 2018-4-5
编辑：James Tursa 2018-4-5
clear all , close all, clc
%defined constant
Ta=20
mw=4
mg=1.1
mb=4
absb=0.95
Ab=1
cpb=460
Asides=0.002016
absg=0.05
tawg=0.85
eg=0.88
Ag=1.18
cpg=840
absw=0.05
taww=0.9
ew=0.96
Aw=1
cpw=4190
rohw=1027
kw=0.613
%hfg=2350000
stefan=5.6697*1.0e-08
V=3
%R=8.3144621
I=600
vis=0.0006527 %viscosity
len=1 %charastaristic length
ki=0.28
Li=0.02
B=4.2.*10.^-2
g=9.81
%Pw=((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))
%Pg=((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))
%eeff=(((1./ew +1./eg)-1).^-1)
%Qcbw=((0.54.*(kw./len).*((((g.*B.*(Tb-Tw).*(rohw.^2).*(len).^3)./(vis.^2)).*(vis.*cpw./kw))).^0.25).*Ab.*(Tb-Tw))
%Ra=(((g.*B.*(Tb-Tw).*(rohw.^2).*(len).^3)./(vis.^2)).*(vis.*cpw./kw))
%hcbw=(0.54.*(kw./len).*((((g.*B.*(Tb-Tw).*(rohw.^2).*(len).^3)./(vis.^2)).*(vis.*cpw./kw))).^0.25)
%Ub=(ki./Li)
%Qloss=((ki./Li).*(Ab+Asides).*(Tb-Ta))
%Qcwg=((0.884.*(Tw-Tg+(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))).*(Tw+273.15)./(268900-((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760)))).^(0.333)).*Aw.*(Tw-Tg))
%hcwg=(0.884.*(Tw-Tg+(Pw-Pg).*(Tw+273.15)./(268900-Pw)).^(0.333))
%Qrwg=((((1./ew +1./eg)-1).^-1).*stefan.*Aw.*((Tw+273.15).^4-(Tg+273.15).^4))
%Qewg=(((0.016237.*((0.884.*(Tw-Tg+(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))).*(Tw+273.15)./(268900-((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760)))).^(0.333)).*Aw.*(Tw-Tg)).*(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))))./(Tw-Tg)).*Aw.*(Tw-Tg))
%hewg=((0.016237.*hcwg.*(Pw-Pg))./(Tw-Tg))
%Qrgs=((eg.*stefan.*((Tg+273.15).^4-((Ta-6.0)+273.15).^4./(Tg-Ta))).*Ag.*(Tg-Ta))
%hrgs=(eg.*stefan.*((Tg+273.15).^4-((Ta-6.0)+273.15).^4./(Tg-Ta)))
%Ts=(Ta-6.0)
%Qcga=((5.7+3.8.*V).*Ag.*(Tg-Ta))
%hcga=(5.7+3.8.*V)
%Tb=((I.*Ab.*absb.*tawg.*taww)-Qcbw-Qloss)./(mb.*cpb)
%Tw=((I.*Aw.*absw.*tawg)+Qcbw-Qcwg-Qrwg-Qewg)./(mw.*cpw)
%Tg=((I.*Ag.*absg)+Qcwg+Qrwg+Qewg-Qrgs-Qcga)./(mg.*cpg)
%defind function
fTb=@(t,Tb,Tw,Tg) ((I.*Ab.*absb.*tawg.*taww)-((0.54.*(kw./len).*((((g.*B.*(Tb-Tw).*(rohw.^2).*(len).^3)./(vis.^2)).*(vis.*cpw./kw))).^0.25).*Ab.*(Tb-Tw))-((ki./Li).*(Ab+Asides).*(Tb-Ta)))./(mb.*cpb)
fTw=@(t,Tb,Tw,Tg) ((I.*Aw.*absw.*tawg)+((0.54.*(kw./len).*((((g.*B.*(Tb-Tw).*(rohw.^2).*(len).^3)./(vis.^2)).*(vis.*cpw./kw))).^0.25).*Ab.*(Tb-Tw))-((0.884.*(Tw-Tg+(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))).*(Tw+273.15)./(268900-((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760)))).^(0.333)).*Aw.*(Tw-Tg))-((((1./ew+1./eg)-1).^-1).*stefan.*Aw.*((Tw+273.15).^4-(Tg+273.15).^4))-(((0.016237.*((0.884.*(Tw-Tg+(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))).*(Tw+273.15)./(268900-((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760)))).^(0.333)).*Aw.*(Tw-Tg)).*(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))))./(Tw-Tg)).*Aw.*(Tw-Tg)))./(mw.*cpw)
fTg=@(t,Tb,Tw,Tg) ((I.*Ag.*absg)+((0.884.*(Tw-Tg+(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))).*(Tw+273.15)./(268900-((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760)))).^(0.333)).*Aw.*(Tw-Tg))+((((1./ew +1./eg)-1).^-1).*stefan.*Aw.*((Tw+273.15).^4-(Tg+273.15).^4))+(((0.016237.*((0.884.*(Tw-Tg+(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))).*(Tw+273.15)./(268900-((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760)))).^(0.333)).*Aw.*(Tw-Tg)).*(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))))./(Tw-Tg)).*Aw.*(Tw-Tg))-((eg.*stefan.*((Tg+273.15).^4-((Ta-6.0)+273.15).^4./(Tg-Ta))).*Ag.*(Tg-Ta))-((5.7+3.8.*V).*Ag.*(Tg-Ta)))./(mg.*cpg)
%initial condition
t(1)=0;
Tb(1)=50;
Tw(1)=40;
Tg(1)=25;
%step size
h=60;
tfinal=3600;
N=ceil(tfinal/h);
%update loop
for i=1:N
%update time
t(i+1)=t(i) + h;
%update Tb Tw Tg
k1Tb=fTb(t(i) ,Tb(i) ,Tw(i) ,Tg(i) );
k1Tw=fTw(t(i) ,Tb(i) ,Tw(i) ,Tg(i) );
k1Tg=fTg(t(i) ,Tb(i) ,Tw(i) ,Tg(i) );
k2Tb=fTb(t(i)+h./2,Tb(i)+h./2.*k1Tb,Tw(i)+h./2.*k1Tw,Tg(i)+h./2.*k1Tg);
k2Tw=fTw(t(i)+h./2,Tb(i)+h./2.*k1Tb,Tw(i)+h./2.*k1Tw,Tg(i)+h./2.*k1Tg);
k2Tg=fTg(t(i)+h./2,Tb(i)+h./2.*k1Tb,Tw(i)+h./2.*k1Tw,Tg(i)+h./2.*k1Tg);
k3Tb=fTb(t(i)+h./2,Tb(i)+h./2.*k2Tb,Tw(i)+h./2.*k2Tw,Tg(i)+h./2.*k2Tg);
k3Tw=fTw(t(i)+h./2,Tb(i)+h./2.*k2Tb,Tw(i)+h./2.*k2Tw,Tg(i)+h./2.*k2Tg);
k3Tg=fTg(t(i)+h./2,Tb(i)+h./2.*k2Tb,Tw(i)+h./2.*k2Tw,Tg(i)+h./2.*k2Tg);
k4Tb=fTb(t(i)+h ,Tb(i)+h .*k3Tb,Tw(i)+h .*k3Tw,Tg(i)+h .*k3Tg);
k4Tw=fTw(t(i)+h ,Tb(i)+h .*k3Tb,Tw(i)+h .*k3Tw,Tg(i)+h .*k3Tg);
k4Tg=fTg(t(i)+h ,Tb(i)+h .*k3Tb,Tw(i)+h .*k3Tw,Tg(i)+h .*k3Tg);
Tb(i+1)=Tb(i)+ h./6.*(k1Tb + 2.*k2Tb + 2.*k3Tb + k4Tb);
Tw(i+1)=Tw(i)+ h./6.*(k1Tw + 2.*k2Tw + 2.*k3Tw + k4Tw);
Tg(i+1)=Tg(i)+ h./6.*(k1Tg + 2.*k2Tg + 2.*k3Tg + k4Tg);
end
%plot
plot(t,Tw)

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Abraham Boayue
Abraham Boayue 2018-4-6
Hey Muhammad, Is it possible for you to attach your three complex equations as they were given? The reason why I'm asking is that most people misrepresent mathematical expressions in their matlab code.
  1 个评论
mohammad abu abbas
mohammad abu abbas 2018-4-6
thank you my brother abraham about your comment ,but my question "is matlab solve very very complex differential equations?"when i run my code there is no error appear on commend window.So this mean that no misrepresent mathematical expressions
clear all , close all, clc %defined constant Ta=20 mw=4 mg=1.1 mb=4 absb=0.95 Ab=1 cpb=460 Asides=0.002016 absg=0.05 tawg=0.85 eg=0.88 Ag=1.18 cpg=840 absw=0.05 taww=0.9 ew=0.96 Aw=1 cpw=4190 rohw=1027 kw=0.613 %hfg=2350000 stefan=5.6697*1.0e-08 V=3 %R=8.3144621 I=600 vis=0.0006527 %viscosity len=1 %charastaristic length ki=0.28 Li=0.02 B=4.2.*10.^-2 g=9.81 %Pw=((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760)) %Pg=((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760)) %eeff=(((1./ew +1./eg)-1).^-1) %Qcbw=((0.54.*(kw./len).*((((g.*B.*(Tb-Tw).*(rohw.^2).*(len).^3)./(vis.^2)).*(vis.*cpw./kw))).^0.25).*Ab.*(Tb-Tw)) %Ra=(((g.*B.*(Tb-Tw).*(rohw.^2).*(len).^3)./(vis.^2)).*(vis.*cpw./kw)) %hcbw=(0.54.*(kw./len).*((((g.*B.*(Tb-Tw).*(rohw.^2).*(len).^3)./(vis.^2)).*(vis.*cpw./kw))).^0.25) %Ub=(ki./Li) %Qloss=((ki./Li).*(Ab+Asides).*(Tb-Ta)) %Qcwg=((0.884.*(Tw-Tg+(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))).*(Tw+273.15)./(268900-((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760)))).^(0.333)).*Aw.*(Tw-Tg)) %hcwg=(0.884.*(Tw-Tg+(Pw-Pg).*(Tw+273.15)./(268900-Pw)).^(0.333)) %Qrwg=((((1./ew +1./eg)-1).^-1).*stefan.*Aw.*((Tw+273.15).^4-(Tg+273.15).^4)) %Qewg=(((0.016237.*((0.884.*(Tw-Tg+(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))).*(Tw+273.15)./(268900-((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760)))).^(0.333)).*Aw.*(Tw-Tg)).*(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))))./(Tw-Tg)).*Aw.*(Tw-Tg)) %hewg=((0.016237.*hcwg.*(Pw-Pg))./(Tw-Tg)) %Qrgs=((eg.*stefan.*((Tg+273.15).^4-((Ta-6.0)+273.15).^4./(Tg-Ta))).*Ag.*(Tg-Ta)) %hrgs=(eg.*stefan.*((Tg+273.15).^4-((Ta-6.0)+273.15).^4./(Tg-Ta))) %Ts=(Ta-6.0) %Qcga=((5.7+3.8.*V).*Ag.*(Tg-Ta)) %hcga=(5.7+3.8.*V) %Tb=((I.*Ab.*absb.*tawg.*taww)-Qcbw-Qloss)./(mb.*cpb) %Tw=((I.*Aw.*absw.*tawg)+Qcbw-Qcwg-Qrwg-Qewg)./(mw.*cpw) %Tg=((I.*Ag.*absg)+Qcwg+Qrwg+Qewg-Qrgs-Qcga)./(mg.*cpg) %defind function fTb=@(t,Tb,Tw,Tg) ((I.*Ab.*absb.*tawg.*taww)-((0.54.*(kw./len).*((((g.*B.*(Tb-Tw).*(rohw.^2).*(len).^3)./(vis.^2)).*(vis.*cpw./kw))).^0.25).*Ab.*(Tb-Tw))-((ki./Li).*(Ab+Asides).*(Tb-Ta)))./(mb.*cpb) fTw=@(t,Tb,Tw,Tg) ((I.*Aw.*absw.*tawg)+((0.54.*(kw./len).*((((g.*B.*(Tb-Tw).*(rohw.^2).*(len).^3)./(vis.^2)).*(vis.*cpw./kw))).^0.25).*Ab.*(Tb-Tw))-((0.884.*(Tw-Tg+(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))).*(Tw+273.15)./(268900-((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760)))).^(0.333)).*Aw.*(Tw-Tg))-((((1./ew+1./eg)-1).^-1).*stefan.*Aw.*((Tw+273.15).^4-(Tg+273.15).^4))-(((0.016237.*((0.884.*(Tw-Tg+(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))).*(Tw+273.15)./(268900-((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760)))).^(0.333)).*Aw.*(Tw-Tg)).*(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))))./(Tw-Tg)).*Aw.*(Tw-Tg)))./(mw.*cpw) fTg=@(t,Tb,Tw,Tg) ((I.*Ag.*absg)+((0.884.*(Tw-Tg+(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))).*(Tw+273.15)./(268900-((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760)))).^(0.333)).*Aw.*(Tw-Tg))+((((1./ew +1./eg)-1).^-1).*stefan.*Aw.*((Tw+273.15).^4-(Tg+273.15).^4))+(((0.016237.*((0.884.*(Tw-Tg+(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))).*(Tw+273.15)./(268900-((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760)))).^(0.333)).*Aw.*(Tw-Tg)).*(((10.^(8.0713-(1730.63./(Tw+233.426)))).*(101325./760))-((10.^(8.0713-(1730.63./(Tg+233.426)))).*(101325./760))))./(Tw-Tg)).*Aw.*(Tw-Tg))-((eg.*stefan.*((Tg+273.15).^4-((Ta-6.0)+273.15).^4./(Tg-Ta))).*Ag.*(Tg-Ta))-((5.7+3.8.*V).*Ag.*(Tg-Ta)))./(mg.*cpg) %initial condition t(1)=0; Tb(1)=50; Tw(1)=40; Tg(1)=25; %step size h=60; tfinal=3600; N=ceil(tfinal/h); %update loop for i=1:N %update time t(i+1)=t(i) + h; %update Tb Tw Tg k1Tb=fTb(t(i) ,Tb(i) ,Tw(i) ,Tg(i) ); k1Tw=fTw(t(i) ,Tb(i) ,Tw(i) ,Tg(i) ); k1Tg=fTg(t(i) ,Tb(i) ,Tw(i) ,Tg(i) ); k2Tb=fTb(t(i)+h./2,Tb(i)+h./2.*k1Tb,Tw(i)+h./2.*k1Tw,Tg(i)+h./2.*k1Tg); k2Tw=fTw(t(i)+h./2,Tb(i)+h./2.*k1Tb,Tw(i)+h./2.*k1Tw,Tg(i)+h./2.*k1Tg); k2Tg=fTg(t(i)+h./2,Tb(i)+h./2.*k1Tb,Tw(i)+h./2.*k1Tw,Tg(i)+h./2.*k1Tg); k3Tb=fTb(t(i)+h./2,Tb(i)+h./2.*k2Tb,Tw(i)+h./2.*k2Tw,Tg(i)+h./2.*k2Tg); k3Tw=fTw(t(i)+h./2,Tb(i)+h./2.*k2Tb,Tw(i)+h./2.*k2Tw,Tg(i)+h./2.*k2Tg); k3Tg=fTg(t(i)+h./2,Tb(i)+h./2.*k2Tb,Tw(i)+h./2.*k2Tw,Tg(i)+h./2.*k2Tg); k4Tb=fTb(t(i)+h ,Tb(i)+h .*k3Tb,Tw(i)+h .*k3Tw,Tg(i)+h .*k3Tg); k4Tw=fTw(t(i)+h ,Tb(i)+h .*k3Tb,Tw(i)+h .*k3Tw,Tg(i)+h .*k3Tg); k4Tg=fTg(t(i)+h ,Tb(i)+h .*k3Tb,Tw(i)+h .*k3Tw,Tg(i)+h .*k3Tg); Tb(i+1)=Tb(i)+ h./6.*(k1Tb + 2.*k2Tb + 2.*k3Tb + k4Tb); Tw(i+1)=Tw(i)+ h./6.*(k1Tw + 2.*k2Tw + 2.*k3Tw + k4Tw); Tg(i+1)=Tg(i)+ h./6.*(k1Tg + 2.*k2Tg + 2.*k3Tg + k4Tg); end %plot plot(t,Tw)

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Data Import and Analysis 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by