MATlAB 数值求解ODE方程

−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];
t2
t2
f (t, x(t)) dt = x (t) dt = x(t2) − x(t1).
t1
Biblioteka Baidu
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.
s1 = f (tk, xk)
h
h
s2 = f tk + 2 , xk + 2 s1
xk+1 = xk + hs2
tk+1 = tk + h
This formula is obtained as follows: First we replace the approximation x(t + h) − x(t)
Lets consider the first example considered above
x (t) = −tx(t), x(0) = 1 on 0 < t < 2.
We build the function file xdot.m given by function dx=xdot(t,x) dx = -t*x;
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
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
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
E(x(b), h) = Ch
where C is a constant that depends on f (t, x). There are many variations on the Euler method. Some of them are given in the first exercise in the
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
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.
to the left hand side with ξ = t + h/2, a = h/2 and g = x, to obtain
h
h
h
h
h
x t + ≈ x(t) + x (t) = x(t) + f t + , y t +
.
2
2
2
2
2
So we finally arrive at
h
h
h
x t + ≈ x(t) + x (t) = x(t) + f (t, x(t)).
2
1.5
1
0.5
0
−0.5
−1
−1.5
−2
0
1
2
3
4
5
For any numerical method the Final Global Error (FGE) is defined as
E(x(b), h) = |x(b) − xN |.
It can be shown that Euler’s method is order one in h, i.e.,
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)).
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 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.
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.
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.
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;
Therefore we can write
x(tj+1) ≈ x(tj) + h f (tj, x(tj)).
With this we can define an iterative scheme for computing approximations xj for x(tj) using the iterative scheme
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
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),
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.
1
1
0.8
0.6
0.4
0.2
0
0
1
2
3
4
5
Exact Solution dotted Line
This method can also be applied to higher dimensional problems as demonstrated in the following
example. Consider the second order initial value problem