數(shù)值積分法仿真第一頁，共六十三頁，編輯于2023年，星期三Overview數(shù)值積分方法的原理是什么？病態(tài)系統(tǒng)的特點和仿真算法選??？算法的穩(wěn)定性分析？第二頁，共六十三頁，編輯于2023年，星期三第一節(jié)數(shù)字仿真原理在連續(xù)系統(tǒng)的仿真中，數(shù)值積分法可分為兩大類：單步法：以龍格-庫塔法為代表多步法：以Adams法為代表數(shù)值積分法的要素：基本特性：穩(wěn)定性空間特性：精度時間特性：速度第三頁，共六十三頁，編輯于2023年，星期三數(shù)值積分基本原理連續(xù)系統(tǒng)的仿真，主要是對一階微分方程（組）的求解可見仿真關鍵是對Qm準確，快速的求解第四頁，共六十三頁，編輯于2023年，星期三步長：將時間t離散t(k)(k=1,2,…n)，相鄰兩點的距離為步長，即h=t(k+1)-t(k)步進法：數(shù)值積分法求近似解根據(jù)初始值y0，按照離散的時間序列步進求解。 t0t1t2t3…tn y0y1y2y3…tn計算格式：由y(k)計算出y(k+1)（k=0,1,…,n)的遞推公式。數(shù)值積分基本名詞第五頁，共六十三頁，編輯于2023年，星期三數(shù)值積分的基本性能基本性能數(shù)值積分算法的性能包含：定性特征：穩(wěn)定性時間粒度：計算速度空間粒度：計算精度不同的數(shù)值積分方法的具有不同的穩(wěn)定性。同一個模型采用不同的積分算法和不同的積分步長h，穩(wěn)定性不同。第六頁，共六十三頁，編輯于2023年，星期三計算速度和計算精度各種數(shù)值積分方法的差分方程是對原微分方程的近似逼近，并且因為計算機的字長有限，存在明顯的截斷誤差。這些誤差都和計算步距h密切相關，所以計算步距是影響計算精度、速度和穩(wěn)定性的重要因素。h取得較大，計算時間少，截斷誤差大；h取得較小，截斷誤差就會減小，但在給定時間范圍內，計算次數(shù)必然增加，使誤差積累增加。第七頁，共六十三頁，編輯于2023年，星期三截斷誤差、累計舍入誤差與步長h截斷誤差、累計舍入誤差與步長h關系如圖。圖中可知，兩種誤差對步距的要求是矛盾的，但兩者之和有一個最小值，步距最好能選在最小值。然而，實際要做到這一點是很困難的。一般只能根據(jù)經(jīng)驗確定一個合理步長區(qū)，通常將步長h限制在系統(tǒng)的最小時間常數(shù)數(shù)量級上。第八頁，共六十三頁，編輯于2023年，星期三引理：泰勒級數(shù):如果f(x)在x0點處任意階可導,則在該鄰域內的n階泰勒公式為：第二節(jié)單步法單步數(shù)值積分法的核心就是泰勒級數(shù)的近似。第九頁，共六十三頁，編輯于2023年，星期三2.1一階歐拉法對于一階微分方程故一般的一階歐拉法遞推形式為：第十頁，共六十三頁，編輯于2023年，星期三一階歐拉法圖示第十一頁，共六十三頁，編輯于2023年，星期三2.22階龍格-庫塔對于一階微分方程第十二頁，共六十三頁，編輯于2023年，星期三2階龍格-庫塔第十三頁，共六十三頁，編輯于2023年，星期三2階龍格-庫塔第十四頁，共六十三頁，編輯于2023年，星期三2階龍格-庫塔故一般的二階龍格-庫塔法遞推形式為：2階龍格-庫塔只取到泰勒級數(shù)展開式中y的二階導數(shù)項，略去了三階以上高階導數(shù)項。其截斷誤差正比于步長h3為紀念提出該方法的德國數(shù)學家C.Runge和M.W.Kutta，稱這種計算方法為二階龍格-庫塔法。第十五頁，共六十三頁，編輯于2023年，星期三2階龍格-庫塔圖示第十六頁，共六十三頁，編輯于2023年，星期三比較第十七頁，共六十三頁，編輯于2023年，星期三高階龍格-庫塔(RK-4)一般在計算精度要求較高的情況下，多使用四階龍格-庫塔法。其計算公式為,其截斷誤差正比于步長h5第十八頁，共六十三頁，編輯于2023年，星期三高階龍格-庫塔(RK-4)第十九頁，共六十三頁，編輯于2023年，星期三單步法的特點單步法以上介紹的幾種數(shù)值積分公式，有一個共同的特點，由于采用了泰勒級數(shù)展開，在本次計算中，僅僅用到前一步的計算結果，而不需要利用更前面各步的結果。這類計算方法稱為單步法。單步法運算有下列優(yōu)點：(1)需要存儲的數(shù)據(jù)量少，占用的存儲空間少。(2)給定初值，就可啟動遞推公式進行運算(自啟動計算能力)(3)容易實現(xiàn)變步長第二十頁，共六十三頁，編輯于2023年，星期三第三節(jié)變步長龍格-庫塔法步長控制是實現(xiàn)高精度的仿真算法的手段之一。實現(xiàn)步長控制涉及：局部誤差估計步長控制策略第二十一頁，共六十三頁，編輯于2023年，星期三3.1誤差估計通常設法尋找一個低一階的龍格-庫塔公式，兩者的結果之差可以設為誤差。為減少計算量，Ki通常要求公用。Runge-Kutta-Merson法(RK34)第二十二頁，共六十三頁，編輯于2023年，星期三Runge-Kutta-Fehlberg(RK45)計算公式為5階6級，誤差估計低階公式為4階五級，具有四階誤差估計和五階精度，稱為RK45法。RK45被公認為對非病態(tài)系統(tǒng)仿真的最有效的方法之一。第二十三頁，共六十三頁，編輯于2023年，星期三Runge-Kutta-Fehlberg(RK45)iaibijcic*i1016/13525/21621/41/40033/83/329/326656/128251408/2565412/131932/2197-7200/21977296/219728561/564302197/410451439/216-83680/513-847/4104-9/50-1/561/2-8/272-3544/25661859/4104-11/402/550第二十四頁，共六十三頁，編輯于2023年，星期三RKF-12第二十五頁，共六十三頁，編輯于2023年，星期三RKS-34（1978，Shamping）第二十六頁，共六十三頁，編輯于2023年，星期三3.2步長控制步長控制策略一般分為：1）加倍-減半法2）最優(yōu)步長法第二十七頁，共六十三頁，編輯于2023年，星期三步長控制：加倍-減半法加倍-減半法第二十八頁，共六十三頁，編輯于2023年，星期三步長控制：最優(yōu)步長法第二十九頁，共六十三頁，編輯于2023年，星期三1.2.2步長控制第三十頁，共六十三頁，編輯于2023年，星期三龍格-庫塔方法的一般形式各種龍格-庫塔法的公式都由兩部分組成，一個是上一步結果，另一個是步長乘以各點導數(shù)的加權和。平均斜率第三十一頁，共六十三頁，編輯于2023年，星期三第三節(jié)線性多步法單步法運算基于泰勒級數(shù)展開法，其特點是：(1)需要存儲的數(shù)據(jù)量少，占用的存儲空間少。(2)給定初值(t0,y0)，就可啟動遞推公式進行運算(自啟動計算能力。(3)容易實現(xiàn)變步長積分?？捎行胶庥嬎闼俣群途戎g的矛盾。多步法的基本原理是多項式擬合利用一個多項式取匹配變量的若干已知值和各階導數(shù)。第三十二頁，共六十三頁，編輯于2023年，星期三線性多步法原理第三十三頁，共六十三頁，編輯于2023年，星期三3.1預報公式第三十四頁，共六十三頁，編輯于2023年，星期三3.1預報公式第三十五頁，共六十三頁，編輯于2023年，星期三3.1預報公式第三十六頁，共六十三頁，編輯于2023年，星期三預報舉例第三十七頁，共六十三頁，編輯于2023年，星期三3.2校正公式第三十八頁，共六十三頁，編輯于2023年，星期三3.2校正公式第三十九頁，共六十三頁，編輯于2023年，星期三3.2校正公式第四十頁，共六十三頁，編輯于2023年，星期三預報-校正舉例第四十一頁，共六十三頁，編輯于2023年，星期三預報-校正舉例第四十二頁，共六十三頁，編輯于2023年，星期三3.3Adams公式根據(jù)前面的分析，我們可以將預報和校正公式統(tǒng)一寫成：第四十三頁，共六十三頁，編輯于2023年，星期三顯式Adams系數(shù)第四十四頁，共六十三頁，編輯于2023年，星期三隱式Adams系數(shù)第四十五頁，共六十三頁，編輯于2023年，星期三3.4多步法的特點與單步法相比，相同精度下，使用過去多步信息，計算量小。隱式法的精度高，穩(wěn)定性好，但在計算y(n+k)時需要用到f[y(n+k),t(n+k)]，只能采用迭代法計算。缺點之一是不能自啟動，需用單步法計算初始值才能啟動計算。第四十六頁，共六十三頁，編輯于2023年，星期三第四節(jié)積分算法的穩(wěn)定性穩(wěn)定的系統(tǒng)采用不同的積分算法，其穩(wěn)定性不同穩(wěn)定性的測試公式為：當?shù)谒氖唔?，共六十三頁，編輯?023年，星期三4.1一階Adams法的穩(wěn)定性分析第四十八頁，共六十三頁，編輯于2023年，星期三一階Adams法的穩(wěn)定性分析第四十九頁，共六十三頁，編輯于2023年，星期三4.2一般算法的穩(wěn)定性分析根據(jù)上例可得數(shù)值積分方法穩(wěn)定域的一般方法。 設系統(tǒng)測試方程為： 而數(shù)值積分公式為： 只有當時，算法穩(wěn)定。各種數(shù)值積分算法的穩(wěn)定域參見書P96圖3.9第五十頁，共六十三頁，編輯于2023年，星期三主要算法的穩(wěn)定性一階、二階Admas法為恒穩(wěn)算法，其他算法條件穩(wěn)定。除恒穩(wěn)法外，其他算法的步長h必須限制在最小時間的數(shù)量級對龍格-庫塔法，階次k增大，穩(wěn)定域略微增大。對Admas法，階次k增大，穩(wěn)定域反而縮小。第五十一頁，共六十三頁，編輯于2023年，星期三第五節(jié)Matlab實現(xiàn)ODE(OrdinaryDifferentaialequation)解法模型描述算法描述算法仿真第五十二頁，共六十三頁，編輯于2023年，星期三微分方程模型描述Lorenz曲線filename:mdLorenz.m

