Numerical Methods For Differential Equations: Euler Method
Numerical Methods For Differential Equations: Euler Method
Numerical Methods For Differential Equations: Euler Method
1
y ( 0) 2 C 2
true solution y (t ) 2e 10t , time constant 1 10 0.1
Matlab:
clear;
r = -10; delta = 0.02; y(1) = 2;
k = 0;
for time = [delta: delta:0.5]
k = k + 1;
y(k+1) = y(k)+r*y(k)*delta;
end
t = [0:delta:0.5];
y_true = 2*exp(-10*t);
plot(t,y,’o’,t,y_true);
Simulation results:
1.8
1.6
1.4
1.2
y trajectory
0.8
0.6
0.4
0.2
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Time(sec.)
There is some noticeable error from simulation results. If we use a step size equal to
5% of the time constant ( t 0.005 ), the error would not be noticeable on the plot.
Example 2: To illustrate the difficulties caused by an oscillating solution, consider the
following equation
y sin t
2
for y(0) = 0 and 0 t 4 .
True solution:
dy
sin t
dt
dy sin dt
y cos t c
y (0) 0 c 1
y (t ) 1 cos t
Matlab:
clear;
delta=2*pi/13; y(1)=0;
k=0;
for time=[delta:delta:4*pi]
k=k+1;
y(k+1)=y(k)+sin(time-delta)*delta;
end
t_true=[0:delta/10:4*pi];
y_true=1-cos(t_true);
t=[0:delta:4*pi];
plot(t,y,’o’,t_true,y_true);
end
3
2
1.8
1.6
1.4
1.2
y trajectory
0.8
0.6
0.4
0.2
0
0 2 4 6 8 10 12 14
Time(sec.)
“x”: t 2 26
“o”: t 2 13
However, very small step sizes require longer run times and can result in a large
accumulated error because of round-off effects. So we seek better algorithms to use
for more challenging applications.
We now consider two methods that are more powerful than the Euler method. The
first is a predictor-corrector method. The second method is the Runge-Kutta family of
algorithms.
Predictor-error method (Modified Euler method)
The Euler method can have a serious deficiency in problems where the variables are
rapidly changing because the method assumes the variables are constant over the time
interval t . One way of improving the method is to use a better approximation to the
right side of the equation
dy
f (t , y )
dt
The Euler approximation is
y (t k 1 ) y (t k ) t f [t k , y (t k )]
Suppose instead we use the average of the right side of equation on the interval
(t k , t k 1 ) . This gives
4
t
y (t k 1 ) y (t k ) f k f k 1
2
where
f k f (t k , y (t k ))
f k 1 f (t k 1 , y (t k 1 ))
So, use Euler method to obtain y (t k 1 ) such that predictor-error method can be
implemented.
In summary,
Euler predictor: x k 1 y k h f (t k , y k )
h
Predictor-error: y k 1 y k f (t k , y k ) f (t k 1 , x k 1 )
2
Example 3: y ry , r 10, y (0) 2 , plots the y trajectory over the range
0 t 0 .5
Step size t 0.02 , which is 20%of the time constant. (the same as Example 1)
Matlab:
clear;
r = -10; delta = 0.02; y(1) = 2;
k = 0;
for time = [delta: delta:0.5]
k = k + 1;
x(k+1)=y(k)+r*y(k)*delta;
y(k+1)= y(k)+(delta/2)*(r*y(k)+r*x(k+1));
end
t = [0:delta:0.5];
y_true = 2*exp(-10*t);
plot(t,y,’o’,t,y_true);
5
2
1.8
1.6
1.4
1.2
y trajectory
0.8
0.6
0.4
0.2
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Time(sec.)
We can find that there is less error than with the Euler method using the same step
size.
Example 4: the same as Example 2
clear;
delta=2*pi/13; y(1)=0;
k=0;
for time=[delta:delta:4*pi]
k=k+1;
y(k+1)=y(k)+(delta/2)*(sin(time-delta)+sin(time));
end
t_true=[0:delta/10:4*pi];
y_true=1-cos(t_true);
t=[0:delta:4*pi];
plot(t,y,’o’,t_true,y_true);
end
6
2
1.5
1
y trajectory
0.5
-0.5
0 2 4 6 8 10 12 14
Time(sec.)
From above simulation results, we can find that the modified Euler method has less
error than with the Euler method, but there is still some noticeable error near the
peaks and valleys where the solution is rapidly changing. This error can be decreased
by decreasing the step size.
Runge-Kutta methods
The Runge-Kutta algorithm lets us solve a differential equation numerically (ie.
approximately) without require the calculation of partial derivatives. Consider the
single variable problem x' = f (t, x) with initial condition x(0) = x0. Suppose that xn is
the value of the variable at time tn.
The Runge-Kutta formula takes xn and tn and calculates an approximation for xn+1 at a
brief time later, tn+h. It uses a weighted average of approximated values of f (t, x) at
several times within the interval (tn, tn+h). The formula is given by
xn+1 = xn + h⁄6 (k1 + 2k2 + 2k3 + k4) where k1 = f (tn, xn)
k2 = f (tn + h⁄2, xn + h⁄2 k1)
k3 = f (tn + h⁄2, xn + h⁄2 k2)
k4 = f (tn + h, xn + h k3)
To run the simulation, we simply start with x0 and find x1 using the formula above.
Then we plug in x1 to find x2 and so on.
The Runge-Kutta algorithm is known to be very accurate and well-behaved for a wide
7
range of problems.
The Runge-Kutta algorithm can be easily extended to a set of first order differential
equations. You wind up with essentially the same formulas shown above, but all the
variables (except for time) are vectors.
To give an idea of how it works, here is an example where we expand the vector
notation. That is, instead of using the highly compact vector notation, we write out all
the components of the vector.
Suppose there are only two variables, x, y and two differential equations
x' = f (t, x, y)
y' = g (t, x, y)
Again we let xn be the value of x at time tn, and similarly for yn. Then the formulas for
the Runge-Kutta algorithm are
xn+1 = xn + h⁄6 (k1 + 2k2 + 2k3 + k4)
yn+1 = yn + h⁄6 (j1 + 2j2 + 2j3 + j4) where k1 = f (tn, xn, yn)
j1 = g (tn, xn, yn)
k2 = f (tn + h⁄2, xn + h⁄2 k1, yn + h⁄2 j1)
j2 = g (tn + h⁄2, xn + h⁄2 k1, yn + h⁄2 j1)
k3 = f (tn + h⁄2, xn + h⁄2 k2, yn + h⁄2 j2)
j3 = g (tn + h⁄2, xn + h⁄2 k2, yn + h⁄2 j2)
k4 = f (tn + h, xn + h k3, yn + h j3)
j4 = g (tn + h, xn + h k3, yn + h j3)
Given starting values x0, y0 we can plug them into the formula to find x1, y1. Then we
can plug in x1, y1 to find x2, y2 and so on.
Derivation of the second-order Runge-Kutta:
Consider
x f (t , x (t )) , t 0, x(0) x 0*
x n 1 x n h a1 f (t n , x n ) a 2 f (t n h, x n f (t n , x n ) (a)
Choose the a1 , a 2 , , such that the right-hand side of above equation approximates
the right-hand side of the following second-order Taylor series with an error o( h 2 ) .
h2
x n 1 x n h f (t n , x n ) f x (t n , x n ) f (t n , x n ) f t (t n , x n ) (b)
2
Review: linearization for a given equilibrium point
x 1 f 1 ( x1 , x 2 )
x 2 f 2 ( x1 , x 2 )
Equilibrium Point ( x1e , x 2 e )
Taylor series:
8
f 1 f 1
f 1 ( x1 , x 2 ) f 1 ( x1e , x 2 e ) ( x1 x1e ) x1 x1e ( x 2 x 2e ) x1 x1e H .O.T
x1 x2 x 2 e x 2 x 2 x2 e
f 2 f 2
f 2 ( x1 , x 2 ) f 2 ( x1e , x 2 e ) ( x11 x1e ) x1 x1 e ( x 2 x 2e ) x1 x1e H .O.T
x1 x2 x2 e x 2 x2 x2 e
x n 1 x n h a1 f (t n , x n ) a 2 f (t n h, x n f (t n , x n )
h a1 f (t n , x n ) a 2 f (t n , x n ) hf t (t n , x n ) f x (t n , x n ) f (t n , x n ) (c)
a1 a 2 hf (t n , x n ) a 2h f t (t n , x n ) a 2 hf x (t n , x n ) f (t n , x n )
2
If set a 2 1 / 2 , then a1 1 / 2, 1, h .
h h
x n 1 x n f (t n , x n ) f (t n h, x n h f (t n , x n )
2 2
Therefore, we have the following formulas for second-order Runge-Kutta method:
h
x n 1 x n k1 k 2
2
k1 f (t n , x n )
k 2 f (t n h, x n h k1 )
Example 5: same as Example 1
Matlab:
clear;
r = -10; delta = 0.02; y(1) = 2;
k = 0;yy(1)=2;
for time = [delta: delta:0.5]
k = k + 1;
deltah=delta/2;
k1=-10*y(k);
9
k2=-10*(y(k)+(deltah*k1));
k3=-10*(y(k)+(deltah*k2));
k4=-10*(y(k)+(delta*k3));
y (k+1)= y(k)+(delta/6)*(k1+2*k2+2*k3+k4);
end
t = [0:delta:0.5];
y_true = 2*exp(-10*t);
plot(t,y,'o',t,y_true);
1.8
1.6
1.4
1.2
y trajectory
0.8
0.6
0.4
0.2
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Time(sec.)
10
k3=sin(time_t+deltah);
k4=sin(time_t+delta);
y (k+1)= y(k)+(delta/6)*(k1+2*k2+2*k3+k4);
end
t_true=[0:delta/10:4*pi];
y_true=1-cos(t_true);
t=[0:delta:4*pi];
plot(t,y,'o',t_true,y_true);
end
Simulation results
1.5
1
y trajectory
0.5
-0.5
0 2 4 6 8 10 12 14
Time(sec.)
11
y f (t , y )
Solver Syntax:
[t,y]=ode23(‘ydot’, tspan, y0)
ydot: name of the function file whose input must be t and y and whose output must be
a column vector representing dy/dt; that is, f(t,y).
tspan: is [t0,tf], where t0 and tf are the desired starting and ending values of the
independent parameter t.
y0: is the initial value of y(t).
Example 7: same as Example 1.
Matlab file: example1.m
function ydot=example1(t,y)
ydot=-10y;
Matlab file: test.m
clear;
[t,y]=ode45(‘example1’,[0,0.4],2);
y_true=2*exp(-10*t);
plot(t,y,’o’,t,y_true)
% note that we need not generate the array t to evaluate y_true, because t is generated
%by the ode45 function.
1.8
1.6
1.4
1.2
y trajectory
0.8
0.6
0.4
0.2
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
Time(sec.)
12
2
1.8
1.6
1.4
1.2
y trajectory
0.8
0.6
0.4
0.2
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
Time(sec.)
x 1 0 1 x 0
x 4 7 1 1 f (t ) Ax bf
x 2 5 x
5 2 5
Suppose that f (t ) sin t .
Matlab: example8.m
function xdot=example8(t,x)
13
xdot(1)=x(2);
xdot(2)=(1/5)*(sin(t)-4*x(1)-7*x(2));
xdot=[xdot(1);xdot(2)];
or
function xdot=example8(t,x)
xdot=[x(2); (1/5)*(sin(t)-4*x(1)-7*x(2))];
or
function xdot=example8(t,x)
A=[0 1;-(4/5) –(7/5)];
b=[0;(1/5)];
xdot=A*x+b*sin(t);
Matlab: test.m
[t,x]=ode45(‘example8’, [0,6],[3,9]);
plot(t,x(:,1);
4
y trajectory
-1
0 1 2 3 4 5 6
Time(sec)
14
x 1 0 1 x1 0
x 2 3 x 1 u
2 2
x
y 1 0 1
x2
x Ax bu
y cx
fuzzy controller:
Control rules:
If e P & ce P, Thenc u P
If e P & ce N , Thenc u Z
If e N & ce P, Thenc u Z
If e N & ce N , Thenc u N
Membership functions:
N P
-4 -2 2 4 e,ce
N Z P
-2 0 2 cu
Objective: check its step response (1) using ode45, (2) using Runge-Kutta, (3) use
Euler method
15
global r u
u=0;
r=1;
[t,x]=ode45('example8', [0,8],[0,0]);
plot(t,x(:,1));
Matlab: fuzz.m
function y=fuzz(e,ce)
% if e=P &ce=P Then cu=P
r1=P(e)*P(ce);
% if e=P &ce=N Then cu=Z
r2=P(e)*N(ce);
% if e=N &ce=P Then cu=Z
r3=N(e)*P(ce);
% if e=N &ce=N Then cu=N
r4=N(e)*N(ce);
y=(r1*2+r2*0+r3*0+r4*(-2))/(r1+r2+r3+r4);
Matlab: N.m
function y=N(x)
if x<=-2, y=1; end;
if x>-2 & x<=2, y=(-1/4)*(x-2);end;
if x>2, y=0; end;
Matlab: P.m
function y=P(x)
if x<=-2, y=0; end;
if x>-2 & x<=2, y=(1/4)*(x+2);end;
if x>2, y=1; end;
16
1.4
1.2
0.8
y trajectory
0.6
0.4
0.2
0
0 1 2 3 4 5 6 7 8
Time(sec)
17
j3=(-2*x11-3*x2+u)+(h/2)*j2;
j4=(-2*x11-3*x2+u)+h*j3;
x2=x2+((h/6)*(j1+2*j2+2*j3+j4));
%================================
e=r-x1;
ce=0-x2;
xx1(i)=x1;
end
tt=[1:mm]*h;
plot(tt,xx1);
0.9
0.8
0.7
0.6
y trajectory
0.5
0.4
0.3
0.2
0.1
0
0 1 2 3 4 5 6 7 8
Time(sec)
18
ce=0-x2;
xx1(i)=x1;
end
tt=[1:mm]*dt;
plot(tt,xx1);
0.9
0.8
0.7
0.6
y trajectory
0.5
0.4
0.3
0.2
0.1
0
0 1 2 3 4 5 6 7 8
Time(sec)
They have the same results in this example with the Euler method and Runge-Kutta
method. However, there is slightly different in using ode45 solver.
Fuzzy rules and the plant model are given the same as in Example 9.
15 5
P
5 5
We simulate it using (1) Euler method, (2) Runge-Kutta method, (3) ode solver.
19
r1=Pa(e)*Pa(ce);
% if e=P &ce=N Then cu=Z
r2=Pa(e)*Na(ce);
% if e=N &ce=P Then cu=Z
r3=Na(e)*Pa(ce);
% if e=N &ce=N Then cu=N
r4=Na(e)*Na(ce);
dd=r1+r2+r3+r4;
eta=[r1/dd r2/dd r3/dd r4/dd];
%========
% pp(4*1)
y=eta*pp;
Matlab: example 10
clear;
mm=12000;
dt=0.001;
x1=0;x2=0;u=0;
r=1;e=r-x1;ce=0-x2;
P=[15 5;5 5];
pp=[0 0 0 0];
b=[0 1]';
eta=[0 0 0 0];
gamad=1;
for i=1:mm
[u,eta]=fuzz(x1,x2,pp');
x1=x1+x2*dt;
x2=x2+(-2*x1-3*x2+u)*dt;
e=r-x1;
ce=0-x2;
ee=[e ce];
%ee(1*2);P(2*2);b(2*1);eta(1*4)
pp=pp+(gamad*ee*P*b*eta)*dt;
xx1(i)=x1;
end
tt=[1:mm]*dt;
plot(tt,xx1);
Matlab:
function y=Na(x)
20
if x<=-2, y=1; end;
if x>-2 & x<=2, y=(-1/4)*(x-2);end;
if x>2, y=0; end;
Matlab:
function y=Pa(x)
if x<=-2, y=0; end;
if x>-2 & x<=2, y=(1/4)*(x+2);end;
if x>2, y=1; end;
0.8
y trajectory
0.6
0.4
0.2
0
0 2 4 6 8 10 12
Time(sec)
21
[u,eta]=fuzz(x1,x2,pp');
k1=x2;
k2=x2+((h/2)*k1);
k3=x2+((h/2)*k2);
k4=x2+(h*k3);
x11=x1;
x1=x1+((h/6)*(k1+2*k2+2*k3+k4));
%================================
j1=(-2*x11-3*x2+u);
j2=(-2*x11-3*x2+u)+(h/2)*j1;
j3=(-2*x11-3*x2+u)+(h/2)*j2;
j4=(-2*x11-3*x2+u)+h*j3;
x2=x2+((h/6)*(j1+2*j2+2*j3+j4));
%===================================
e=r-x1;
ce=0-x2;
ee=[e ce];
%ee(1*2);P(2*2);b(2*1);eta(1*4)
n1=(gamad*ee*P*b*eta);
n2=(gamad*ee*P*b*eta)+(h/2)*n1;
n3=(gamad*ee*P*b*eta)+(h/2)*n2;
n4=(gamad*ee*P*b*eta)+h*n3;
pp=pp+((h/6)*(n1+2*n2+2*n3+n4));
xxx1(i)=x1;
end
tt=[1:mm]*h;
plot(tt,xxx1);
22
1
0.8
y trajectory
0.6
0.4
0.2
0
0 2 4 6 8 10 12
Time(sec)
23
P=[15 5;5 5];
eta=[0 0 0 0];
wdot=[0 0 0 0];
gamad=1;
u=1;
to=0;
r=1;
[t,x]=ode23('example8', [0,12],[0,0]);
plot(t,x(:,1));
Matlab:
function [y,eta]=fuzz(e,ce,pp)
global eta
% if e=P &ce=P Then cu=P
r1=Pa(e)*Pa(ce);
% if e=P &ce=N Then cu=Z
r2=Pa(e)*Na(ce);
% if e=N &ce=P Then cu=Z
r3=Na(e)*Pa(ce);
% if e=N &ce=N Then cu=N
r4=Na(e)*Na(ce);
dd=r1+r2+r3+r4;
eta=[r1/dd r2/dd r3/dd r4/dd];
%========
% pp(4*1)
y=eta*pp;
24
1
0.8
y trajectory
0.6
0.4
0.2
0
0 2 4 6 8 10 12
Time(sec)
1.2
0.8
0.6
y trajectory
0.4
0.2
-0.2
-0.4
-0.6
0 2 4 6 8 10 12
Time(sec)
25
Matlab:
function y=fuzzs(e)
% if e=P &ce=P Then cu=P
r1=Pa(e(1))*Pa(e(2));
% if e=P &ce=N Then cu=Z
r2=Pa(e(1))*Na(e(2));
% if e=N &ce=P Then cu=Z
r3=Na(e(1))*Pa(e(2));
26
% if e=N &ce=N Then cu=N
r4=Na(e(1))*Na(e(2));
%pp=[ .1 .1 .1 .1]';
dd=r1+r2+r3+r4;
eta=[r1/dd r2/dd r3/dd r4/dd];
%========
y=eta;
27