MATLAB 數值微積分與微分方程式求解

1 MATLAB x = a x = b f(x) ( )baf x dx (a) (b) f(x) trapzEx:>>x=[0 10 20 30 40];>>y=[ ]; >>area=trapz(x,y) % area = = trapz(x, y)x : x y : f(x) ( ) f(x) quad ,quadlI = quad(@fun, a, b)I = quadl(@fun, a, b)( )( )fun function m-file a b Ex:1. edit y=fun(x)y=exp(-x).*cos(x);2. ( MATLAB Command Window)area=quadl(@fun,0,1) area=quadl( exp(-x).)

2 *cos(x) ,0,1)NOTE: (.* ./ .^)10cos( )xex dx ( trapz ) x2 diff Ex:>>x=0 :1;>>y=[ ];>>dx=diff(x);>>dy=diff(y);>>dydx=diff( y)./diff(x) Ex:[]( ) sin( ),( ) ?0,f xxf xx == >> x = linspace(0,pi,20);>> y = sin(x);>> d = diff(y)./diff(x);% backward or forward difference>> dc = (y(3:end)-y(1:end-2))./(x(3:end)-x(1:end -2));% central difference>> dy = cos(x);% >> plot(x, dy, x(2:end), d,'o', x(1:end-1), d,'x', x(2:end-1), dc,'^')>> xlabel('x').

3 Ylabel('Derivative') (increment function) yi yi+1 h (step size) (one-step method) i dydt=f t,y()hyyii +=+1 ti dydtti=f ti,yi()(),iif t y = (Euler's method) y t h ()1,iiiiyyf t y h+=+ (RK) (Runge-Kutta methods) n ( ) a p q k =a1k1+a2k2+ +ankn k1=f ti,yi()k2=f ti+p1h,yi+q11k1h()k3=f ti+p2h,yi+q21k1h+q22k2h() kn=f ti+pn 1h,yi+qn 1,1k1h+qn 1,2k2h+ +qn 1,n 1kn 1h() RK ()hkkkkyyii432112261++++=+()()1213243,,2 2,22,iiiiiiiikf t yhhkf tykhhkf tykkf th yk h= =++ =++ =++(4 ) n t dy1dt=f1t,y1,y2, ,yn()dy2dt=f2t,y1,y2, ,yn() dyndt=fnt,y1,y2, ,yn() (adaptive step-size control)

4 (step halving) RK (embedded RK method) RK MATLAB : ode45 , ode23 ode : t y dy/dt ()( ),dy tf t ydt=[t, y] = ode45 (@fun, tspan, y0)functiondydt = fun(t, y)dydt(1) = < Insert a function of t and/or y here. >dydt(2) = ..dydy = dydt ; % t ( )y ( t ) fun ode function m-file tspan ([ ])y0 MATLAB Ex:1.

5 Edit dydt=fun(t,y)dydt(1) = y(1)+y(2)*exp(-t);dydt(2) = -y(1)*y(2)+cos(t);dydt = dydt ;2. ode45 ode23 % MATLAB command window[t,y]= ode45 (@fun, [0 10],[0 1]) %[0 10]: t , [0 1]: x y1=y(:,1);y2=y(:,2);plot(t,y1,'r',t,y2, g')1121221 2(0) 0,(0) 1cos( )tdyyy edtyydyy ytdt =+ == = + MATLAB ode function m-file ode : t , y, p1, p2 dy/dt [t, y] = ode45 (@fun, tspan, y0, options, p1, p2)functiondydt = fun(t, y, p1, p2)dydt(1) = < Insert a function of t and/or y here. >dydt(2) = ..dydt = dydt ; % t y fun ode function m-file tspan ([ ])y0 options ode45 [ ] p1, p2 MATLAB Ex:Step 1.

6 Step 2. 0)0(,1)0(,)sin()cos(= ==+ + yyetytyeytt21221sin( )cos( )ttyyyyyt ee yt y == == 1212221,(0) 1,(0) 0sin( )cos( )ttyyyyyt ee yt y = == = 11yyyy = =MatlabODE MATLAB Stiff ODE ode23tbode23tode23sode15sStiff ODEode113ode23ode45 Non-stiff ODE (stiff system) y(0)=0 dydt= 1000y+3000 2000e ty=3 1000t tExercise ode45 t=0 20 y1 y2 Case 1: a = , b= Case 2.

7 A = , b= 111 212221 ,(0)2,(0) yy ydtyydyb yy ydt = == =+ 11220110001001 = xxxx1(0) 1x=2(0) 0x= (DAE) (Differential-Algebraic Equation, DAE) ()()()()()()111222121211221212, ,,,, ,,, , ,,,0, ,,,0, ,,, 0, ,,,NNnnNnNnNn mNyf x y yyyfx y yyyfx y yyfx y yyfx y yyfx y yy+++====== ()[][]1212,( )( )( )( )000 TNTNn nxyyyfff === = M yFyyFIM i i i ix : Use solver for stiff ODE: ode15s, ode23toptions = odeset( Mass ,M, MassSingular , yes )0(0),1, 2,,iiyyiN== (N = n+m)(Mass matrix)