3 Differential equation using 4th order runge kutta

Question

Yinka Olatunji 2021-10-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1572773-3-differential-equation-using-4th-order-runge-kutta

评论： Bjorn Gustavsson 2021-10-28

Hi, pls how do i solve a three coupled differential equations using 4th order runge kutta? My area is in advanced kinetics and it seams the invoke ode 45 is for two differential eqn.How do I change it to 3??

4 个评论
显示 2更早的评论隐藏 2更早的评论

Yinka Olatunji 2021-10-28

i have 2 and 3 differential equations, the code i use for 2 goes thus

function diffeqs =ode_Q10(t,var)
t = var(1);
Xe = var(2); 
Te=var(3);
tht = 300;
ko = 4.48e6;
E = 62800;
Rg = 8.314;
Tf = 298;
cf = 3;
Hxn = -2.09e8;
De = 1000;
cp = 4190;
diffeqs (1,1) = (-Xe./tht) + ko.*exp(-E./(Rg.*Te)).*(1-Xe); %dxe/dt
diffeqs(2,1) = ((Tf-Te)./tht) - cf.*(Hxn./De.*cp).*ko.*exp(-E./(Rg.*Te)).*(1-Xe); %dTe/dt
end

and i have the script i run as

% script
range=[1:20];
ICs=[0, 0,298]; %initial cond
[tsol,varsol]=ode45(@ode_Q10,range,ICs);

but unfortunately, my Te values was just 0 apart from Te(1)

I dont know where i made a mistake for it to give me zero values for Te

Yinka Olatunji 2021-10-28

Secondly, for the 3 coupled ode, i have this

function diffeqs=ode_sys(z,var)
z=var(1);
x1=var(2);
x2=var(3);
Te=1302;%degree Rankine
diffeqs(1,1) = (1201.6592 +(96705.4341.*((0.8-x1-x2).*(0.2-x1-x2))./(1-x2).^2)+((2359.5424.*(0.8-x1-x2).*(0.2-x1-x2))./(1-x2).^2))./7.4; %dt/dz
diffeqs(2,1) = 0.0287*(824000.*exp(-13737.37/Te).*(((0.8-x1-x2).*(0.2-x1-x2))./(1-x2).^2)); %dx1/dz
diffeqs(3,1) = 0.0287.*(46.8.*exp(-3464.65/Te).*(((0.8-x1-x2).*(0.2-x1-x2))./(1-x2).^2)); %dx2/dz
end

and my script goes thus,

% script
range=[1:11];
ICs=[4,0,0]; %initial cond
[zsol,varsol]=ode45(@ode_sys,range,ICs);

the values I obtained for z are too large and seems incorrect.Can anyone spot any mistake in my eqns?

my "varsol" result are in 3 column. How do i compute the command of "fprint" to put them in a table?

I willreally appreciate any help.

请先登录，再进行评论。

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

Answer 1

Bjorn Gustavsson 2021-10-28

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1572773-3-differential-equation-using-4th-order-runge-kutta#answer_818728

编辑：Bjorn Gustavsson 2021-10-28

You have a function for a 2-D system of ODEs. That then is messed up:

function diffeqs =ode_Q10(t,var)
% t = var(1); % This overwrites the value of the independent input-variable t. 
            % does this mean that you have some other differential equation
            % in a variable "T" (just to differentiate if from the
            % time-variable "t")? If so then lets introduce it here:
T = var(1);
Xe = var(2); 
Te=var(3);
tht = 300;
ko = 4.48e6;
E = 62800;
Rg = 8.314;
Tf = 298;
cf = 3;
Hxn = -2.09e8;
De = 1000;
cp = 4190;
diffeqs(1,1) = sin(t)*T; % Here you better insert the ODE for the "T" variable.
diffeqs (2,1) = (-Xe./tht) + ko.*exp(-E./(Rg.*Te)).*(1-Xe); %dxe/dt
diffeqs(3,1) = ((Tf-Te)./tht) - cf.*(Hxn./De.*cp).*ko.*exp(-E./(Rg.*Te)).*(1-Xe); %dTe/dt
end

This way you will now have a function for 3 ODEs for dT/dt, dXe/dt and dTe/dt. This would then be integrateable with ode45.

HTH

8 个评论
显示 6更早的评论隐藏 6更早的评论

Bjorn Gustavsson 2021-10-28

The function-head doesn't need to change as far as I can tell, it should be the first version. The chage made is to return the 3 derivatives of "T", Xe, and Te, in the diffeqs array that now is 3-by-1, the function still takes 2 input-arguments: t - the time and the array var for "T", Xe and Te.

Now I see that you've written yet another ode-function. Just a QD-comment on that one:

function diffeqs=ode_sys(z,var)
% z=var(1); % Do not overwrite the independent variable with one of the dependent variables! 
t = var(1); % As per comment below 
x1=var(2);
x2=var(3);
Te=1302;%degree Rankine
diffeqs(1,1) = (1201.6592 +(96705.4341.*((0.8-x1-x2).*(0.2-x1-x2))./(1-x2).^2)+((2359.5424.*(0.8-x1-x2).*(0.2-x1-x2))./(1-x2).^2))./7.4; %dt/dz
diffeqs(2,1) = 0.0287*(824000.*exp(-13737.37/Te).*(((0.8-x1-x2).*(0.2-x1-x2))./(1-x2).^2)); %dx1/dz
diffeqs(3,1) = 0.0287.*(46.8.*exp(-3464.65/Te).*(((0.8-x1-x2).*(0.2-x1-x2))./(1-x2).^2)); %dx2/dz
end

Here you seem to want to solve 3 coupled ODEs with respect to z: dt/dz, dx1/dz and dx2/dz. If you overwrite the independent variable z with the values of t that ode45 will try I don't have energy to untangle what the function ends up calculating - but it is not doing what you want it to...

HTH