Introduction to scientific Computing A Matrix Vector Approach Using Matlab Written by Charles FVan Loan 陈文斌 Wbchen(fudan. edu. cn 复日大学
Introduction to Scientific Computing -- A Matrix Vector Approach Using Matlab Written by Charles F.Van Loan 陈 文 斌 Wbchen@fudan.edu.cn 复旦大学
Chapter 9 The Initial value Problem Ba asic concepts The Runge-Kutta Methods The adams methods
Chapter 9 The Initial Value Problem • Basic concepts • The Runge-Kutta Methods • The Adams Methods
()=-5 Solutions to y()=-5 y(t) 0.6 0.2 0.150.2 0.35
-0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Solutions to y'(t) = -5 y(t) y'(t) = −5t
Euler method Five Steps of Euler Method (y'=-5y, y(o)=1) n+1 n+h,f(tn,yu) 16 0.8 9105 0.0 0.1 0.15 0.2 0.25 0.35
-0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Five Steps of Euler Method (y'=-5y, y(0)=1) ( , ) n 1 n n n n y = y + h f t y + Euler Method
Theorem 9 Assume that, for n=0: N, a function y, (t) exists that solves the ivp y(t=f(t, y(t), y(t,=y where(to,yo),…,(tN,y) are given and t<t1<…< Define the global error by gn=yo(tn)-y and the local truncatio n error by LTEM=Ym(tn)-y If for all tElto, tN]and none of the trajector ies {(t,yn():≤t≤N} =0 intersect. then for n=1: N gn≤∑|LTEk k=1
= − = = = = − = − = = = 1 0 0 N 1 0 0 0 N 0 1 | | | | intersect, then for 1: {( , ( )) : } 0 : for all [ , t ] and none of the trajectories 0 ( , ) If ( ) . and the local truncatio n error by Define the global error by ( ) where ( , y ),...,( , y ) are given and . ( )), ( ) solves the IVP Theorem 9 Assume that, for 0 : , a function (t) exists that k n k n y n n n n n n n N N n n n g LTE n N t y t t t N n n t t y f t y f LTE y t y g y t y t t t t t y'(t) f(t,y t y t y n N y
Euler Solution(h=0. 1)of y'=-5y, y(o)=1 y"=-5.0y,y(0)=1 0= computed solution 04 -0.05 0.0 0.15 0.3 0.3 Showtrunc m
-0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 y' = -5.0y, y(0) = 1 * = exact solution o = computed solution Euler Solution (h=0.1) of y'=-5y, y(0) = 1 Showtrunc.m
h=(t max t0)/ h LTEM cm)-y∑E=M2N=m oM.h max M,≤tol 2N n=ceil(((tmax-tO)2*M2)/(2*tol)+1
h = (tmax −t 0 )/ N 2 | | 2 2 h LTEn M M h h t t y t y LTE M N N n N n 2 1 max 0 2 max 2 2 2 | ( ) | | | = − − = = M tol N t t − 2 2 max 0 2 ( ) n = ceil(((tmax-t0)^2*M2)/(2*tol))+1;
Fixed h Euler frror for v'=-V. v(0)=1 tol=01 vals,ovals FixedEuler(fl, 1, 0, 5, 1, tol) Showfixedeuler m 15 ol=0010.n=1251 be taken in regions where the solution is smoo er step Sizes can It is better to determine h adaptively so that long
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 2 4 6 8 x 10-4 Fixed h Euler Error for y'=-y, y(0) = 1 tol = 0.010, n = 1251 Showfixedeuler.m tol=.01; [tvals,yvals] = FixedEuler('f1',1,0,5,1,tol); It is better to determine h adaptively so that longer step sizes can be taken in regions where the solution is smooth
Euler in 3-Digit Arithmetic 0.045 . o4 Note that the error gets worse as h get smaller because the step size are in the neighborhood of 0.035 unit roundoff. 0.03 0.025 0.02 h=1/180 0.015 h=1/160 0.01 0.005 h=1/140 0.2 0.3 04 0.6 0.7 0.8 0.9
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 Euler in 3-Digit Arithmetic h = 1/140 h = 1/160 h = 1/180 Note that the error gets worse as h get smaller because the step size are in the neighborhood of unit roundoff
Stability y(t)=-10y(t) yn+1=(1-10h)yn 1-10h1→h<1/5
Stability y'(t) = −10y(t) n n y (1 10h)y +1 = − |1−10h |1 h 1/5