CP502 – Advanced Fluid Mechanics Example 1
Let us consider the equation given by (2).
Matlab – Numerical Solutions of Ordinary Differential
Equations (Initial Value Problems) Step 1
Numerical Simulation using Matlab numerical routines To use the routine ode45 we need to create a Matlab file
There are number of Matlab routines to simulate a (m-file) that would compute derivative dy/dt for any given
dynamical system describe by the set of simultaneous first y and t. We want to compute the solution from t=0 to
order differential equations t=10.
f ( y, u , t ) with y(t0)=yo (1) Thus we create the following Matlab function file.
where y can be a vector of n elements if the system is of nth function dy=ex1(t,y)
order. u is the input to the system and t is the time. f is a a=-2;b=1;
vector of n functions. %computing derivatives for given y %
First Order Differential Equations dy=a*y+b;
However, we can begin the study by considering a simple and save this file in ex1.m. You may notice that we have
first order differential equation of the form selected that a=-2 and b=1.
ay b f ( y) with y(t0)=y0 (2) Step 2
dt The simplest format of ode45 is given by
where we need to find the solution from t=t0 to t=tf. [t,y] = ode45(odefun,tspan,y0)
Constants a and b will depend on the system, and the
initial value y0 we can select. where odefun is the function that provides the derivatives,
tspan a vector indicating the initial time, time step and
The analytical solution of this equation is a continuous final time, and yo is the initial value.
function y(t) defined in the domain
[t0, tf]. The solution to the differential equation can be computed
by using ode45, by typing
The numerical solution of this equation is given by a set >>[t,y]=ode45(@ex1,[0:0.1:10],5);
of numerical values
yn= y(t0 + nT) for n=0,1,2, ……..,p and T is a step size Note that the 0.1 is the step size we have selected, and the initial
decided by us. value of y is 5. the output variable of the ode45 is the vectors t
and y which contains the time instants ans the solution values
The solution can be visualized by
yn Further Investigations
- Now you can investigate the numerical method by trying
different step sizes T and comparing the solutions you get
by plotting them. As you know the small step sizes should
give more accurate results.
- You can also investigate the behaviour of the solution of
the differential equation by computing the solution for
t0 t1 t2 T tp different initial values.
- By giving different values for a and b we get different
Fig. 1 Analytical solution y(t) and numerical differential equations and hence different solutions. Try
solution yn. different values for a and b. When changing one keep the
other constant (also keep the initial value constant). Try
Matlab functions to solve the equation both positive and negative values for a and b.
To obtain the solution of the above equation we can use - Above we have specified the step size T. Matlab ode45
the Matlab functions ode23 or ode45. To get familiar with has a special feature where it can decide the step size at the
how to use them we shall do a few examples. end of each iteration. Give the command
>>[t,y]=ode45(@ex1,[0 10],5); and Matlab will use
variable step sizes. Plot the solution and compare with solution
for constant step size.
Exercise 1 system by creating a function and passing the function
The single tank system dynamics name to the ode45.
be described by qi
dh Create the function with the file name sys1.m
A h 1 qi
%to make dh a column vector
dx = zeros(2,1);
qi is the flow in, A the tank x-sectional %computing derivatives
area h is the water level. Take k=1.1 and dx(1)=x(2);
A=10 and qi= 2 and find the solution for dx(2)= -(k/m)*x(1)-(b/m)*x(2);
different initial water levels.
Investigate the behaviour of the tank for different k and A. Step 2:
You can also change the in flow qi. You may have to The simplest format of ode45 is given by
change the final time tf. [t,x] = ode45(odefun,tspan,x0)
Use Matlab help to learn more about ode45 and the where odefun is the function that provides the derivatives,
other ODE solvers available in Matlab. tspan a vector indicating the initial time and final time,
and yo is the initial value of vector x.
Second Order Differential Equations
For the above system, to use ode45, type
(a) System without an input (u=0) To begin let us >>[t,x]=ode45(@sys1,[0:0.1:10],[5 0]);
consider a system without an input u. The system is simulates from t=0 to t=10 with initial values of x1=5
described by the equation and x2=0. Step size is 0.1.
dy You can also try to plot x1 against x2.
f ( y, t ) with y(0)=yo >> plot(x(:,1),x(:,2))
dt This is called the phase plane plot. It clearly shows how
the system comes to the equilibrium from different initial
Please note that if the initial condition yo=0 then the points. Use the hold command
solution y(t) is zero for all t. In the following two examples >> hold
you can change the initial conditions and observe how the and draw the phase plane curves for different initial values.
system behaves and reach the equilibrium point.
Note the response of the system. Change the values of m,b
Example 2. and k, and try to obtain different responses of the system.
A spring, mass and damper system without an external Also observe the variation in the time step, where the step
force acting on it is described by the equation size is not constant.
mx bx kx 0 with the initial values of x’ and x
specified. This second order equation can be expressed as
a system of first order coupled differential equations as Example 3
follows (taking x1=x) The previous system is a linear system. A well known non-
linear system is the Van der Pols oscillator described by
dx1 the differential equation
x2 d 2x dx
( x 2 1) x 0
dx2 b k dt dt
x 2 x1 This can be written as two differential equations as follows
dt m m (taking x1=x)
with the initial values of x1 and x2 specified. Note that x1 x2
is the displacement and x2 is the velocity. Initially we can dt
displace the mass and let it free to move, or give an initial dx2
( x12 1) x2 x1
velocity too. dt
To simulate this system (taking m=1, b=.5 and k=3) we Write a function to pass the derivatives to the ode45 and
use the Matlab function ode45. simulate the system for different initial values of x1 and x2.
Draw a phase plane curves for this system. What is the
Step 1: difference of this oscillator compared to a simple
At each time step ode45 needs to compute the derivatives pendulum.
as needed in Runge-Kutta methods, and this is given to the
(b) System with an input ( u (t ) 0 ) Since the input generated is not available it is generated
by the following commands.
A system with a driving force or input function can be >>for i=1:length(t),u(i)=feval(‘finp’,t(i)); end
described by >> plot(t,u,t,x)
f ( y, u , t ) with y(0)=yo Investigate the behaviour of the system for different input
dt functions such as, impulse, step, random inputs etc., by
creating different input function files or having them on
The input u(t) can be a constant or a function of t. The the same file and select the appropriate one.
Matlab function ode45 can be used as before with a slight
modification. At the end of this lab you should know
1. To create a system of differential equations for a given
Example 4 higher order non-linear differential equation.
2. To write the Matlab functions to provide the derivatives
System is same as in example 1, but driven by an input. to the ode45
The respective equation is given by 3. To provide different input functions to ode45
mx bx kx f (t ) where f(t) is the input function 4. simulate a system using ode45
or the driving force. The respective system of equation is 5. plot the phase plane curves of the system and recognize
given by the behaviour of the system..
x2 Exercise 2
dx2 b k Simulate the second order interactive tank system
f (t ) x2 x1 described by
dt m m
A1 Fi a12 (h1 h2 )
Let us give a sinusoidal input to the system dt
f(t)= A sin(wt) dh
A2 2 a12 (h1 h2 ) a 2 h2
Step 1 dt
Create a Matlab function for the input function as follows Indicate clearly how you would handle the situation
during the simulation if h1<h2. Study the behaviour of
Function u=finp(t) this system for different types of inputs Fi. Can you
A=1 select A1, A2, a12 and a2 such the system would exhibit
w=2 damped oscillations for a step input?
u=A*sin(w*t) Typical set of values to begin with A1=1, A2=0.5, ,
and store in finp.m
Create a function to provide the derivatives to the ode45 in
every iteration as follows (similar to in example 1 but with
a change to accommodate the input function)
%evaluate the external function %value
for a given value of t
%to make dh a column vector
dx = zeros(2,1);
dx(2)= u -(k/m)*x(1)-(b/m)*x(2);
As before run the ode45 as follows
>>[t,x]=ode45(@sys2,[0 10],[5 0]);