微分方程數(shù)值解實(shí)驗(yàn)報(bào)告

1、-. z.微分方程數(shù)值解法課程設(shè)計(jì)報(bào)告班級(jí)：_*： _*：_成績(jī)：2017年 6月 21 日摘要自然界與工程技術(shù)中的很多現(xiàn)象，可以歸結(jié)為微分方程定解問題。其中，常微分方程求解是微分方程的重要根底內(nèi)容。但是，對(duì)于許多的微分方程，往往很難得到甚至不存在準(zhǔn)確的解析表達(dá)式，這時(shí)候，數(shù)值解提供了一個(gè)很好的解決思路。，針對(duì)于此，本文對(duì)常微分方程數(shù)值解法進(jìn)展了簡(jiǎn)單研究，主要討論了一些常用的數(shù)值解法，如歐拉法、改良的歐拉法、RungeKutta方法、Adams法以及橢圓型方程、拋物型方程的有限差分方法等，通過具體的算例，結(jié)合MATLAB求解畫圖，初步給出了一般常微分方程數(shù)值解法的求解過程。同時(shí)，通過對(duì)各種方法

2、的誤差分析，讓大家對(duì)各種方法的特點(diǎn)和適用*圍有一個(gè)直觀的感受。關(guān)鍵詞：微分方程數(shù)值解、MATLAB 目錄HYPERLINK l _Toc313386949摘要2 TOC o 1-3 h z u HYPERLINK l _Toc313386949目錄3HYPERLINK l _Toc313386950第一章常微分方程數(shù)值解法的根本思想與原理4HYPERLINK l _Toc3133869541.1常微分方程數(shù)值解法的根本思路4HYPERLINK l _Toc3133869561.2用matlab編寫源程序41HYPERLINK l _Toc313386957.3常微分方程數(shù)值解法應(yīng)用舉例及結(jié)果5

3、HYPERLINK l _Toc313386950第二章常系數(shù)擴(kuò)散方程的經(jīng)典差分格式的根本思想與原理6HYPERLINK l _Toc3133869542.1常系數(shù)擴(kuò)散方程的經(jīng)典差分格式的根本思路6HYPERLINK l _Toc3133869562.2用matlab編寫源程序72HYPERLINK l _Toc313386957.3常系數(shù)擴(kuò)散方程的經(jīng)典差分格式的應(yīng)用舉例及結(jié)果8HYPERLINK l _Toc313386952第三章橢圓型方程的五點(diǎn)差分格式的根本思想與原理10HYPERLINK l _Toc3133869543.1橢圓型方程的五點(diǎn)差分格式的根本思路103HYPERLINK l

4、 _Toc313386956.2用matlab編寫源程序103HYPERLINK l _Toc313386957.3橢圓型方程的五點(diǎn)差分格式的應(yīng)用舉例及結(jié)果12HYPERLINK l _Toc313386961第四章總結(jié)12HYPERLINK l _Toc313386961參考文獻(xiàn)12第一章常微分方程數(shù)值解法的根本思想與原理1.1常微分方程數(shù)值解法的根本思路常微分方程數(shù)值解法(numerical methods forordinary differential equations)計(jì)算數(shù)學(xué)的一個(gè)分支.是解常微分方程各類定解問題的數(shù)值方法.現(xiàn)有的解析方法只能用于求解一些特殊類型的定解問題，實(shí)用上

