第3節(jié)常微分方程數(shù)值解及其matlab實現(xiàn)

4.3

常微分方程數(shù)值解及其MATLAB實現(xiàn)方程無解析解數(shù)值解解的性質(zhì)

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

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

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

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

較多的累積誤差

向前Euler公式的

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

在上選取區(qū)間的右端點，則由方程可得（4-19）此公式稱為向后Euler公式4.3.1數(shù)值解的Euler法圖4.4向后Euler方法求近似解結(jié)合向前和向后Euler公式,可以得到更高精度的梯形公式：隱式公式，需要迭代方法求解得到以下公式進行迭代：梯形公式3．改進的Euler法或例4.9分別用向前Euler公式和改進的Euler公式求方程（4-26）在區(qū)間上步長為h=0.1的數(shù)值解。解：方程（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）的三種數(shù)值解4.3.2數(shù)值解的Runge-Kutta法Euler方法是使用差商代替導(dǎo)數(shù)，并進行迭代的方法。事實上，根據(jù)微分中值定理可以得到（4-30）結(jié)合應(yīng)有（4-31）

Runge-Kutta方法Euler方法的基本思想，即差商代替導(dǎo)數(shù)，很自然的想法是在區(qū)間內(nèi)多取幾個點，將它們的斜率加權(quán)平均作為導(dǎo)數(shù)近似值，這就是Runge-Kutta方法的思想.1.基本思想:2.2階R-K公式:3.4階R-K公式:公式復(fù)雜，如何編程實現(xiàn)？？？在Matlab中，應(yīng)用R-K法求解方程數(shù)值解的函數(shù)：[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函數(shù)可以由指定的M-文件給出;[t_0,t_f]為自變量取值區(qū)間;options可給出誤差限(缺省時相對誤差1e-3,絕對誤差1e-6);具體形式為：options=odeset(‘reltol’,rt,’abstol’,at).

4.3.3數(shù)值解的MATLAB實現(xiàn)了解解例質(zhì)量m=1kg，線長l=1m，重力加速度g=9.8m/s2，初始角度θ0=π/12，初始速度0，阻力系數(shù)λ=0.1，求數(shù)學擺方程的數(shù)值解。4.3.3數(shù)值解的MATLAB實現(xiàn)數(shù)學擺模型將方程化為一階方程組已知質(zhì)點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系統(tǒng)數(shù)值解曲線4.3.3數(shù)值解的MATLAB實現(xiàn)在命令窗口中運行

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

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

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

首先編寫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數(shù)值解的MATLAB實現(xiàn)

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

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

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

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

重新運行[t,x]=