# Solving four ordinary dif equations simultaneously using ode45 matlab

• MATLAB
• J-bonz
In summary, Jack is trying to replicate an example problem from an "intro to combustion book" for school, but there is no MATLAB class offered and he is struggling to get the example code to run. He is also attempting to use it as a model for a problem he needs to solve involving four equations simultaneously. LeBrad suggests breaking the code into two separate pieces, one for the function and one for the solver, and makes some changes to the code to help it run correctly.
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

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
J-bonz said:
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
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);
saved as rigid.m, and then you need another file which contains
Code:
 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.

Thanks, I will give that a try LeBrad

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);

I made some changes that should make it work for you.
Code:
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.

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

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]);

Jack

J-bonz said:
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]);
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.

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

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

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

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

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.
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);

>>

thanks

## 1. What is ode45 in Matlab and how does it solve differential equations?

Ode45 is a built-in function in Matlab that is used to numerically solve ordinary differential equations (ODEs) of the form dy/dt = f(t,y). It uses a combination of fourth and fifth order Runge-Kutta methods to approximate the solution of the ODE at a specified time interval. This allows for the solution of a wide range of ODEs, including systems of ODEs.

## 2. How do I input multiple equations into ode45?

In order to solve multiple ODEs simultaneously, you will need to create a function file that contains all of the equations in vector form. This function file can then be called in the ode45 function, with the initial conditions and time interval specified.

## 3. Can ode45 handle stiff systems of ODEs?

Yes, ode45 is designed to handle both non-stiff and stiff systems of ODEs. When solving stiff systems, it automatically adjusts the time step to maintain accuracy and stability.

## 4. What are the limitations of ode45 for solving ODEs?

Ode45 is limited to solving systems of ODEs that can be expressed in the form dy/dt = f(t,y). It may also have difficulty handling higher order ODEs or systems with discontinuities or singularities.

## 5. How do I interpret the output of ode45?

Ode45 returns a time vector and a solution vector for each of the equations in the system. These vectors can be used to plot the solution of the ODEs over the specified time interval. Additionally, ode45 also provides the option to output the derivative at each time step, which can be useful for further analysis of the solution.

• MATLAB, Maple, Mathematica, LaTeX
Replies
5
Views
1K
• MATLAB, Maple, Mathematica, LaTeX
Replies
6
Views
2K
• MATLAB, Maple, Mathematica, LaTeX
Replies
5
Views
2K
• MATLAB, Maple, Mathematica, LaTeX
Replies
8
Views
2K
• MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
1K
• MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
2K
• Engineering and Comp Sci Homework Help
Replies
1
Views
1K
• Introductory Physics Homework Help
Replies
5
Views
3K
• Engineering and Comp Sci Homework Help
Replies
2
Views
944
• MATLAB, Maple, Mathematica, LaTeX
Replies
2
Views
1K