4.3

常微分方程數值解及其MATLAB實現方程無解析解數值解解的性質

考察方程隨（幾）個參數的變化，以及解是如何變化的引入新變量替換

首先把高階微分方程（組）化簡為一階微分方程組m+n維的一階微分方程組解例將Vanderpol方程化為一階微分方程組說明：常微分方程數值解的提法一般是針對一個一階方程（組）和初始條件來說的，比如有方程：數值解求方程數值解的常規方法Euler法變步長Euler法Runge-Kutta法Runge-Kutta-Felhberg法等

本書主要介紹Euler法和Runge-Kutta法Euler方法在小區間上用差商來代替方程中導數，Euler方法可分為向前和向后Euler方法.如果

中的t在上選取區間的左端點，則有向前Euler公式：1.基本思想：t0t1t2近似解xth2.幾何意義：解析解h大小影響計算速度與精度迭代次數的增加會帶來

較多的累積誤差

向前Euler公式的

精度并不很高提高精度的辦法:使用變步長改進Euler公式等1階精度1、向前Euler法（一下修改？？？）若f(t,x(t))中的t在上選取區間的左端點tn，則由方程可得（4-18）此公式稱為向前Euler公式4.3.1數值解的Euler法t0t1t2近似解xt圖4.3向前Euler方法求近似解2、向后Euler法若中的

在上選取區間的右端點，則由方程可得（4-19）此公式稱為向后Euler公式4.3.1數值解的Euler法圖4.4向后Euler方法求近似解結合向前和向后Euler公式,可以得到更高精度的梯形公式：隱式公式，需要迭代方法求解得到以下公式進行迭代：梯形公式3．改進的Euler法或例4.9分別用向前Euler公式和改進的Euler公式求方程（4-26）在區間上步長為h=0.1的數值解。解：方程（4-26）的向前Euler公式的形式：改進Euler公式：（4-28）（4-27）方程精確解：（4-29）（4-27）（4-28）（4-29）（4-27）（4-28）（4-29）0.10.90000.90500.90480.60.53140.54940.54880.20.81000.81900.81870.70.47830.49720.49660.30.72900.74120.74080.80.43050.45000.44930.40.65610.67080.67030.90.38740.40720.40660.50.59050.60710.606510.34870.36850.3679表4.1方程（4-26）的三種數值解4.3.2數值解的Runge-Kutta法Euler方法是使用差商代替導數，并進行迭代的方法。事實上，根據微分中值定理可以得到（4-30）結合應有（4-31）

Runge-Kutta方法Euler方法的基本思想，即差商代替導數，很自然的想法是在區間內多取幾個點，將它們的斜率加權平均作為導數近似值，這就是Runge-Kutta方法的思想.1.基本思想:2.2階R-K公式:3.4階R-K公式:公式復雜，如何編程實現？？？在Matlab中，應用R-K法求解方程數值解的函數：[t,x]=ode45（或ode23）(Fun,[t_0,t_f],x_0)[t,x]=ode45(Fun,[t_0,t_f],x_0,options)[t,x]=ode45(Fun,[t_0,t_f],x_0,options,p_1,…）Matlab函數可以由指定的M-文件給出;[t_0,t_f]為自變量取值區間;options可給出誤差限(缺省時相對誤差1e-3,絕對誤差1e-6);具體形式為：options=odeset(‘reltol’,rt,’abstol’,at).

4.3.3數值解的MATLAB實現了解解例質量m=1kg，線長l=1m，重力加速度g=9.8m/s2，初始角度θ0=π/12，初始速度0，阻力系數λ=0.1，求數學擺方程的數值解。4.3.3數值解的MATLAB實現數學擺模型將方程化為一階方程組已知質點M：首先編寫M-文件lorenzeq.m

functionxdot=lorenzeq(t,x)xdot=[-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3)];解例畫出Lorenz系統數值解曲線4.3.3數值解的MATLAB實現在命令窗口中運行

>>t_final=100;x0=[0;0;1e-10];[t,x]=ode45('lorenzeq',[0,t_final],x0);在左側workspace中看到輸出的數值解如下：4.3.3數值解的MATLAB實現續上例應用plot（t，x）畫圖得到三條曲線如下圖：續上例4.3.3數值解的MATLAB實現t-x1曲線t-x2曲線t-x3曲線也可繪制三維曲線（x1-x2-x3曲線）如下：>>figure;plot3(x(:,1),x(:,2),x(:,3));得到了Loren系統具有混沌性質的解曲線.應用comet3(x(:,1),x(:,2),x(:,3));可看動態圖像，

可以看到Lorenz系統具有混沌性的解。續上例4.3.3數值解的MATLAB實現解例針對上述方程，也可在編程時將參數設計在程序中，通過主窗口下賦不同值，求不同參數下微分方程組

解，此時[]用于占位，不可省略。

首先編寫m-文件lorenzleq.mfunctionxdot=lorenzleq(t,x,flag,beta,rho,sigma)xdot=[-beta*x(1)+x(2)*x(3);-rho*x(2)+rho*x(3);-x(1)*x(2)+sigma*x(2)-x(3)];4.3.3數值解的MATLAB實現

在窗口中可對不同變量賦值，

>>t_final=100;x0=[0;0;1e-10];b1=8/3;r1=10;s1=28;（注意：這里的參數名不要求與M-文件中一致）

[t,x]=ode45('lorenzleq',[0,t_final],x0,[],b1,r1,s1);plot(t,x（:,1）)；plot(t,x（:,2）)plot(t,x（:,3）)

得到與上例相同的圖像，變化參數可以不用修改程序得到不同的方程的數值解4.3.3數值解的MATLAB實現例如改變參數值為b1=8;r1=5;s1=5.

重新運行[t,x]=