Solving four ordinary dif equations simultaneously using ode45 matlab

[J-bonz]

I am trying to replicate an example problem that was found in an "into to combustion book"... And of course, yes it is for school. Unfortunately there is no MATLAB class offered but, they expect us to know it... I tried the example in MATLAB help found under the ode45 and couldn't even get this to run. I thought if i could get it to run, I would use this as a model for the problem i am working on... Solving four equations simultanously

Example problem code:
-----------------------------------------------------------
function dy = rigid(t,y)

dy = zeros(3,1);

dy(1) = y(2) * y(3);

dy(2) = -y(1) * y(3);

dy(3) = -0.51 * y(1) * y(2);

options = odeset('RelTol', 1e-4, 'AbsTol', [1e-4 1e-4 1e-5]);

[T,Y] = ode45(@rigid, [0 12], [0 1 1], options);

plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:3),'.');
-----------------------------------------------------------

This Code Shows an error in the "dy(1) = y(2) * y(3);" line.

The equations i really need to solve for are as follows
-----------------------------------------------------------
dF/dt = -4.7316(10^8) * (F^0.1) * (O^1.65) * exp(-15098/T)
F(0) = .023918

dO/dt = 16 * (dF/dt)
O(0) = .382681

dP/dt = -17 * (dF/dt)
P(0) = 0

dT/dt = -3715.19* (dF/dt)
T(0) = 753
-----------------------------------------------------------

Any Help would greatly be appreceated, its due monday...
and it seems everything else is.
Thanks,
Jack

 

[LeBrad]

For the ode solvers in matlab, you need two pieces. If you have the equation
dy/dt = f(t,y)
you need to create a MATLAB function which evalutes f(t,y) when given y and t, and then you need a separate piece which calls ode45, and ode45 takes in as an argument the function name of f(t,y). Your function rigid(t,y) is what evaluates f(t,y), but then you need to pass that as an argument to ode45. So what you have

function dy = rigid(t,y)

dy = zeros(3,1);

dy(1) = y(2) * y(3);

dy(2) = -y(1) * y(3);

dy(3) = -0.51 * y(1) * y(2);

options = odeset('RelTol', 1e-4, 'AbsTol', [1e-4 1e-4 1e-5]);

[T,Y] = ode45(@rigid, [0 12], [0 1 1], options);

plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:3),'.');

should be broken up into two separate things. First you need the function

function dy = rigid(t,y)

dy = zeros(3,1);

dy(1) = y(2) * y(3);

dy(2) = -y(1) * y(3);

dy(3) = -0.51 * y(1) * y(2);

saved as rigid.m, and then you need another file which contains

 options = odeset('RelTol', 1e-4, 'AbsTol', [1e-4 1e-4 1e-5]);

[T,Y] = ode45(@rigid, [0 12], [0 1 1], options);

plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.');

Run this second file and it will declare the options you want, it will call matlab's ode solver which in turn will call your function rigid, and then it will plot the results. Also, you have a type in your original thing in the last line plot(...), you need a comma between the colon and three.

 

[J-bonz]

Thanks, I will give that a try LeBrad

 

[J-bonz]

I was able to get the problem to work, know i am trying to use it as a model for solving the poblem below. For some reason I get an error in the first line of the second function expression. I do agree something is wrong there, it's just MATLAB help is not helping, Any thoughts or suggestions?
Jack

%problem solving for
% dF/dt = -4.7316*(10^8)* (F^0.1)*(O^1.65)*exp(-15098/T)
% F(0) = .023918
%
% dO/dt = 16 * dF/dt
% O(0) = .382681
%
% dP/dt = -17 * dF/dt
% P(0) = 0
%
% dT/dt = -3715.19 * dF/dt
% T(0) = 753
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Begin myodee.m
function myodee
[t,y] = ode45(@rigid, [0 3.100], [.023918 0.38261 0 753]);
%t= the scaler time
%y= the colum vector
%ode45 is the solver
%@rigid is the function handle calling function below
%[0 3.100] is the time to be evaluated from t0 to tf
%[.023918...etc] is the initial conditions
odeprint(t,y);
%odeprint to print results on screen
%plot(t,y(:,1),'-',t,y(:,2),'-.',t,y(:,3),'.'); in case to plot


function (t,[dF,dO,dP,dT]) = rigid ([F O P T])
%[dF,dO,dP,dT] = returns outputs to y
%rigid(t,[F,O,P,T]) accepts the inputs (initial conditions)
dF = -4.7316*(10^8)*(F^0.1)*(Ox^1.65)* exp(-15098/T);
dO = 16 * dF;
dP = -17 * dF;
dT = -3715.9 * dF;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Problem Modelled after

%Begin myode.m
%function myode % Added this line
% Moved these three lines up to the top of the file
%options = odeset('RelTol', 1e-4, 'AbsTol', [1e-4 1e-4 1e-5]);
%[T,Y] = ode45(@rigid, [0 12], [0 1 1], options);
%plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.');

%function dy = rigid(t,y)
%dy = zeros(3,1);
%dy(1) = y(2) * y(3);
%dy(2) = -y(1) * y(3);
%dy(3) = -0.51 * y(1) * y(2);

 

[LeBrad]

I made some changes that should make it work for you.

