資源描述:
《隱式euler求解一階常微分方程》由會員上傳分享,免費在線閱讀,更多相關(guān)內(nèi)容在行業(yè)資料-天天文庫。
1、《MATLAB語言及應(yīng)用》大作業(yè)姓名:學(xué)號:學(xué)院:班級:題目編號:2013年10月13隱式Euler求解一階常微分方程。一、隱式Euler的數(shù)學(xué)理論微分方程的本質(zhì)特征是方程中含有導(dǎo)數(shù)項,數(shù)值解法的第一步就是設(shè)法消除其導(dǎo)數(shù)值,這個過程稱為離散化。實現(xiàn)離散化的基本途徑是用向前差商來近似代替導(dǎo)數(shù),這就是歐拉算法實現(xiàn)的依據(jù)。歐拉(Euler)算法是數(shù)值求解中最基本、最簡單的方法,但其求解精度較低,一般不在工程中單獨進(jìn)行運算。所謂數(shù)值求解,就是求問題的解y(x)在一系列點上的值y(xi)的近似值yi。對于常微分方程:dy
2、/dx=f(x,y),x∈[a,b]y(a)=y0可以將區(qū)間[a,b]分成n段,那么方程在第xi點有y'(xi)=f(xi,y(xi)),再用向前差商近似代替導(dǎo)數(shù)則為:(y(xi+1)-y(xi))/h=f(xi,y(xi)),在這里,h是步長,即相鄰兩個結(jié)點間的距離。因此可以根據(jù)xi點和yi點的數(shù)值計算出yi+1來:yi+1=yi+h*f(xi,yi),i=0,1,2,L這就是歐拉公式,若初值yi+1是已知的,則可依據(jù)上式逐步算出數(shù)值解y1,y2,L。為簡化分析,人們常在yi為準(zhǔn)確即yi=y(xi)的前提下估
3、計誤差y(xi+1)-yi+1,這種誤差稱為局部截斷誤差。如果一種數(shù)值方法的局部截斷誤差為O(h^(p+1)),則稱它的精度是p階的,或稱之為p階方法。歐拉格式的局部截斷誤差為O(h^2),由此可知歐拉格式僅為一階方法。二、隱式Euler的算法和流程圖算法:用向后差商[y(x(i+1))-y(x(i))]/h替代方程y’(x(i+1))=g[x(i+1),y(x(i+1))]中的導(dǎo)數(shù)項y’(x(i+1)),再離散化,即可導(dǎo)出下列格式y(tǒng)(i+1)=y(i)+hg(x(i+1),y(i+1)),i=0,1,2。。。
4、因此是隱式公式,一般要用迭代法求解,迭代公式通式為:y^(0)(i+1)=y(n)+hg(x(i),y(i));y^(k+1)(i+1)=y(i)+hg(x(i+1),y^(k)(i+1))(k=0,1,2….)流程圖:三、隱式Euler的Matlab實現(xiàn)建立Euler_2.m文件:functionE2=Euler_2(fun,x0,y0,xN,N)%向后Euler公式,其中,%fun為一階微分方程的函數(shù)%x0,y0為初始條件%xN為取值范圍的一個端點%h為區(qū)間步長%N為區(qū)間的個數(shù)%x為xN構(gòu)成的向量%y為yN
5、構(gòu)成的向量x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;forn=1:N%用迭代法求y(n+1)x(n+1)=x(n)+h;z0=y(n)+h*feval(fun,x(n),y(n));fork=1:3z1=y(n)+h*feval(fun,x(n+1),z0);ifabs(z1-z0)<1e-3break;endz0=z1;endy(n+1)=z1;endT=[x',y']四、隱式Euler的算例實現(xiàn)例:d(y)/d(x)=-y,y(0)=
6、1,(0≤x≤1)的近似解,并與精確值比較。建立f.m文件functionz=f(x,y)z=-y;>>Euler_2('f',0,1,1,10)T=01.00000.10000.90910.20000.82640.30000.75120.40000.68280.50000.62070.60000.56420.70000.51290.80000.46620.90000.42381.00000.3852>>