Finding difficult in solving these coupled seven ODE's.. Initially tried to solve by runge kutta method... I want to solve by through ODE45 solver because solving by rk method requires more lines of execute,

1 次查看（过去 30 天）
ruban ramalingam
ruban ramalingam 2018-11-25
clear;
%Runge kutta for four coupled system of equations%
%depsi_dt=A/(1-omega)^(b*sigma(1-hc)*sinh[(B*sigma*(1-H))]
%dh_dt=(h*sigma)*(1-(H/Hstar)*(A/(1-omega)*sinh[(B*sigma*(1-H))]
%domega_dt=D*A/(1-omega)*sinh[(B*sigma*(1-H))]
%Define function handles
%Define function handles
ftf=@(t,tf,tg1,tg2,tag,tab,tpot,tload) (1.4915-1.077E-10*(tf^4-8.03E+8)-0.0408*(tf-305)+9.7E-11*(tf^4-tg1^4));
ftg1=@(t,tf,tg1,tg2,tag,tab,tpot,tload) (6.54E-3-3.2E-12*(tg1^4-8.03E+8)-1.5E-3*(tg1-305)-4.5E-14*(tf^4-tg1^4)+((1.26E-3*(tg2-tg1)^1.27)/((tg1+tg2)^0.31))+5.88E-12*(tg2^4-tg1^4));
ftg2=@(t,tf,tg1,tg2,tag,tab,tpot,tload) (7.104E-3+((0.295*((tag-tg2)^1.304)/((tag+tg2)^0.304))+2.087E-8*(tab^4-tg2^4)+0.261E-8*(tpot^4-tg2^4)-((1.26E-3*(tg2-tg1)^1.27))/((tg1+tg2)^0.31))-5.88E-12*(tg2^4-tg1^4));
ftag=@(t,tf,tg1,tg2,tag,tab,tpot,tload) (0.06*(tab-tag)^1.25-((1.3169*(tag-tg2)^1.304)/((tag+tg2)^0.304)))-6.37E-3*(tag-tpot)-((0.32*(tag-tpot)^1.166)/((tag+tpot)^0.166))-((1.12E-4*(tab-tag)*(tag-305))/(0.8703*(tab-tag)+0.606));
ftab=@(t,tf,tg1,tg2,tag,tab,tpot,tload) (-2.12E-11*(tab^4-tg2^4)-6.16E-4*(tab-tag)+8.103E-3-4.622E-4*(tab-305));
ftpot=@(t,tf,tg1,tg2,tag,tab,tpot,tload) (-3.536E-10*(tpot^4-tg2^4)+38.26*(tab-tpot)+8.635E-4*(tag-tpot)+((0.0442*(tag-tpot)^1.166)/((tag+tg2)^0.166))-((17.168*(tpot-tload)^1.25)/((tpot+tload)^0.25))+0.1646);
ftload=@(t,tf,tg1,tg2,tag,tab,tpot,tload) ((1.1835*(tpot-tload)^1.25)/((tpot+tload)^0.25));
%Initial conditions
t(1)=0;
tf(1)=305;
tg1(1)=304;
tg2(1)=305;
tag(1)=306;
tab(1)=306;
tpot(1)=304;
tload(1)=304;
%step size
h=10;
s=h/2;
tfinal=40;
N=ceil(tfinal/h);
%update loop
for i=1:N-1
%update time
t(i+1)=t(i)+h;
% Runge Kutta Method
K1tf=h*ftf(t(i) ,tf(i) ,tg1(i) ,tg2(i) ,tag(i) ,tab(i) ,tpot(i) ,tload(i));
K1tg1=h*ftg1(t(i) ,tf(i) ,tg1(i) ,tg2(i) ,tag(i) ,tab(i) ,tpot(i) ,tload(i));
K1tg2=h*ftg2(t(i) ,tf(i) ,tg1(i) ,tg2(i) ,tag(i) ,tab(i) ,tpot(i) ,tload(i));
K1tag=h*ftag(t(i) ,tf(i) ,tg1(i) ,tg2(i) ,tag(i) ,tab(i) ,tpot(i) ,tload(i));
K1tab=h*ftab(t(i) ,tf(i) ,tg1(i) ,tg2(i) ,tag(i) ,tab(i) ,tpot(i) ,tload(i));
K1tpot=h*ftpot(t(i) ,tf(i) ,tg1(i) ,tg2(i) ,tag(i) ,tab(i) ,tpot(i) ,tload(i));
K1tload=h*ftload(t(i) ,tf(i) ,tg1(i) ,tg2(i) ,tag(i) ,tab(i) ,tpot(i) ,tload(i));
K2tf=h*ftf(t(i)+s ,tf(i)+1/2*K1tf ,tg1(i)+1/2*K1tg1 ,tg2(i)+1/2*K1tg2 ,tag(i)+1/2*K1tag ,tab(i)+1/2*K1tab ,tpot(i)+1/2*K1tpot ,tload(i)+1/2*K1tload);
K2tg1=h*ftg1(t(i)+s ,tf(i)+1/2*K1tf ,tg1(i)+1/2*K1tg1 ,tg2(i)+1/2*K1tg2 ,tag(i)+1/2*K1tag ,tab(i)+1/2*K1tab ,tpot(i)+1/2*K1tpot ,tload(i)+1/2*K1tload);
K2tg2=h*ftg2(t(i)+s ,tf(i)+1/2*K1tf ,tg1(i)+1/2*K1tg1 ,tg2(i)+1/2*K1tg2 ,tag(i)+1/2*K1tag ,tab(i)+1/2*K1tab ,tpot(i)+1/2*K1tpot ,tload(i)+1/2*K1tload);
K2tag=h*ftag(t(i)+s ,tf(i)+1/2*K1tf ,tg1(i)+1/2*K1tg1 ,tg2(i)+1/2*K1tg2 ,tag(i)+1/2*K1tag ,tab(i)+1/2*K1tab ,tpot(i)+1/2*K1tpot ,tload(i)+1/2*K1tload);
K2tab=h*ftab(t(i)+s ,tf(i)+1/2*K1tf ,tg1(i)+1/2*K1tg1 ,tg2(i)+1/2*K1tg2 ,tag(i)+1/2*K1tag ,tab(i)+1/2*K1tab ,tpot(i)+1/2*K1tpot ,tload(i)+1/2*K1tload);
K2tpot=h*ftpot(t(i)+s ,tf(i)+1/2*K1tf ,tg1(i)+1/2*K1tg1 ,tg2(i)+1/2*K1tg2 ,tag(i)+1/2*K1tag ,tab(i)+1/2*K1tab ,tpot(i)+1/2*K1tpot ,tload(i)+1/2*K1tload);
K2tload=h*ftload(t(i)+s ,tf(i)+1/2*K1tf ,tg1(i)+1/2*K1tg1 ,tg2(i)+1/2*K1tg2 ,tag(i)+1/2*K1tag ,tab(i)+1/2*K1tab ,tpot(i)+1/2*K1tpot ,tload(i)+1/2*K1tload);
K3tf=h*ftf(t(i)+s ,tf(i)+1/2*K2tf ,tg1(i)+1/2*K2tg1 ,tg2(i)+1/2*K2tg2 ,tag(i)+1/2*K2tag ,tab(i)+1/2*K2tab ,tpot(i)+1/2*K2tpot ,tload(i)+1/2*K2tload);
K3tg1=h*ftg1(t(i)+s ,tf(i)+1/2*K2tf ,tg1(i)+1/2*K2tg1 ,tg2(i)+1/2*K2tg2 ,tag(i)+1/2*K2tag ,tab(i)+1/2*K2tab ,tpot(i)+1/2*K2tpot ,tload(i)+1/2*K2tload);
K3tg2=h*ftg2(t(i)+s ,tf(i)+1/2*K2tf ,tg1(i)+1/2*K2tg1 ,tg2(i)+1/2*K2tg2 ,tag(i)+1/2*K2tag ,tab(i)+1/2*K2tab ,tpot(i)+1/2*K2tpot ,tload(i)+1/2*K2tload);
K3tag=h*ftag(t(i)+s ,tf(i)+1/2*K2tf ,tg1(i)+1/2*K2tg1 ,tg2(i)+1/2*K2tg2 ,tag(i)+1/2*K2tag ,tab(i)+1/2*K2tab ,tpot(i)+1/2*K2tpot ,tload(i)+1/2*K2tload);
K3tab=h*ftab(t(i)+s ,tf(i)+1/2*K2tf ,tg1(i)+1/2*K2tg1 ,tg2(i)+1/2*K2tg2 ,tag(i)+1/2*K2tag ,tab(i)+1/2*K2tab ,tpot(i)+1/2*K2tpot ,tload(i)+1/2*K2tload);
K3tpot=h*ftpot(t(i)+s ,tf(i)+1/2*K2tf ,tg1(i)+1/2*K2tg1 ,tg2(i)+1/2*K2tg2 ,tag(i)+1/2*K2tag ,tab(i)+1/2*K2tab ,tpot(i)+1/2*K2tpot ,tload(i)+1/2*K2tload);
K3tload=h*ftload(t(i)+s ,tf(i)+1/2*K2tf ,tg1(i)+1/2*K2tg1 ,tg2(i)+1/2*K2tg2 ,tag(i)+1/2*K2tag ,tab(i)+1/2*K2tab ,tpot(i)+1/2*K2tpot ,tload(i)+1/2*K2tload);
K4tf=h*ftf(t(i)+h ,tf(i)+h*K3tf ,tg1(i)+h*K3tg1 ,tg2(i)+h*K3tg2 ,tag(i)+h*K3tag ,tab(i)+h*K3tab ,tpot(i)+h*K3tpot ,tload(i)+h*K3tload);
K4tg1=h*ftg1(t(i)+h ,tf(i)+h*K3tf ,tg1(i)+h*K3tg1 ,tg2(i)+h*K3tg2 ,tag(i)+h*K3tag ,tab(i)+h*K3tab ,tpot(i)+h*K3tpot ,tload(i)+h*K3tload);
K4tg2=h*ftg2(t(i)+h ,tf(i)+h*K3tf ,tg1(i)+h*K3tg1 ,tg2(i)+h*K3tg2 ,tag(i)+h*K3tag ,tab(i)+h*K3tab ,tpot(i)+h*K3tpot ,tload(i)+h*K3tload);
K4tag=h*ftag(t(i)+h ,tf(i)+h*K3tf ,tg1(i)+h*K3tg1 ,tg2(i)+h*K3tg2 ,tag(i)+h*K3tag ,tab(i)+h*K3tab ,tpot(i)+h*K3tpot ,tload(i)+h*K3tload);
K4tab=h*ftab(t(i)+h ,tf(i)+h*K3tf ,tg1(i)+h*K3tg1 ,tg2(i)+h*K3tg2 ,tag(i)+h*K3tag ,tab(i)+h*K3tab ,tpot(i)+h*K3tpot ,tload(i)+h*K3tload);
K4tpot=h*ftpot(t(i)+h ,tf(i)+h*K3tf ,tg1(i)+h*K3tg1 ,tg2(i)+h*K3tg2 ,tag(i)+h*K3tag ,tab(i)+h*K3tab ,tpot(i)+h*K3tpot ,tload(i)+h*K3tload);
K4tload=h*ftload(t(i)+h ,tf(i)+h*K3tf ,tg1(i)+h*K3tg1 ,tg2(i)+h*K3tg2 ,tag(i)+h*K3tag ,tab(i)+h*K3tab ,tpot(i)+h*K3tpot ,tload(i)+h*K3tload);
%% Repeating iteration step by step
tf(i+1)=tf(i)+1/6*(K1tf+2*K2tf+2*K3tf+K4tf);
tg1(i+1)=tg1(i)+1/6*(K1tg1+2*K2tg1+2*K3tg1+K4tg1);
tg2(i+1)=tg2(i)+1/6*(K1tg2+2*K2tg2+2*K3tg2+K4tg2);
tag(i+1)=tag(i)+1/6*(K1tag+2*K2tag+2*K3tag+K4tag);
tab(i+1)=tab(i)+1/6*(K1tab+2*K2tab+2*K3tab+K4tab);
tpot(i+1)=tpot(i)+1/6*(K1tpot+2*K2tpot+2*K3tpot+K4tpot);
tload(i+1)=tload(i)+1/6*(K1tload+2*K2tload+2*K3tload+K4tload);
end
%plot the solution
% plot(t,tf)
% hold on
% plot(t,tg1)
% hold on
% plot(t,tg2)
% hold on
% plot(t,tag)
% hold on
% plot(t,tpot)
% hold on
plot(t,tload)
Please help me to sort out this problem... Help will be appreciated.

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

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics and Optimization 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by