5、許多很有價(jià)值的常微分方程的解不能用初等函數(shù)來表示，常常需要求其數(shù)值解.所謂數(shù)值解，是指在求解區(qū)間內(nèi)一系列離散點(diǎn)處給出真解的近似值.這就促成了數(shù)值方法的產(chǎn)生與開展.1.2用matlab編寫源程序龍格庫塔法：M文件:function d*=Lorenz(t,*) %r=28,sigma=10,b=8/3d*=-10*(*(1)-*(2);-*(1)*(3)+28*(1)-*(2);*(1)*(2)-8*(3)/3;運(yùn)行程序：*0=1,1,1;t,y=ode45(Lorenz,0,100,*0);subplot(2,1,1) %兩行一列的圖第一個(gè)plot(t,y(:,3)*label(time);y

6、label(z);%畫z-t圖像subplot(2,2,3)%兩行兩列的圖第三個(gè)plot(y(:,1),y(:,2)*label(*);ylabel(y); %畫*-y圖像subplot(2,2,4)plot3(y(:,1),y(:,2),y(:,3)*label(*);ylabel(y);zlabel(z);%畫*yz圖像歐拉法：h=0.010;a=16;b=4;c=49.52;*=5;y=10;z=10;Y=;for i=1:800 *1=*+h*a*(y-*); y1=y+h*(c*-*z-y); z1=z+h*(*y-b*z); *=*1; y=y1; z=z1; Y(i,:)=* y

7、 z;endplot3(Y(:,1),Y(:,2),Y(:,3);1.3常微分方程數(shù)值解法的應(yīng)用舉例及結(jié)果應(yīng)用舉例：a=10,b=8/3,0r+,當(dāng)1r24.74時(shí)，c和都變成不穩(wěn)定的，此時(shí)存在混沌和奇怪吸引子。運(yùn)行結(jié)果：龍格庫塔法：歐拉法：第二章常系數(shù)擴(kuò)散方程的經(jīng)典差分格式的根本思想與原理2.1 常系數(shù)擴(kuò)散方程的經(jīng)典差分格式的根本思路用有限差分法解常系數(shù)擴(kuò)散方程有加權(quán)隱式差分格式其中，當(dāng)時(shí)為Crank-Nicolson格式，當(dāng)時(shí)為向后差分格式，當(dāng)時(shí)為向前差分格式。加權(quán)隱式格式穩(wěn)定的條件是，當(dāng)，無限制，當(dāng)。加權(quán)隱式格式是兩層隱式格式，用第n層計(jì)算第n+1層節(jié)點(diǎn)值的時(shí)候，要解線性方程組。2.2

8、用matlab編寫源程序M文件：function M = chase(a,b,c,f)% 追趕法求解三對(duì)角矩陣方程，A*=f% a是對(duì)角線下邊一行的元素% b是對(duì)角線元素% c是對(duì)角線上邊一行的元素% M是求得的結(jié)果，以列向量形式保存n = length(b);beta = ones(1,n-1); y = ones(1,n); M = ones(n,1);for i = (n-1):(-1):1 a(i+1) = a(i);end% 將a矩陣和n對(duì)應(yīng)beta(1) = c(1)/b(1);for i = 2:(n-1) beta(i) = c(i)/( b(i)-a(i)*beta(i-1)

9、 );endy(1) = f(1)/b(1);for i = 2:n y(i) = (f(i)-a(i)*y(i-1)/(b(i)-a(i)*beta(i-1);endM(n) = y(n);for i = (n-1):(-1):1 M(i) = y(i)-beta(i)*M(i+1);end endM文件：function output = diffuse_equation(a0,t_ma*,h,tao,D,a1,b1,c1,a2,b2,c2)% 一維擴(kuò)散方程的有限差分法，采用隱式六點(diǎn)差分格式(Crank-Nicolson)% a0: *的最大值% t:_ma*: t的最大值% h: 空間步

10、長% tao: 時(shí)間步長% D：擴(kuò)散系數(shù)% a1,b1,c1是*=0邊界條件的系數(shù)；a2,b2,c2是*=a0邊界條件的系數(shù)* = 0:h:a0;n = length(*);t = 0:tao:t_ma*;k = length(t); P = tao * D/h2;P1 = 1/P + 1;P2 = 1/P - 1;u = zeros(k,n);%初始條件u(1,:) = e*p(*);%求A矩陣的對(duì)角元素dd = zeros(1,n);d(1,1) = b1*P1+h*a1;d(2:(n-1),1) = 2*P1;d(n,1) = b2*P1+h*a2;%求A矩陣的對(duì)角元素下面一行元素ee

11、= -ones(1,n-1);e(1,n-1) = -b2;%求A矩陣的對(duì)角元素上面一行元素ff = -ones(1,n-1);f(1,1) = -b1;R = zeros(k,n);%求R%追趕法求解for i = 2:k R(i,1) = (b1*P2-h*a1)*u(i-1,1)+b1*u(i-1,2)+2*h*c1; for j = 2:n-1 R(i,j) = u(i-1,j-1)+2*P2*u(i-1,j)+u(i-1,j+1); end R(i,n) = b2*u(i-1,n-1)+( b2*P2-h*a2)*u(i-1,n)+2*h*c2; M = chase(e,d,f,R(

12、i,:); u(i,:) = M; plot(*,u(i,:); a*is(0 a0 0 t_ma*); pause(0.1)endoutput = u% 繪圖比擬解析解和有限差分解*,T = meshgrid(*,t);Z = e*p(-pi.*pi.*T).*sin(pi.*);surf(*,T,Z),*label(*),ylabel(t),zlabel(u),title(解析解);%colormap(gray(1);%使圖向變?yōu)楹谏玣iguresurf(*,T,u),*label(*),ylabel(t),zlabel(u),title(有限差分解);%colormap(gray(1);

13、%使圖向變?yōu)楹谏\(yùn)行程序：% 一維擴(kuò)散方程的有限差分法clear,clc;%定義初始常量a1 = 1; b1 = 1; c1 = 0; a2 = 1;b2 = -1; c2 = 0;a0 = 1.0; t_ma* = 8; D = 0.1; h = 0.1; tao = 0.1;%調(diào)用擴(kuò)散方程子函數(shù)求解u = diffuse_equation(a0,t_ma*,h,tao,D,a1,b1,c1,a2,b2,c2);2.3常系數(shù)擴(kuò)散方程的經(jīng)典差分格式的應(yīng)用舉例及結(jié)果應(yīng)用舉例：考慮常系數(shù)擴(kuò)散方程的初邊值問題其中，取,為時(shí)間步長，為網(wǎng)格比，對(duì)不同的時(shí)間步長，計(jì)算當(dāng)時(shí)初邊值問題的解u(0.4,0.4)

14、,并且與準(zhǔn)確解比擬，分析比擬結(jié)果。運(yùn)行結(jié)果：第三章橢圓型方程的五點(diǎn)差分格式的根本思想與原理3.1橢圓型方程的五點(diǎn)差分格式的根本思路對(duì)Laplace方程的第一邊值問題利用taylor展開可得逼近它的五點(diǎn)差分格式的差分逼近其中分別為軸和軸步長，邊界條件可以由離散可得，當(dāng)時(shí)有。注意五點(diǎn)格式計(jì)算節(jié)點(diǎn)是由邊界的節(jié)點(diǎn)，計(jì)算內(nèi)部節(jié)點(diǎn)，計(jì)算時(shí)需要聯(lián)立大型方程組，該方程組可以用迭代法求解。3.2用matlab編寫源程序M文件：function p e u * y k=wudianchafenfa(h,m,n,kma*,ep) % g-s迭代法解五點(diǎn)差分法問題 %kma*為最大迭代次數(shù) %m,n為*,y方向的網(wǎng)格

15、數(shù)，例如2-0/0.01=200; %e為誤差，p為準(zhǔn)確解 syms temp; u=zeros(n+1,m+1); *=0+(0:m)*h; y=0+(0:n)*h; for i=1:n+1 u(i,1)=sin(pi*y(i); u(i,m+1)=e*p(1)*e*p(1)*sin(pi*y(i); endfor i=1:n for j=1:m f (i,j)=(pi*pi-1)*e*p(*(j)*sin(pi*y(i); endendt=zeros(n-1,m-1);for k=1:kma* for i=2:n for j=2:m temp=h*h*f(i,j)/4+(u(i,j+1)+

16、u(i,j-1)+u(i+1,j)+u(i-1,j)/4; t(i,j)=(temp-u(i,j)*(temp-u(i,j); u(i,j)=temp; end end t(i,j)=sqrt(t(i,j); if kkma* break; end if ma*(ma*(t)ep break; end end for i=1:n+1 for j=1:m+1 p(i,j)=e*p(*(j)*sin(pi*y(i); e(i,j)=abs(u(i,j)-e*p(*(j)*sin(pi*y(i); end end運(yùn)行程序：p e u * y k=wudianchafenfa(0.1,20,10,10000,1e-6); surf(*,y,u); *label(*);ylabel(y);zlabel(u); title(五點(diǎn)差分法解橢圓型偏微分方程);3.3橢圓型方程的五點(diǎn)差分格式的應(yīng)用舉例及結(jié)果應(yīng)用舉