常微分方程課程設(shè)計(jì)PPT學(xué)習(xí)教案

1、會(huì)計(jì)學(xué)1常微分方程課程設(shè)計(jì)常微分方程課程設(shè)計(jì)q 自牛頓發(fā)明微積分以來，微分方程在描述事物運(yùn)自牛頓發(fā)明微積分以來，微分方程在描述事物運(yùn)動(dòng)規(guī)律上已發(fā)揮了重要的作用。實(shí)際應(yīng)用問題通過動(dòng)規(guī)律上已發(fā)揮了重要的作用。實(shí)際應(yīng)用問題通過數(shù)學(xué)建模所得到的方程，絕大多數(shù)是微分方程。數(shù)學(xué)建模所得到的方程，絕大多數(shù)是微分方程。q 由于實(shí)際應(yīng)用的需要，人們必須求解微分方程。然而能夠求得解析解的微分方程十分有限，絕大多數(shù)微分方程需要利用數(shù)值方法來近似求解。q 本實(shí)驗(yàn)主要研究如何用 Matlab 來計(jì)算微分方程（組）的數(shù)值解，并重點(diǎn)介紹一個(gè)求解微分方程的基本數(shù)值解法Euler折線法。第1頁/共25頁q 考慮一維經(jīng)典初值問題

2、00 , , ( ,)() , dyf x yy xyxa bdx u 基本思想：用差商代替微商根據(jù) Talyor 公式，y(x) 在點(diǎn) xk 處有2()()() ()()kkky xy xxxyxOx 1kkhxx 11()()()()( )kkkkky xy xy xy xdyO hdxhhx 21()()()()kkky xy xhyxO h 第2頁/共25頁q 具體步驟：等距剖分：0121nnaxxxxxb 步長：1 0 1 21() /, , ,kkhxxknnba u 分割求解區(qū)間u 差商代替微商1()()kkky xy xdydxhx 1 ()()()kkky xy xh yx

3、得方程組：0011 ()(,)kkkkkkyy xyyh f xyxxh 分割求解區(qū)間，差商代替微商，解代數(shù)方程 為分割點(diǎn) 0nkkx k = 0, 1, 2, ., n-1yk 是 y (xk) 的近似第3頁/共25頁例：用 Euler 法解初值問題22 0 201 , ( )dyxydxyxy 取步長 h = (2 - 0)/n = 2/n，得差分方程00110 1 2 ,(,)(/)kkkkkkkkkkxyyyh f xyyh yxyxxh 當(dāng) ，即 n=5 時(shí)，Matlab 源程序見 解：第4頁/共25頁clearf=sym(y+2*x/y2);a=0; b=2;h=0.4;n=(b-

4、a)/h+1; % n=(b-a)/h;x=0; y=1;szj=x,y;for i=1:n-1 % i=1:n y=y+h*subs(f,x,y,x,y);%第一次代初值 x=x+h; szj=szj;x,y;endszjplot(szj(:,1),szj(:,2),or-)第5頁/共25頁解析解：1 3352233/xyex 解析解近似解第6頁/共25頁q 為了減小誤差，可采用以下方法：u 讓步長 h 取得更小一些；u 改用具有較高精度的數(shù)值方法：q 龍格-庫塔方法Runge-Kutta (龍格-庫塔) 方法u 是一類求解常微分方程的數(shù)值方法u 有多種不同的迭代格式第7頁/共25頁q 用得

5、較多的是 四階R-K方法（教材第 79 頁）00111234 (22)/6(),kkkkyy xxxhyyh LLLL 12132432222(,)(/ ,/ )(/ ,/ )(,)kkkkkkkkLf xyLf xhyhLLf xhyhLLf xh yhL 其中第8頁/共25頁clear;f=sym(y+2*x/y2);a=0; b=2; h=0.4;n=(b-a)/h+1; % n=(b-a)/h;x=0; y=1; szj=x,y;for i=1:n-1 % i=1:n l1=subs(f,x,y,x,y); l2=subs(f,x,y,x+h/2,y+l1*h/2); l3=subs(

6、f,x,y,x+h/2,y+l2*h/2); l4=subs(f,x,y,x+h,y+l3*h); y=y+h*(l1+2*l2+2*l3+l4)/6; x=x+h; szj=szj;x,y;endplot(szj(:,1),szj(:,2), dg-)第9頁/共25頁第10頁/共25頁第11頁/共25頁q 用 Maltab自帶函數(shù) 解初值問題u 求解析解：dsolveu 求數(shù)值解： ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb第12頁/共25頁q dsolve 的使用y=dsolve(eq1,eq2, . ,cond1,cond2, .

7、,v)其中 y 為輸出， eq1、eq2、.為微分方程，cond1、cond2、.為初值條件，v 為自變量。例 1：求微分方程 的通解，并驗(yàn)證。22xdyxyxedx y=dsolve(Dy+2*x*y=x*exp(-x2),x) syms x; diff(y)+2*x*y - x*exp(-x2)y = 1/2*exp(-x2)*x2+exp(-x2)*C1第13頁/共25頁q 幾點(diǎn)說明l 如果省略初值條件，則表示求通解；l 如果省略自變量，則默認(rèn)自變量為 t dsolve(Dy=2*x,x); dy/dx = 2xdsolve(Dy=2*x); dy/dt = 2xl 若找不到解析解，則返

