pk = xk + hf (tk, xk), tk+1 = tk + h, h
xk+1 = xk + 2 [f (tk, xk) + f (tk+1, pk)] .
4
1.3 Matlab Builtin ODE Solvers
In addition there are many other methods for approximating solutions to ordinary differential equations, but due to a lack of time left in the semester I will just introduce you to Matlabs built-in Runge-Kutta solver ode45 and show you how it works.
s1 = f (tk, xk)
h
h
s2 = f tk + 2 , xk + 2 s1
xk+1 = xk + hs2
tk&s obtained as follows: First we replace the approximation x(t + h) − x(t)
next section. They are all of higher order than Euler’s method.
1.2 Modified Euler’s Methods
Once again we consider (1) and apply the fundamental theorem of calculus to get
t2
t2
f (t, x(t)) dt = x (t) dt = x(t2) − x(t1).
t1
t1
This implies that
t2
x(t2) = x(t1) + f (t, x(t)) dt.
t1
Now we consider two possible ways to approximate the integral on the right above.
x (t) = e−1/2 t2 .
We build the m-file eul1d.m and save it to disk (use cd T:\ to change to the T drive)
% Euler’s method 1 clear all defn=inline(’-t*x’,’t’,’x’); a=0; b=5; N=100; xa=1; h=(b-a)/N; t=a+h*(0:N); LT=length(t); x(1)=xa;
−tx(t)
x (t) =
, x(0) = 1, x (0) = 1, on 0 < t < 5.
(x (t)2 + 1)
Now we build an m-file eul2d.m to solve the problem.
clear % Euler’s method 2d clear defn=inline(’[y;-t*x/(y^2+1)]’,’t’,’x’,’y’); a=0; b=5; N=200; h=(b-a)/N; t=a+h*(0:N); LT=length(t); xa=1; xpa=1; w(:,1)=[xa;xpa];
x = f (t, x), t ∈ (a, b)
(1)
x(a) = xa
(2)
we observe that defining tj = a + hj, for j = 0, 1, 2, · · · , N where h = (b − a)/N we have
x(tj+1) − h
x(tj )
≈
f (tj, x(tj)).
for j=2:LT w(:,j)=w(:,j-1)+h*defn(t(j-1),w(1,j-1),w(2,j-1));
end x=w(1,:); plot(t,x,’LineWidth’,2)
grid on 2
We plot the result of executing this Matlab code in the following figure.
x = f (t, x, x ), t ∈ [a, b]
x(a) = xa x (a) = xpa
If we set y = x , then the sxstem can be written as
d dt
x y
=
y f (t, x, y)
x(a) = xa
y(a) = xpa
The numerical approximation can now be carrried out just as above. As we have already mentioned it is very difficult to find exact solutions to most interesting problems. The next example is a second order initial value problem. Consider
Math 4330 Sec. 1, Matlab Assignment # 4 , April 26, 2006 Name
1 Numerical Solution of ODEs Using Matlab
1.1 Euler’s Method
Euler’s one step method is undoubtedly the simplest method for approximating the solution to an ordinary differential equation. It goes something like this: Given a first order initial value problem
The basic syntax for using ode45 is the following
[t,x] = ode45(’xdot’,tspan,x0)
where xdot.m is a function file containing the the right hand side of the differential equation, tspan is a vector of time values at which you want to obtain the solution and x0 is the initial condition.
x (t) = h
3
by the more accurate approximation (see the figure)
h x(t + h) − x(t)
x t+ =
2
h
Then we apply Taylor’s formula of order one
g(ξ) = g(a) + g (a)(ξ − a) + O(((ξ − a)2),
xj+1 = xj + f (tj, xj), j = 0, 2, · · · , (N − 1).
As an example let us consider the IVP
x (t) = −tx(t), x(0) = 1 on 0 < t < 5.
The exact solution is given by
The Midpoint Rule:
The midpoint method uses Euler to step halfway across the interval, evaluates the function at this
intermediate point, then uses that slope to take the actual step.
2
2
2
h
h
h
x(t + h) ≈ x(t) + f t + , x(t) + f (t, x(t)) .
2
2
2
The Trapezoid Rule: For the trapezoid rule we approximate the integral with step size h = t2 − t1 using
h x(t2) ≈ x(t1) + 2 [f (t1, x(t1)) + f (t2, x(t2))].
Unfortunately, we do not know x(t2) on the right side, so we use, for example, Euler’s method to approximate it: x(t2) ≈ x(t1) + hf (t1, x(t1)). We thus obtain a numerical procedure by denoting x(tk) by yk for each k and setting
for j=2:LT x(j)=x(j-1)+h*defn(t(j-1),x(j-1));
end xact_sol=exp(-t.^2/2); plot(t,x,’b’,t,xact_sol,’r--’,’LineWidth’,2)
Execute the file in the Workspace by typing the first name of the file without the .m extension. We plot the result of executing this Matlab code, on 0 ≤ t ≤ 5 with initial condition x(0) = 1, in the following figure.