Ordinary Differential Equations

1
ODE
Ordinary Differential Equations
2

• ????????????
• ??x1ltx2lt ltxn
• ?? ???

3

• ? ????(???)
• yi1yih f(xi,yi) (i 0,1, , n-1)
• ???????
• ???????
• ? ???????

4

• ???????????
• ??????f(x,y)??,?????O(h3)

5
???
function t,y Heun(ode,tspan,h,y0) t
(tspan(1)htspan(end))' n length(t) y
y0ones(n,1) for i2n k1
feval(ode,t(i-1),y(i-1)) k2
feval(ode,t(i),y(i-1)hk1) y(i)
y(i-1)h(k1k2)/2 end
6

• ? ?????(Runge-Kutta)
• ????????
• ?????? f (xi,yi) ??,?????O(h2)

7

• ??????????
• ????? f (x, y) ??,?????O(h5)

0
1/2 1/2
1/2 0 1/2
1 0 0 1
1/6 2/6 2/6 1/6
8
???????????????
???4? Runge-Kutta ???????
9
n ???Runge-Kutta ?????????
??
Adams ????(Adams-Bashforth ??)??? k1 ? k1
????? ???(k2),
???(k3),
Adams ????(Adams-Moulton ??)??? k1 ? k2
????? ???(k2),
Adams ??-????(Adams-Bashforth-Moulton
??) ?????????????????
10
include ltstdio.hgt include ltstdlib.hgt include
ltmath.hgt define TRUE 1 main() int nstep_pr,
j, k float h, hh, k1, k2, k3, k4, t_old,
t_limit, t_mid, t_new, t_pr, y, ya, yn double
fun() printf( "\n Fourth-Order Runge-Kutta
Scheme \n" ) while(TRUE) printf(
"Interval of t for printing ?\n" )
scanf( "f", t_pr ) printf( "Number of
steps in one printing interval?\n" )
scanf( "d", nstep_pr ) printf( "Maximum
t?\n" ) scanf( "f", t_limit ) y
1.0 / Setting the initial value of the
solution / h t_pr/nstep_pr
printf( "hg \n", h ) t_new 0
/ Time is initialized. / hh
h/2 printf( "------------------------------
--------\n" ) printf( " t
y\n" ) printf( "-----------------------
---------------\n" ) printf( " 12.5f
15.6e \n", t_new, y )
11
do for( j 1 j lt nstep_pr j )
t_old t_new t_new
t_new h yn y t_mid
t_old hh yn y k1
hfun( yn, t_old ) ya yn k1/2
k2 hfun( ya, t_mid ) ya yn
k2/2 k3 hfun( ya, t_mid ) ya
yn k3 k4 hfun( ya, t_new )
y yn (k1 k22 k32 k4)/6
printf( " 12.5f 15.6e \n",
t_new, y ) while( t_new lt t_limit )
printf( "-------------------------------------
-\n" ) printf( " Maximum t limit exceeded
\n" ) printf( "Type 1 to continue, or 0 to
stop.\n" ) scanf( "d", k ) if( k
! 1 ) exit(0)
double fun(y, t) float y, t float fun_v
fun_v -y return( fun_v )
12
? ????? ??????????????,??xi??,???????y(
xi1)????? ? ??xi???h??????y(xi1)? ???,?
??xi??? h/2 ???? ?????y(xi1)???????????e??????
??,?
??y(xi1)????????????
13
??? ??2????h???????h?????h??????,????h?????????
???????,?????h??????????? ?????h?,2???????????????
??????????,???????h?
14
Matlab ??
Ode23 ???, ???, ???Runge-Kutta,??? Ode45???, ???,
???Runge-Kutta,????,??? Ode113???, ???,
?????(1-13)Adams PECE ??, ?????? Ode15s ??,
???,??Gears (?BDF)??, ????. ??ode45??,
????????,???? Ode23s ??, ???, ??2?Rosenbrock?,
????, ???ode15s ?????????. Ode23t ????,
??????,?????????,?????????. Ode23tb ??, TR-BDF2,
?R-K?????????,????Gear ?. ????,
??????????????,?ode15s?.
15
Matlabs ode23 (Bogacki, Shampine)
16
Runge-Kutta-Fehlberg??(RKF45)
4?Runge-Kutta??
5?Runge-Kutta??
??????
Matlabs ode45 is a variation of RKF45
17
Van der Pol
function xdot vdpol(t,x) xdot(1)
x(2) xdot(2) -(x(1)2 -1)x(2) -x(1) xdot
xdot' VDPOL must return a column vector.
xdot x(2) -(x(1)2 -1)x(2) -x(1)
xdot 0 , 1 -1 ,-(x(1)2 -1) x t0
0 tf 20 x0 0 0.25 t,x
ode45(_at_vdpol,t0,tf,x0) plot(t,x) figure(101) p
lot(x(,1),x(,2))
18
Lorenz??? function xdot lorenz(t,x) xdot
-8/3, 0, x(2) 0, -10, 10 -x(2), 28,
-1x x0 0,0,eps' t,x
ode23('lorenz',0,100,x0) plot3(x(,1),x(,2),x(
,3)) plot(x(,1),x(,2))
19
(No Transcript)
20
????
??????(Gears method) ??Runge-Kutta?
function yp examstiff(t,y) yp -2, 1 998,
-998y 2sin(t)999(cos(t)-sin(t)) y0
23 tic,t,y ode113('examstiff',0,10,y0)t
oc tic,t,y ode45('examstiff',0,10,y0)toc ti
c,t,y ode23('examstiff',0,10,y0)toc tic,t,
y ode23s('examstiff',0,10,y0)toc tic,t,y
ode15s('examstiff',0,10,y0)toc tic,t,y
ode23t('examstiff',0,10,y0)toc tic,t,y
ode23tb('examstiff',0,10,y0)toc
21
????ODE
????ODE(Matlab7)
Weissinger??
function f weissinger(t,y,yp) f ty2yp3 -
y3yp2 t(t21)yp - t2y t0 1 y0
sqrt(3/2) yp0 0 guess y0,yp0
decic(_at_weissinger,t0,y0,1,yp0,0)
?????????y0?? t,y ode15i(_at_weissinger,1,10,y0
,yp0) ytrue sqrt(t.20.5) plot(t,ytrue,t,y,'o
')
22
function yp ddefun(t,y,Z) yp zeros(2,1)
define lags1,3 yp(1) Z(1,2)2
Z(2,1)2 yp(2) y(1) Z(2,1) function y
ddehist(t) y zeros(2,1) y(1) 1 y(2) t-2
??????
Sol dde23(ddefun,lags,ddehist,tspan)
lags 1,3 sol dde23(_at_ddefun,lags,_at_ddehist,0
,1) hold on plot(sol.x,sol.y(1,),'b-') plot(s
ol.x,sol.y(2,),'r-')
?? 
23
????? ????????
????
bvp4c
24
?????????? ????????(11)??????????????????
(IVP1)
(IVP2)
?????????????
???????????
25
???? y dsolve('D2y -a2y','x') ??? y
dsolve('D2y -a2y','y(0)1','Dy(pi/a)0','x')
x,y dsolve('Dx4x-2y','Dy2x-y','t')
Pdetool ????,????,????,??,??,??????
??
???
??????
26
(No Transcript)
27
Stokes ??
Q1-P0?????
28
(No Transcript)
29
(No Transcript)
30
Navier-Stokes ??
31
MAC????
32
(No Transcript)
33
(No Transcript)
34
(No Transcript)
35
????(A Mach 3 Wind Tunnel with a
Step) ????????????????????1.0,?3.0,???0.2,???????
??0.6???????
?????????
36
Sod?? Sod??????????????,????????,?????????????,??
???,?????????? Euler???
37
A picture is worth a thousand words.

-
Anonymous
Make it right before you make it faster. - Brain
W. Kernighan, P. J. Plauger, The Elements of
Programming Style(1978)