8、回其積分形式。l 微分方程中用 D 表示對(duì) 自變量 的導(dǎo)數(shù)，如：Dy y； D2y y； D3y y第14頁/共25頁例 2：求微分方程 在初值條件 下的特解，并畫出解函數(shù)的圖形。0 xxyye y=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x) ezplot(y);12( )ye 第15頁/共25頁例3：求微分方程組 在初值條件 下的特解，并畫出解函數(shù)的圖形。530tdxxyedtdyxydt x,y=dsolve(Dx+5*x+y=exp(t),Dy-x-3*y=0, . x(0)=1, y(0)=0, t)ezplot(x,y,0,1.3);0010|tt

9、xy 第16頁/共25頁例：x,y=dsolve(Dx+5*x=0,Dy-3*y=0, . x(0)=1, y(0)=1,t) r = dsolve(Dx+5*x=0,Dy-3*y=0, . x(0)=1, y(0)=1,t)這里返回的 r 是一個(gè) 結(jié)構(gòu)類型 的數(shù)據(jù)r.x %查看解函數(shù) x(t)r.y %查看解函數(shù) y(t)只有很少一部分微分方程（組）能求出解析解。大部分微分方程（組）只能利用數(shù)值方法求數(shù)值解。dsolve的輸出個(gè)數(shù)只能為一個(gè) 或 與方程個(gè)數(shù)相等。第17頁/共25頁T,Y = solver(odefun,tspan,y0)其中 y0 為初值條件，tspan為求解區(qū)間；Matla

10、b在數(shù)值求解時(shí)自動(dòng)對(duì)求解區(qū)間進(jìn)行分割，T (向量) 中返回的是分割點(diǎn)的值(自變量)，Y (向量) 中返回的是解函數(shù)在這些分割點(diǎn)上的函數(shù)值。solver 為Matlab的ODE求解器（可以是 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb）沒有一種算法可以有效地解決所有的 ODE 問題，因此MATLAB 提供了多種ODE求解器，對(duì)于不同的ODE，可以調(diào)用不同的求解器。第18頁/共25頁odefun 為顯式常微分方程，可以用命令 inline 定義，或在函數(shù)文件中定義，然后通過函數(shù)句柄調(diào)用。fun=inline(-2*y+2*x2+2*x,x,y);

11、x,y=ode23(fun,0,0.5,1);注：也可以在 tspan 中指定對(duì)求解區(qū)間的分割，如：x,y=ode23(fun,0:0.1:0.5,1); %此時(shí) x=0:0.1:0.5T,Y = solver(odefun,tspan,y0) 求初值問題 的數(shù)值解，求解范圍為 0,0.5222201( )dyyxxdxy 例 4：第19頁/共25頁如果需求解的問題是高階常微分方程，則需將其化為一階常微分方程組，此時(shí)需用函數(shù)文件來定義該常微分方程組。122212112101 00 7/()( ),( ),dxdtxdxdtxxxxx 令 ，則原方程可化為12,dyxy xdt 求解 Ver d

12、er Pol 初值問題2221001 00 7()( ),( ),d ydyyydtdtyy 例 5：第20頁/共25頁l 先編寫函數(shù)文件 function xprime=verderpol(t,x)global mu;xprime=x(2); mu*(1-x(1)2)*x(2) - x(1);l 再編寫腳本文件 ，在命令窗口直接運(yùn)行該文件。 clear;global mu;mu=7;y0=1;0;t,x=ode45(verderpol,0,40,y0); plot(t,x(:,1),r-);第21頁/共25頁q Matlab 函數(shù)u 求解析解（通解或特解），用 dsolveu 求數(shù)值解（特解），用