3 Differential equation using 4th order runge kutta

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??

[Jan]

Simply insert the 3rd equation.

Currently the question is too general to be answered explicitly.

[James Tursa]

Please post your differential equations and the code you have tried and we can help fix it.

[Yinka Olatunji]

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]

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]

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

[Yinka Olatunji]

and i have to change my function heading to

function diffeqs =ode_Q10(t,var)

or

function diffeqs =ode_Q10(T,var)

[Bjorn Gustavsson]

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

[Yinka Olatunji]

ok thank you so much for your response. They've helped me alot.

However, i have one more code. I have an ode that have rh inside

function diffeqs=ode_Q28(Xh,var)
Xh=var(1); 
t=var(2);

Fh = 0.25;
Fw = 1.0185;
No = 4.0217e3;
Vo = 0.185;
Vf = 3.333e-5;
K1 = 0.0793e-7;
K2 = 0.027e-7;

rh (1,1) = (K1.*Fh.*t.*(1-Xh).*(No - Fh.*t.*Xh)- K2.*Fh.*t.*Xh.*(Fw.*t + Fh + Xh))./((Vo + Vf.*t).^2);
diffeqs(2,1) = (Fh.*t)./(rh.*(Vo + Vf.*t) - Fh.*Xh); %dt/dx

end

and the script

% script

range=[1:10];
ICs=[0,0]; %initial cond

[Xhsol,varsol]=ode45(@ode_Q28,range,ICs);

after running the script, the values generated for "varsol" was "NAN".

I'm wondering where my mistake is.

thank you

[Bjorn Gustavsson]

NO. You will have to first read the documentation for ode45 properly, have a look at a couple of examples and then come back with additional questions. Again you write:

function diffeqs=ode_Q28(Xh,var)
Xh=var(1); 
t=var(2);

When you do this you explicitly overwrite the independent variable Xh with the first dependent variable in var. That is wrong.

A 2-D ode-system should look something like this:

function dxdtdydt = ode_Q28(t,xy)
% t will be given automagically during the ode-integration, not necessarily
% in equal steps in sequence.
x = xy(1); % extract the variables if beneficial for reading code
y = xy(2);
k1 = sqrt(2);
k2 = pi/3;
dxdt = -k1*y; 
dydt = -k2*x;

dxdtdydt = [dxdt; dydt];

Do never again overwrite the independent variable with one of the dependent variables.

[Yinka Olatunji]

Thank you so much.

I really appreciate your time

[Bjorn Gustavsson]

My pleasure, and welcome to matlab.

If you're new to matlab you will get faster up to general speed if you look through the on-ramp material - focusing on your current topics (obviously). That way you will get a better "lay of the land" than any spot-help you can get from "answers".

[Yinka Olatunji]

Yes I'm new to matlab. I will love to get the material. Any link to it?

[Bjorn Gustavsson]

Sure: on-ramp_1