第二部分常用算法與MATLAB語言編程1．算法的概念2.方程求根3．矩陣求逆4．?dāng)?shù)字積分（歐拉，梯形，龍格庫塔）5．?dāng)?shù)字積分的精度及穩(wěn)定性常用算法

方程求根矩陣求逆數(shù)值積分歐拉法t0t1ty(t)hyy(t0)y(t1)圖2-5-1t0t1tf(t)ff(t0)圖2-5-2h設(shè)一階微分方程

如果t1離t0很近，即h很小，曲線y(t1)可用切線來近似，其切線方程重復(fù)上述近似過程，當(dāng)t=t2時(shí)，則有一般近似公式

如果令，稱為計(jì)算步矩，則

(1)這就是歐拉法數(shù)字積分的遞推計(jì)算公式。如果hn是固定不便的，稱定步長積分；如果hn是變化的，稱變步長積分。由歐拉法數(shù)字積分的遞推計(jì)算公式可看出，只要給出方程的初值（t0,y0）以及相應(yīng)的步距，逐步進(jìn)行遞推就可獲得微分方程的近似數(shù)字解。歐拉法的計(jì)算是十分簡單的，其計(jì)算的截?cái)嗾`差正比于h2

，由此，要獲得高精度解，必須減小步距，但這使得計(jì)算次數(shù)增加，又由于計(jì)算機(jī)的字長有限，h減小得過小，將引起由于舍入誤差累計(jì)而產(chǎn)生的累積誤差過大，所以此方法的精度提高有限，實(shí)際應(yīng)用中較少采用。%%%%歐拉法積分程序x0=0;y0=0;h=0.0001;xm=5;xx=[];x=x0;y=y0;whilex<=xmx=x+h;y=y+fx(x)*h;

%x=x+h;xx=[xx;x,y];endplot(xx(:,1),xx(:,2))holdonx=0:0.1:xm;plot(x,(x.^2)./4,'--r')functiony=fx(x)y=x/2;梯形法（預(yù)報(bào)——校正法）四階龍格-庫塔（Runge-Kutta）法1、Runge-Kutta法推導(dǎo)實(shí)際上，可以由對dy/dx=求導(dǎo)得到高階微商就更復(fù)雜了。為了提高精度，用r階展開計(jì)算公式：一般常用是計(jì)算4階函數(shù)值，得到每步截?cái)嗾`差為的四階Runge-Kutta法其遞推公式為:2、根據(jù)四階龍格-庫塔法的遞推公式：已知系統(tǒng)的狀態(tài)方程為 采用四階龍格-庫塔法進(jìn)行求解和仿真，其求解步驟和方法如下:：1、由，可知；2、仿真算法3、由時(shí)刻的狀態(tài)為，得到

取不斷遞推，便可得到所需時(shí)刻各點(diǎn)的狀態(tài)變量和輸出量。functiony=runk(x1,a,b,c,d)globalts

ymvspttsh=.5*ts;h=[0tsh

tsh

ts];k0=zeros(size(a(:,1)));u1=(ymv(sp)+x1)/2;v=[ymv(sp)u1u1x1];ymv(sp)=x1;sp=sp+1;n=length(k0);y0=[ymv(sp:(sp+n-1))]';k1=a*(y0+h(1)*k0)+b*v(1);k2=a*(y0+h(2)*k1)+b*v(2);k3=a*(y0+h(3)*k2)+b*v(3);k4=a*(y0+h(4)*k3)+b*v(4);y1=y0+[k1,k2,k3,k4]*([1;2;2;1]*ts/6);y=c*y1+d*x1;ymv(sp:(sp+n-1))=y1';sp=sp+n;MATLAB中的ode45（）函數(shù)可實(shí)現(xiàn)四階/五階龍格-庫塔算法，其調(diào)用格式為：[t,y]=ode45('f',tspa,x0)其中：f為定義的常微分方