常微分方程數(shù)值解法

1、常微分方程數(shù)值解法【作用】微分方程建模是數(shù)學建模的重要方法，因為許多實際問題的數(shù)學描述將導致求解微分方程的定解問題。把形形色色的實際問題化成微分方程的定解問題，大體上可以按以下幾步：1. 根據(jù)實際要求確定要研究的量(自變量、未知函數(shù)、必要的參數(shù)等)并確定坐標系。2. 找出這些量所滿足的基本規(guī)律(物理的、幾何的、化學的或生物學的等等)。3. 運用這些規(guī)律列出方程和定解條件。 基本模型1.發(fā)射衛(wèi)星為什么用三級火箭2.人口模型3.戰(zhàn)爭模型4.放射性廢料的處理通常需要求出方程的解來說明實際現(xiàn)象，并加以檢驗。如果能得到解析形式的解固然是便于分析和應用的，但是我們知道，只有線性常系數(shù)微分方程，并且自由項是

2、某些特殊類型的函數(shù)時，才可以得到這樣的解，而絕大多數(shù)變系數(shù)方程、非線性方程都是所謂“解不出來”的于是對于用微分方程解決實際問題來說，數(shù)值解法就是一個十分重要的手段。1. 改進Euler 法：2. 龍格庫塔（RungeKutta）方法：【源程序】1. 改進Euler 法：function x,y=eulerpro(fun,x0,x1,y0,n);%fun為函數(shù)，(x0，x1)為x區(qū)間，y0為初始值，n為子區(qū)間個數(shù)if nargin5,n=50;endh=(x1-x0)/n;x(1)=x0;y(1)=y0;for i=1:nx(i+1)=x(i)+h;y1=y(i)+h*feval(fun,x(i

3、),y(i);y2=y(i)+h*feval(fun,x(i+1),y1);y(i+1)=(y1+y2)/2;end調(diào)用 command窗口f=inline(-2*y+2*x2+2*x)x,y=eulerpro(f,0,0.5,1,10)求解函數(shù)y=2y+2+2x ，(0 x 0.5)， y(0) = 12. 龍格庫塔（RungeKutta）方法：t,y=solver(F,tspan，y0)這里solver為ode45，ode23，ode113，輸入?yún)?shù) F 是用M文件定義的微分方程y= f (x, y)右端的函數(shù)。tspan=t0,tfinal是求解區(qū)間，y0是初值。注：ode45和ode2

4、3變步長的，采用Runge-Kutta算法。ode45表示采用四階-五階Runge-Kutta算法，它用4階方法提供候選解，5階方法控制誤差，是一種自適應步長（變步長）的常微分方程數(shù)值解法，其整體截斷誤差為(x)5。解決的是Nonstiff(非剛性)常微分方程。ode45是解決數(shù)值解問題的首選方法，若長時間沒結(jié)果，應該就是剛性的，可換用ode23試試例如：odefun=(t,y) (y+3*t)/t2; %定義函數(shù)tspan=1 4; %求解區(qū)間y0=-2; %初值t,y=ode45(odefun,tspan,y0);plot(t,y) %作圖title(t2y=y+3t,y(1)=-2,1t

5、4)legend(t2y=y+3t)xlabel(t)ylabel(y)% 精確解dsolve(t2*Dy=y+3*t,y(1)=-2) ans = (3*Ei(1,-1/t)+(-3*exp(-1)*Ei(1,-1)-2)/exp(-1)*exp(-1/t)求解高階常微分方程需要求解的高階常微分方程：求解的關(guān)鍵是將高階轉(zhuǎn)為一階,odefun的書寫.F(y,y,y.y(n-1),t)=0用變量替換,y1=y,y2=y.注意odefun方程定義為行向量程序:function Testode45tspan=3.9 4.0; %求解區(qū)間y0=2 8; %初值t,x=ode45(odefun,tspan,y0);plot(t,x(:,1),-o,t,x(:,2),-*)legend(y1,y2)title(y =-t*y + et*y +3sin2t)xlabel(t)ylabel(y)endfunction y=odefun(t,x)y=zeros(2,1); % 列向量y(1)=x(2);y(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t);end【原理】改進Euler 法：式稱為由Euler 公式和梯形公式得到的預測校正系統(tǒng)，也叫改進Euler 法。為了編程方便，將式子改寫