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.

is rung kutta method solve very complex equation such equation in this file

[James Tursa]

Please post your code as text in the Question section and format it with the { } Code button.

[mohammad abu abbas]

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)

[ayman alkezza]

Peace, mercy and blessings of allah

My brother Mohammed Abbas I need this code if you are properly running .

Answer 1

[Abraham Boayue]

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.

[mohammad abu abbas]

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)