function myodee
tspan = [0 3.1]; y0 = [.023918; 0.38261; 0; 753];
[t,y] = ode45(@rigid,tspan,y0);
%t= the scaler time
%y= the colum vector
%ode45 is the solver
%@rigid is the function handle calling function below
%[0 3.100] is the time to be evaluated from t0 to tf
%[.023918...etc] is the initial conditions
odeprint(t,y);
%odeprint to print results on screen
plot(t,y(:,1),'-',t,y(:,2),'-.',t,y(:,3),'.'); %in case to plot


function dy = rigid(t,y)
%[dF,dO,dP,dT] = returns outputs to y
%rigid(t,[F,O,P,T]) accepts the inputs (initial conditions)
F = y(1); 
Ox = y(2);
P = y(3);
T = y(4);
dF = -4.7316*(10^8)*(F^0.1)*(Ox^1.65)* exp(-15098/T);
dO = 16 * dF;
dP = -17 * dF;
dT = -3715.9 * dF;
dy = [dF; dO; dP; dT];

A couple things to note here.

First, the odefunc rigid which you pass to ode45 should take in two vectors (t and y) and return one vector (dy). If you try to pass the vector as [a b c d] in the function definition (first line of the function)
function (t,[dF,dO,dP,dT]) = rigid (t,[F O P T])
, it will just confuse matlab, so just pass it as y,
function dy = rigid(t,y)
and then pull it apart in the next line like
a = y(1);
b = y(2);
...
Similarly, define dy = [da db dc dd] and then return the vector dy.

It is, however, ok to write
[t,y] = ode45(@rigid, [0 3.100], [.023918 0.38261 0 753]);
when you are calling a function. But it just makes things confusing so I defined tspan and y0.

And you can't have spaces between a function name and it's argument, e.g.
rigid (a,b) should be rigid(a,b).

Let me know if this doesn't work for you or if something I changed isn't clear.

 

[J-bonz]

Thank You very much, I will give it a try.... and report back
Thanks again
Jack

 

[J-bonz]

Wow, I tried it and i worked, I just had a few more q's

1. how do increase the significant digets that are printed out on the screen by the odeprint function?

2. I changed the time scale to (0 to 3.1 ms), The question I had is between 3 and 3.1 ms there is a major /defining it as a "stiff" problem. How can I focus in on this area without initial conditions? I already added this code below to increase the precision.

options = odeset('RelTol',1e-6,'AbsTol',[1e-6 1e-6 1e-6 1e-6]);

Thanks for all your help.

Jack

 

[LeBrad]

1. how do increase the significant digets that are printed out on the screen by the odeprint function?

add 'format long' somewhere near the top of your first function.

2. I changed the time scale to (0 to 3.1 ms), The question I had is between 3 and 3.1 ms there is a major /defining it as a "stiff" problem. How can I focus in on this area without initial conditions? I already added this code below to increase the precision.

options = odeset('RelTol',1e-6,'AbsTol',[1e-6 1e-6 1e-6 1e-6]);

Are you saying you only want to look at y(t) between t=3ms and t=3.1ms? If that's what you want, you need to know the initial condition y(3ms). If you want, you could run ode45 twice, once on t=[0,3ms], and then again from t=[3ms,3.1ms] using the final value from the first run as the intial condition for the second part.

If you just want smaller time steps near 3ms, you can define the time discretization for ode45. When tspan is [0,3.1] as it is now, ode45 picks how to discretize it, but if you define tspan = [0:0.01:3.1] it will chop up the time into into 0.01 size pieces. If you want to go from 0 to 3.1 but you want to go slower from 1 to 2, you can define
tspan = [[0:0.1:1], [1:0.001:2], [2:0.1:3.1]]
Then you have timesteps of 0.1 from 1 to 2, timesteps of 0.001 from 1 to 2, and timesteps of 0.1 from 2 to 3.1. So you can do something like that to alter the timestep size in important regions.

 

[J-bonz]

LeBrad,
That sounds like it will help, I will have to give it a try tommorow
Jack

 

[yagyesh]

I have a doubt...I want to pass additional parameters in ode 45 ...how can i do it...! Pllz reply...!

 

[pramod123]

Suppose my function F itself is a function of O or P or T...Then how to code it...?
Please help me asap...

 

[pramod123]

Suppose my function F itself is a function of O or P or T...Then how to code it...?
Please help me asap...

 

[tapon.mahbub]

Hi,I am trying to solve a system of differential equations consisting 3 equations with 3 variables.one of the given eqations is of second order. I converted it to first order.everything was ok compare to the first example here given. but unfortunately I found many errors.
would you please help me to solve it...
here are the MATLAB codes...

function dx=test(r,x)
dx=zeros(4,1);
dx(1)=(1/(6*r))*(-exp(2*x(1))*(1+k*x(3)^2)+3-r^2*k*x(4)^2);
dx(2)=(1/(2*r))*(((3+k*x(3)^2)-exp(-2*x(1))*k*r^2*x(4)^2)*1/3*exp(-2*x(1))-1);
dx(4)=exp(2*x(1))*(exp(-2*x(1))*x(4)*r^2*(dx(1)-dx(2)-2/r)+2*x(3));

I saved namely test.m file
then I wrote another script file as follows...

k=1;
x0=[0 -1 0 1];
[r,x]=ode45('test',[0 10],x0);
plot(r,x(:,1))




here the reply from the matlab...

? Undefined function or variable 'k'.

Error in ==> test at 3
dx(1)=(1/(6*r))*(-exp(2*x(1))*(1+k*x(3)^2)+3-r^2*k*x(4)^2);

Error in ==> odearguments at 109
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> solve at 3
[r,x]=ode45('test',[0 10],x0);

>>

please help me to solve it!
thanks