1、Matlab解微分方程,1,除了上述的已知ODE外，還須有起始條件y0=y(x0)才能解方程式，即是在x=x0時，y(x)=y0。上述各個方程式 的解析解 (analytical solution) 如下：,2,阮奇-庫達 (Runge-Kutta) 方法是最通用的解 ODE 的方法，它可以依計算精確度的要求有低階到高階的各個 計算式，,3,MATLAB應(yīng)用阮奇-庫達法的函數(shù)有ode23, ode45，其中ode23是同時以二階及三階阮奇-庫達法求解，而ode45 則是以四階及五階阮奇-庫達法求解。其語法為ode23(dy,x0,xn,y0)，其中 dy 為ODE中的等式右邊的函數(shù)（如之 前介

2、紹的），x0, xn 是要解ODE的區(qū)間 x0, xn 的二個端點，y0是起始值 (y0=y(x0)。而ode45的語法 與ode23相同。,4,例一、要在區(qū)間 2, 4 解以下的 ODE： % m-function, g1.m function dy=g1(x,y) dy=3*x.2; x,num_y=ode23(g1,2,4,0.5); anl_y=x.3-7.5; plot(x,num_y, b:,x,anl_y, r-) title(Solution of g1) xlabel(x), ylabel(y=f(x), grid,Y = 3x2,5,例二、要在區(qū)間 0, 5 解以下的 OD

3、E： % m-function, g2.m function dy=g2(x,y) dy=-0.131*y; x,num_y=ode23(g2,0,5,4); anl_y=4*exp(-0.131*x); plot(x,num_y,x,anl_y,o) title(Solution of g2) xlabel(x), ylabel(y=f(x), grid,Y = -0.13y,6,例三、要在區(qū)間 0, 2 解以下的 ODE： % m-function, g3.m function dy=g3(x,y) dy=2*x*cos(y)2; x,num_y=ode23(g3,0,2,pi/4); a

4、nl_y=atan(x.*x+1); plot(x,num_y, b:,x,anl_y, r-) title(Solution of g3) xlabel(x), ylabel(y=f(x), grid,Y = 2x cos2y,7,例四、要在區(qū)間 0, 3 解以下的 ODE： % m-function, g4.m function dy=g4(x,y) dy=3*y+exp(2*x); x,num_y=ode23(g4,0,3,3); anl_y=4*exp(3*x)-exp(2*x); plot(x,num_y,x,anl_y,o) title(Solution of g4) xlabel

5、(x), ylabel(y=f(x), gri,Y = 3y + e2x,8,如果將上述方法改以 ode45 計算，你可能無法察覺出其與ode23的解之間的差異，那是因為我們選的 ODE 代表的函數(shù)分布變化平緩，所以高階方法就顯現(xiàn)不出其優(yōu)點。不過以二者在計算的誤差上做比較，ode45 的誤差量級會比 ode23要小。以下是幾個例子：,9,% m-function, g1.m function dy=g1(x,y) dy=3*x.2; % m-file, odes1.m ; % Solve an ode using ode23 and ode45 clg x1,num_y1=ode23(g1,2

6、,4,0.5); anl_y1=x1.3-7.5; error_1=abs(anl_y1-num_y1)./abs(anl_y1); % ode23 的計算誤差,10,x2,num_y2=ode45(g1,2,4,0.5); anl_y2=x2.3-7.5; % 注意 x2 個數(shù)與 x1 不一定相同 error_2=abs(anl_y2-num_y2)./abs(anl_y2); % ode45 的計算誤差hold on subplot(2,2,1) plot(x1,num_y1,x1,anl_y1,o) title(ODE23 solution), ylabel(y),11,subplot(

7、2,2,2) plot(x1,error_y1) % 注意二種方法的誤差 title(ODE23 error), ylabel(y) % ode23 的誤差的量級為 1.e-16 subplot(2,2,3) plot(x2,num_y2,x2,anl_y2,o) title(ODE45 solution), ylabel(y) subplot(2,2,4) plot(x1,error_y2) title(ODE45 error), ylabel(y) % ode45 的解沒有誤差 hold off,12,另一個例子： % m-function, g5.m function dy=g5(x,y

8、) dy=-y+2*cos(x); % m-file, odes1.m % Solve an ode using ode23 and ode45 clg,13,x1,num_y1=ode23(g5,0,5,1); anl_y1=sin(x1)+cos(x1); error_1=abs(anl_y1-num_y1)./abs(anl_y1); x2,num_y2=ode45(g5,0,5,1); anl_y2=sin(x2)+cos(x2); error_2=abs(anl_y2-num_y2)./abs(anl_y2); hold on,14,subplot(2,2,1) plot(x1,nu

9、m_y1,x1,anl_y1,o) title(ODE23 solution), ylabel(y) subplot(2,2,2) plot(x1,error_1) % 注意二種方法的誤差 title(ODE23 error), ylabel(y) % ode23 的誤差的量級為 1.e-4,15,subplot(2,2,3) plot(x2,num_y2,x2,anl_y2,o) title(ODE45 solution), ylabel(y) subplot(2,2,4) plot(x2,error_2) title(ODE45 error), ylabel(y) % ode45 的誤差的

10、量級為 1.e-6 hold off,16,高階常微分方程式,一個高階常微分方程式可以利用變數(shù)改變(change of variables) 方式改寫成一組聯(lián)立的一階常微分方程式。 例如以下的 n 階方程式,17,以下即是解二階 ODE 的解法: u_prime(1) = u(1)*(1-u(2)2) - u(2); u_prime(2) = u(1); u(1)=y=y*(1-y2)-y; y=u(1); function u_prime =eqns2(x,u) u_prime = u(1)*(1-u(2)2) - u(2); u(1); initial = 0 0.25; x,num_y = ode23(eqns2,0,20,initial); subplot(2,1,1), plot(x,num_y(:,1) title(1st derivative of y