functiondx=mdLorenz(t,x)dx=[-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3)];end第五十三頁，共六十三頁，編輯于2023年，星期三數(shù)值積分算法描述matlab中的數(shù)值積分算法函數(shù)的格式如下:function[tout,yout]=solver(ModelName,tspan,x0,option)第五十四頁，共六十三頁，編輯于2023年，星期三數(shù)值積分算法描述%一階Euler算法,filename:svEulerfunction[tout,yout]=svEuler(odeFcn,tspan,y0)t0=tspan(1);t1=tspan(2);iflength(tspan)<3,h=(t1-t0)/1000;elseh=tspan(3);tout=[t0:h:t1]';N=length(y0);M=length(tout)-1;tout=[t0:h:t1]';yout=[y0';zeros(M,N)];fori=1:Mk1=h*feval(odeFcn,tout(i),y0);y0=y0+k1;yout(i+1,:)=y0';endend第五十五頁，共六十三頁，編輯于2023年，星期三數(shù)值積分算法描述function[tout,yout]=svRungeKutta4(odeFcn,tspan,y0)t0=tspan(1);t1=tspan(2);iflength(tspan)<3,h=(t1-t0)/1000;elseh=tspan(3);tout=[t0:h:t1]';N=length(y0);M=length(tout)-1;tout=[t0:h:t1]';yout=[y0';zeros(M,N)];fori=1:Mk1=h*feval(odeFcn,tout(i),y0);k2=h*feval(odeFcn,tout(i)+h/2,y0+0.5*k1);k3=h*feval(odeFcn,tout(i)+h/2,y0+0.5*k2);k4=h*feval(odeFcn,tout(i)+h,y0+k3);y0=y0+(k1+2*k2+2*k3+k4)/6;yout(i+1,:)=y0';endend第五十六頁，共六十三頁，編輯于2023年，星期三微分方程模型描述在Matlab文件中調用方法為:[t,y]=svEuler(@eqLorenz,[0,100],x0);第五十七頁，共六十三頁，編輯于2023年，星期三ode45ODE45Solvenon-stiffdifferentialequations,mediumordermethod.[TOUT,YOUT]=ODE45(ODEFUN,TSPAN,Y0)withTSPAN=[T0TFINAL]integratesthesystemofdifferentialequationsy'=f(t,y)fromtimeT0toTFINALwithinitialconditionsY0.ODEFUNisafunctionhandle.ForascalarTandavectorY,ODEFUN(T,Y)mustreturnacolumnvectorcorrespondingtof(t,y).EachrowinthesolutionarrayYOUTcorrespondstoatimereturnedinthecolumnvectorTOUT.ToobtainsolutionsatspecifictimesT0,T1,...,TFINAL(allincreasingoralldecreasing),useTSPAN=[T0T1...TFINAL].第五十八頁，共六十三頁，編輯于2023年，星期三ode45ODE45Solvenon-stiffdifferentialequations,mediumordermethod.

[TOUT,YOUT]=ODE45(ODEFUN,TSPAN,Y0,OPTIONS)solvesasabovewithdefaultintegrationpropertiesreplacedbyvaluesinOPTIONS,anargumentcreatedwiththeODESETfunction.SeeODESETfordetails.Commonlyusedoptionsarescalarrelativeerrortolerance'RelTol'(1e-3bydefault)andvectorofabsoluteerrortolerances'AbsTol'(allcomponents1e-6bydefault).Ifcertaincomponentsofthesolutionmustbenon-negative,useODESETtosetthe'NonNegative'propertytotheindicesofthese第五十九頁，共六十三頁，編輯于2023年，星期三ode45ODE45Solvenon-stiffdifferentialequations,mediumordermethod.

ODE45cansolveproblemsM(t,y)*y'=f(t,y)withmassmatrixMthatisnonsingular.UseODESETtosetthe'Mass'propertytoafunctionhandleMASSifMASS(T,Y)returnsthevalueofthemassmatrix.Ifthemassmatrixisconstant,thematrixcanbeusedasthevalueofthe'Mass'option.IfthemassmatrixdoesnotdependonthestatevariableYandthefunctionMASSistobecalledwithoneinputargumentT,set'MStateDependence'to'none'.ODE15SandODE23Tcansolveproblemswithsingularmassmatrices.第六十頁，共六十三頁，編輯于2023年，星期三ode45

ODE45Solvenon-stiffdifferentialequations,mediumordermethod.Example[t,y]=ode45(@vdp1,[020],[20]);plot(t,y(:,1));solvesthesystemy'=vdp1(t,y),usingthedefaultrelativeerrortolerance1e-3andthedefaultabsolutetoleranceof1e-6foreachcomponent,andplotsthefirstcomponentofthesolution.

ClasssupportforinputsTSPAN,Y0,andtheresultofODEFUN(T,Y):float:double,sin