個人資料整理 僅限學(xué)習(xí)使用數(shù)學(xué)與信息科學(xué)系實(shí)訓(xùn)報常微分方程初值問題數(shù)值解法及其MATLAB實(shí)現(xiàn)姓名：學(xué)號：專業(yè)：信息與計算科學(xué)年級：2018級指導(dǎo)教師：朱耀生完成時間：2018年11月25—2018年12月6實(shí)驗(yàn)?zāi)康?1/16個人資料整理 僅限學(xué)習(xí)使用(1>研究滿足給定方程的可微函數(shù)的數(shù)值方法(2>培養(yǎng)matlab編程和上機(jī)調(diào)試能力實(shí)驗(yàn)基本原理和內(nèi)容 :根據(jù)給定的初始條件，確定常微分方程惟一解的問題叫常微分方程初值問題。大多數(shù)實(shí)際問題難以求得解讀解，必須將微分問題離散化 ， 用 數(shù) 值 方 法 求 其 近 似 解 。一階常微分方程的初值問題的提法是，求出函數(shù) ，使?jié)M足條件(1>b5E2RGbCAP利用數(shù)值方法解問題 (1>時，通常假定解存在且惟一，解函數(shù)及右端函數(shù)具有所需的光滑程度。數(shù)值解法的基本思想是 :先取自變量一系列離散點(diǎn),把微分問題(1>離散化，求出離散問題的數(shù)值解 ,并以此作為微分問題解的近似。例如取步長 >0,以剖分區(qū)間[,]，令 =+ ,把微分方程離散化成一個差分方程。以 (>表微分方程初值問題的解,以 表差分問題的解， 就是近似解的誤差，稱為全局誤差。因此，設(shè)計各種離散化模型 ,求出近似解,估計誤差以及研究數(shù)值方法的穩(wěn)定性和收斂性等構(gòu)成了數(shù)值解法的基本內(nèi)容 p1EanqFDPw1.歐拉法和后退歐拉法原理:2/16個人資料整理 僅限學(xué)習(xí)使用function [x,y]=meuler(df,xspan,y0,h>%用途：改進(jìn)歐拉公式解常微分方程 y'=f(x,y>,y(x0>=y0%格式：[xy]=meuler(df,a,b,y0,h>df 為函數(shù)f(x,y>,xspan 為求解DXDiTa9E3d%區(qū)間[x0,xn],y0 為初值y(x0>,h為步長,[xy] 返回節(jié)點(diǎn)和數(shù)值解矩陣x=xspan(1>:h:xspan(2> 。y(1>=y0。forn=1:(length(x>-1>k1=feval(df,x(n>,y(n>> 。y(n+1>=y(n>+h*k1 。k2=feval(df,x(n+1>,y(n+1>> 。y(n+1>=y(n>+h*(k1+k2>/2 。end例題給定的初值問題y′=-y+x+2，0=<x=<1Y(0>=-1,取精確解y(x>=exp(-x>+x后退歐拉法，步長h=0.003,h=0.1求在節(jié)點(diǎn)k=1+0.1k(k=1,2,3 10>處的數(shù)值解若>>h=0.1。3/16個人資料整理 僅限學(xué)習(xí)使用y=-1。x=1。fori=1:20k1=h*oulei_wf(x,y>。k2=h*oulei_wf(x+h,y+k1>。y=y+0.5*k1+0.5*k2x=x+h。z=ouleij_q(x>t=y-zendz=4,z=3.7000,y=-0.6150z=1.4329,z=1.4329,t=-2.0479z=3.7150,z=3.4435,y=-0.2571z=1.5012,z=1.5012,t=-1.7583z=3.4571,z=3.2114,y=0.0763z=1.5725,z=1.5725,t=-1.4962z=3.2237,z=3.0013,y=0.3876z=1.6466,z=1.6466,t=-1.2590z=3.0124,z=2.8112,y=0.6788z=1.7231,z=1.7231,t=-1.0444z=2.8212,z=2.6391,y=0.9518z=1.8019,z=1.8019,t=-0.85014/16個人資料整理 僅限學(xué)習(xí)使用z=2.6482,z=2.4834,y=1.2084z=1.8827,z=1.8827,t=-0.6743z=2.4916,z=2.3425,y=1.4501z=1.9653,z=1.9653,t=-0.5152實(shí)驗(yàn)結(jié)果的分析與研究1.對于歐拉法，步長越小，精度越高，而產(chǎn)生的誤差越小。總對于歐拉法，步長越小，精度越高，而產(chǎn)生的誤差越小。對于歐拉法 體來說，歐拉法的優(yōu)點(diǎn)是形式簡單 計算方便，是形式簡單，體來說，歐拉法的優(yōu)點(diǎn)是形式簡單，計算方便，缺點(diǎn)是總的運(yùn)算精度比較低。 的增大，誤差值也越來越大。算精度比較低。而且隨著 x 的增大，誤差值也越來越大。根據(jù)歐拉公式的截斷誤差計算，歐拉法是一階方法。歐拉公式的截斷誤差計算，歐拉法是一階方法。 RTCrpUDGiT2.對于改進(jìn)歐拉法， 對于改進(jìn)歐拉法， 其基本特征與歐拉法相似，也是步長越小， 對于改進(jìn)歐拉法 其基本特征與歐拉法相似，也是步長越小，精度越高，誤差越小。優(yōu)點(diǎn)是精度相對歐拉法來說較高，精度越高，誤差越小。優(yōu)點(diǎn)是精度相對歐拉法來說較高，截斷誤差為O(h^3>，缺點(diǎn)是比歐拉法計算量大。根據(jù)歐拉改進(jìn)法 ，缺點(diǎn)是比歐5/16個人資料整理 僅限學(xué)習(xí)使用拉法計算量大。 公式的截斷誤差計算，歐拉改進(jìn)法是二階方法。公式的截斷誤差計算，歐拉改進(jìn)法是二階方法。 5PCzVD7HxA2．龍格-庫塔方法原理：對于龍哥庫塔三階：利用Yn+1=Yn+h/6(K1+4K2+K3>K1=f<Xn,Yn）K2=f<Xn+1/2*h,Yn+(h/2>*K1）K3=f<Xn+h,Yn-K1+2*h*k2）若>>clearallh=0.1。y=-1。x=1。fori=1:25。k1=oulei_wf(x,y>。k2=oulei_wf(x+h/2,y+k1*(h/2>> 。k3=oulei_wf(x+h,y-k1*h+k2*(2*h>> 。y=y+(k1+4*k2+k3>*(h/6>。x=x+h。z=ouleij_q(x>。t=abs(y-z>。6/16個人資料整理 僅限學(xué)習(xí)使用A=[xyzt]endz=4,z=3.8500,z=3.7300,A=1.1000-0.61451.43292.0474z=3.7145z=3.5788z=3.4702z=1.5012A=1.2000-0.25621.50121.7574z=3.4562z=3.3334z=3.2351z=1.5725A=1.30000.07761.57251.4950z=3.2224z=3.1113z=3.0224z=1.6466A=1.40000.38911.6466 1.2575z=3.0109z=2.9104z=2.8299z=1.7231A=1.50000.68041.72311.0427z=2.8196z=2.7286z=2.6558z=1.8019A=1.60000.95361.80190.8483z=2.6464z=2.5641z=2.4982z=1.8827A=1.70001.21031.8827 0.6724z=2.4897z=2.4152z=2.3556z=1.9653A=1.80001.45211.96530.5132z=2.3479z=.2805z=2.2266z=2.0496A=1.90001.68032.04960.3692例題17/16個人資料整理 僅限學(xué)習(xí)使用給定的初值問題y′=-y+x+2，0=<x=<1Y(0>=-1,取精確解y(x>=exp(-x>+x龍格庫塔法三階方法,h=0.1求在節(jié)點(diǎn)k=1+0.1k(k=1,2,3 10>處的數(shù)值解例題2給定的初值問題y′=-y+x+2，0=<x=<1Y(0>=-1,取精確解y(x>=exp(-x>+x龍格庫塔四階法,h=0.1求在節(jié)點(diǎn)k=1+0.1k(k=1,2,3 10>處的數(shù)值解原理：function [x,y]=m4rkodes(df,xspan,y0,h>%用途:4階經(jīng)典龍格庫塔公式解常微分方程組 y'=f(x,y>,y(x0>=y0%格式:[x,y]=m4rkodes(df,xspan,y0,h>%df為向量函數(shù)f(x,y>表達(dá)式,xspan為求解區(qū)間[x0,xn],%y0為初值向量,h為步長,x返回節(jié)點(diǎn),y返回數(shù)值解向量x=xspan(1>:h:xspan(2> 。y=zeros(length(y0>,length(x>> 。y(:,1>=y0(:> 。8/16個人資料整理 僅限學(xué)習(xí)使用forn=1:(length(x>-1>k1=feval(df,x(n>,y(:,n>> 。k2=feval(df,x(n>+h/2,y(:,n>+h/2*k1> 。k3=feval(df,x(n>+h/2,y(:,n>+h/2*k2> 。k4=feval(df,x(n+1>,y(:,n>+h*k3> 。y(:,n+1>=y(:,n>+h*(k1+2*k2+2*k3+k4>/6 。endx=x'。y=y'。若>>h=0.1。y=-1。x=1。fori=1:20k1=h*oulei_wf(x,y>。k2=oulei_wf(x+h/2,y+k1*(h/2>>。k3=oulei_wf(x+h/2,y+k2*(h/2>>。k4=oulei_wf(x+h,y+h*k3>。y=y+(k1+2*k2+2*k3+k4>*(h/6> 。x=x+h。z=ouleij_q(x>。t=abs(y-z>。A=[xyzt]end9/16個人資料整理 僅限學(xué)習(xí)使用z=4z=3.8500z=3.8575z=3.7142z=1.4329A=1.1000-0.61451.43292.0407z=3.7145z=3.5788z=3.5856z=3.4560z=1.5012A=1.2000-0.25621.50121.7574z=3.4562z=3.3334z=3.3395z=3.2222z=1.5725A=1.30000.07751.57251.4950z=3.2225z=3.1113z=3.1169z=3.0108z=1.6466z=1.40000.38901.64661.2576z=3.0110z=2.9104z=2.9154z=2.8194z=1.7231A=1.50000.68041.72311.0427z=2.8196z=2.7286z=2.7332z=2.6463z=1.8019A=1.60000.95361.80190.8483z=2.6464z=2.5641z=2.5682z=2.4896z=1.8827A=1.70001.21021.88270.6724z=2.4898z=2.4153z=2.4190z=2.3479z=1.9653A=1.80001.45201.96530.5133z=2.3480z=2.2806z=2.2840z=2.2196z=2.0496A=1.90001.68032.04960.3693z=2.2197z=2.1587z=2.1618z=2.1035z=2.1353A=2.00001.89642.13530.2390z=2.1036z=2.0485z=2.0512z=1.9985z=2.2225A=2.10002.10142.22250.121110/16個人資料整理 僅限學(xué)習(xí)使用z=1.9986z=1.9487z=1.9512z=1.9035z=2.3108A=2.20002.29642.31080.0144結(jié)果的分析與總結(jié)對于龍格—庫塔法，優(yōu)點(diǎn)是精度更高，同樣的步長下精度比對于龍格—對于龍格 庫塔法，優(yōu)點(diǎn)是精度更高， 歐拉法高的多，歐拉法高的多，誤差更小 ，截斷誤差為 O(h^5>，故龍格—庫，故龍格—塔法是四階方法。缺點(diǎn)是每步都要計算四次微分值。塔法是四階方法。缺點(diǎn)是每步都要計算四次微分值。 jLBHrnAILg3.亞當(dāng)斯法%程序function[x,y]=m4adams(df,xspan,y0,h>%用途：4階亞當(dāng)斯預(yù)報-校正格式解常微分方程 y'=f(x,y>,y(x0>=y0%格式：[x,y]=m4adams(df,xspan,y0,h>%df為函數(shù)f(x,y>表達(dá)式,xspan為求解區(qū)間[x0,xn],%y0為初值,h為步長,x返回節(jié)點(diǎn),y返回數(shù)值解x=xspan(1>:h:xspan(2> 。[x1,y]=m4rungekutta(df,[x(1>,x(4>],y0,h> 。forn=4:(length(x>-1>p=y(n>+h/24*(55*feval(df,x(n>,y(n>>-59*feval(df,x(n-1>,y(n-1>> ... xHAQX74J0X11/16個人資料整理 僅限學(xué)習(xí)使用+37*feval(df,x(n-2>,y(n-2>>-9*feval(df,x(n-3>,y(n-3>>>。LDAYtRyKfEy(n+1>=y(n>+h/24*(feval(df,x(n-2>,y(n-2>>-5*feval(df,x(n-1>,y(n-1>>...Zzz6ZB2Ltk+19*feval(df,x(n>,y(n>>+9*feval(df,x(n+1>,p>> 。end綜合實(shí)例：市場價格模型對于純粹的市場經(jīng)濟(jì)來說 ,商品市場價格取決于市場供需之間的關(guān)系,市場價格能促使商品的供給與需求相等 (這樣的價格稱為(靜態(tài)>均衡價格>.也就是說,如果不考慮商品價格形成的動態(tài)過程 ,那么商品的市場價格應(yīng)能保證市場的供需平衡 ,但是,實(shí)際的市場價格不會恰好等于均衡價格,而且價格也不會是靜態(tài)的 ,應(yīng)是隨時間不斷變化的動態(tài)過程.dvzfvkwMI1試建立描述市場價格形成的動態(tài)過程的數(shù)學(xué)模型假設(shè)在某一時刻,商品的價格為,它與該商品的均衡價格間有差別,此時,存在供需差,此供需差促使價格變動.對新的價格,又有新的供需差,如此不斷調(diào)節(jié),就構(gòu)成市場價格形成的動態(tài)過程,假設(shè)價格的變化率 與需求和供給之差成正比 ,并記 為需求函數(shù),為供給函數(shù)<為參數(shù)）,于是rqyn14ZNXI12/16個人資料整理 僅限學(xué)習(xí)使用其中 為商品在 時刻的價格, 為正常數(shù).若設(shè) , ,則上式變?yōu)棰倨渲?均為正常數(shù),其解為.下面對所得結(jié)果進(jìn)行討論 :<1）設(shè) 為靜態(tài)均衡價格,則其應(yīng)滿足,即 ,于是得 ,從而價格函數(shù) 可寫為,令 ,取極限得這說明,市場價格逐步趨于均衡價格 .又若初始價格 ,則動態(tài)價格就維持在均衡價格 上,整個動態(tài)過程就化為靜態(tài)過程； EmxvxOtOco<2）由于,13/16個人資料整理 僅限學(xué)習(xí)使用所以,當(dāng) 時, , 單調(diào)下降向 靠攏；當(dāng) 時, ,單調(diào)增加向 靠攏.這說明：初始價格高于均衡價格時 ,動態(tài)價格就要逐步降低,且逐步靠近均衡價格。否則 ,動態(tài)價格就要逐步升高 .因此,式①在一定程度上反映了價格影響需求與供給 ,而需求與供給反過來又影響價格的動態(tài)過程 ,并指出了動態(tài)價格逐步向均衡價格靠攏的變化趨勢.SixE2yXPq5總結(jié)：通過這次與眾不同的計算方法課程的學(xué)習(xí)，我和合作同學(xué)在此間都收獲了很多。從我們的角度將從以下幾點(diǎn)談?wù)搶@個學(xué)期計算方法學(xué)習(xí)做一個總結(jié)。6ewMyirQFL1、首先從學(xué)習(xí)收獲的角度。從我們的交流過程中可以得到幾點(diǎn)。首先，我們對運(yùn)用計算解決數(shù)學(xué)問題有了一定的了解，以前大家學(xué)習(xí)高數(shù)的時候整天研究做題，但是對高數(shù)中的積分微分等方法的原理沒有一個基礎(chǔ)的理解。但是這次，我們能深刻的理解到計算的基本原理，這樣對知識的掌握有一個新的高度。其次，除了計算原理理解之外，在計算方法的學(xué)習(xí)和編程實(shí)現(xiàn)中也漸漸的明白了一個理論的公式計算<如同常微分中復(fù)雜的公式，高難度的運(yùn)算）如何轉(zhuǎn)化成程序去實(shí)現(xiàn)。在這個過程中，明白了計算機(jī)運(yùn)行，或者說是設(shè)計一個算法的基礎(chǔ)是什么，簡單的說是借用計算機(jī)告訴可重復(fù)的運(yùn)算來彌補(bǔ)其中其離散的缺點(diǎn)。這樣雖然不能帶到醉精確的數(shù)值但是可以得到一個你需要的數(shù)值。在這個過程中我們從教材上，從課堂上14/16個人資料整理 僅限學(xué)習(xí)使用去體會是怎么把一個問題去離散化，如何去用迭代等方式去逼近一個問題等等。通過學(xué)習(xí)我們基本上認(rèn)識到了怎么將一個現(xiàn)實(shí)的問題去用編程思想去考慮。這不僅僅體現(xiàn)在我們算法的理解上面，更是我們對編程的一次新的理解。kavU42VRUs2、其次從授課方式角度。