微分方程數(shù)值解實驗報告

微分方程數(shù)值解法課程設(shè)計報告班級：姓名：學(xué)號：成績：2017年6月21日摘要自然界與工程技術(shù)中的很多現(xiàn)象，可以歸結(jié)為微分方程定解問題。其中，常微分方程求解是微分方程的重要基礎(chǔ)內(nèi)容。但是，對于許多的微分方程，往往很難得到甚至不存在精確的解析表達(dá)式，這時候，數(shù)值解提供了一個很好的解決思路。，針對于此，本文對常微分方程數(shù)值解法進(jìn)行了簡單研究，主要討論了一些常用的數(shù)值解法，如歐拉法、改進(jìn)的歐拉法、Runge—Kutta方法、Adams法以及橢圓型方程、拋物型方程的有限差分方法等，通過具體的算例，結(jié)合MATLAB求解畫圖，初步給出了一般常微分方程數(shù)值解法的求解過程。同時，通過對各種方法的誤差分析，讓大家對各種方法的特點和適用范圍有一個直觀的感受。關(guān)鍵詞：微分方程數(shù)值解、MATLAB目錄摘要 2目錄 3第一章常微分方程數(shù)值解法的基本思想與原理 41.1常微分方程數(shù)值解法的基本思路 41.2用matlab編寫源程序 41.3常微分方程數(shù)值解法應(yīng)用舉例及結(jié)果 5第二章常系數(shù)擴(kuò)散方程的經(jīng)典差分格式的基本思想與原理 62.1常系數(shù)擴(kuò)散方程的經(jīng)典差分格式的基本思路 62.2用matlab編寫源程序 72.3常系數(shù)擴(kuò)散方程的經(jīng)典差分格式的應(yīng)用舉例及結(jié)果 8第三章橢圓型方程的五點差分格式的基本思想與原理 103.1橢圓型方程的五點差分格式的基本思路 103.2用matlab編寫源程序 103.3橢圓型方程的五點差分格式的應(yīng)用舉例及結(jié)果 12第四章總結(jié) 12參考文獻(xiàn) 12第一章常微分方程數(shù)值解法的基本思想與原理1.1常微分方程數(shù)值解法的基本思路常微分方程數(shù)值解法(numericalmethodsforordinarydifferentialequations)計算數(shù)學(xué)的一個分支.是解常微分方程各類定解問題的數(shù)值方法.現(xiàn)有的解析方法只能用于求解一些特殊類型的定解問題，實用上許多很有價值的常微分方程的解不能用初等函數(shù)來表示，常常需要求其數(shù)值解.所謂數(shù)值解，是指在求解區(qū)間內(nèi)一系列離散點處給出真解的近似值.這就促成了數(shù)值方法的產(chǎn)生與發(fā)展.1.2用matlab編寫源程序龍格庫塔法：M文件:functiondx=Lorenz(t,x)%r=28,sigma=10,b=8/3

dx=[-10*(x(1)-x(2));-x(1)*x(3)+28*x(1)-x(2);x(1)*x(2)-8*x(3)/3];運行程序：x0=[1,1,1];

[t,y]=ode45('Lorenz',[0,100],x0);

subplot(2,1,1)%兩行一列的圖第一個

plot(t,y(:,3))xlabel('time');ylabel('z');%畫z-t圖像subplot(2,2,3)

%兩行兩列的圖第三個

plot(y(:,1),y(:,2))xlabel('x');ylabel('y');%畫x-y圖像subplot(2,2,4)

plot3(y(:,1),y(:,2),y(:,3))

xlabel('x');ylabel('y');zlabel('z');%畫xyz圖像歐拉法：h=0.010;a=16;b=4;c=49.52;x=5;y=10;z=10;Y=[];fori=1:800x1=x+h*a*(y-x);y1=y+h*(c*x-x*z-y);z1=z+h*(x*y-b*z);x=x1;y=y1;z=z1;Y(i,:)=[xyz];endplot3(Y(:,1),Y(:,2),Y(:,3));1.3常微分方程數(shù)值解法的應(yīng)用舉例及結(jié)果應(yīng)用舉例：a=10,b=8/3,0<r<+∞,當(dāng)1<r<24.74時，Lorenz方程有兩個穩(wěn)定的不動點c(，，r-1)和（-，-，r-1），一個穩(wěn)定的不動點0=（0，0，0），當(dāng)r>24.74時，c和都變成不穩(wěn)定的，此時存在混沌和奇怪吸引子。運行結(jié)果：龍格庫塔法：歐拉法：第二章常系數(shù)擴(kuò)散方程的經(jīng)典差分格式的基本思想與原理2.1常系數(shù)擴(kuò)散方程的經(jīng)典差分格式的基本思路用有限差分法解常系數(shù)擴(kuò)散方程有加權(quán)隱式差分格式其中，當(dāng)時為Crank-Nicolson格式，當(dāng)時為向后差分格式，當(dāng)時為向前差分格式。加權(quán)隱式格式穩(wěn)定的條件是，當(dāng)，無限制，當(dāng)。加權(quán)隱式格式是兩層隱式格式，用第n層計算第n+1層節(jié)點值的時候，要解線性方程組。2.2用matlab編寫源程序M文件：functionM=chase(a,b,c,f)%追趕法求解三對角矩陣方程，Ax=f%a是對角線下邊一行的元素%b是對角線元素%c是對角線上邊一行的元素%M是求得的結(jié)果，以列向量形式保存n=length(b);beta=ones(1,n-1);y=ones(1,n);M=ones(n,1);fori=(n-1):(-1):1a(i+1)=a(i);end%將a矩陣和n對應(yīng)beta(1)=c(1)/b(1);fori=2:(n-1)beta(i)=c(i)/(b(i)-a(i)*beta(i-1));endy(1)=f(1)/b(1);fori=2:ny(i)=(f(i)-a(i)*y(i-1))/(b(i)-a(i)*beta(i-1));endM(n)=y(n);fori=(n-1):(-1):1M(i)=y(i)-beta(i)*M(i+1);endendM文件：functionoutput=diffuse_equation(a0,t_max,h,tao,D,a1,b1,c1,a2,b2,c2)%一維擴(kuò)散方程的有限差分法，采用隱式六點差分格式(Crank-Nicolson)%a0:x的最大值%t:_max:t的最大值%h:空間步長%tao:時間步長%D：擴(kuò)散系數(shù)%a1,b1,c1是（x=0）邊界條件的系數(shù)；a2,b2,c2是（x=a0）邊界條件的系數(shù)x=0:h:a0;n=length(x);t=0:tao:t_max;k=length(t);P=tao*D/h^2;P1=1/P+1;P2=1/P-1;u=zeros(k,n);%初始條件u(1,:)=exp(x);%求A矩陣的對角元素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矩陣的對角元素下面一行元素ee=-ones(1,n-1);e(1,n-1)=-b2;%求A矩陣的對角元素上面一行元素ff=-ones(1,n-1);f(1,1)=-b1;R=zeros(k,n);%求R%追趕法求解fori=2:kR(i,1)=(b1*P2-h*a1)*u(i-1,1)+b1*u(i-1,2)+2*h*c1;forj=2:n-1R(i,j)=u(i-1,j-1)+2*P2*u(i-1,j)+u(i-1,j+1);endR(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(i,:));u(i,:)=M';plot(x,u(i,:));axis([0a00t_max]);pause(0.1)endoutput=u%繪圖比較解析解和有限差分解[X,T]=meshgrid(x,t);Z=exp(-pi.*pi.*T).*sin(pi.*X);surf(X,T,Z),xlabel('x'),ylabel('t'),zlabel('u'),title('解析解');%colormap(gray(1));%使圖向變?yōu)楹谏玣iguresurf(X,T,u),xlabel('x'),ylabel('t'),zlabel('u'),title('有限差分解');%colormap(gray(1));%使圖向變?yōu)楹谏\行程序：%一維擴(kuò)散方程的有限差分法clear,clc;%定義初始常量a1=1;b1=1;c1=0;a2=1;b2=-1;c2=0;a0=1.0;t_max=8;D=0.1;h=0.1;tao=0.1;%調(diào)用擴(kuò)散方程子函數(shù)求解u=diffuse_equation(a0,t_max,h,tao,D,a1,b1,c1,a2,b2,c2);2.3常系數(shù)擴(kuò)散方程的經(jīng)典差分格式的應(yīng)用舉例及結(jié)果應(yīng)用舉例：考慮常系數(shù)擴(kuò)散方程的初邊值問題其中，取,為時間步長，為網(wǎng)格比，對不同的時間步長，計算當(dāng)時初邊值問題的解u(0.4,0.4),并且與精確解比較，分析比較結(jié)果。運行結(jié)果：第三章橢圓型方程的五點差分格式的基本思想與原理3.1橢圓型方程的五點差分格式的基本思路對Laplace方程的第一邊值問題利用taylor展開可得逼近它的五點差分格式的差分逼近其中分別為軸和軸步長，邊界條件可以由離散可得，當(dāng)時有。注意五點格式計算節(jié)點是由邊界的已知節(jié)點，計算內(nèi)部節(jié)點，計算時需要聯(lián)立大型方程組，該方程組可以用迭代法求解。3.2用matlab編寫源程序M文件：function[peuxyk]=wudianchafenfa(h,m,n,kmax,ep)%g-s迭代法解五點差分法問題%kmax為最大迭代次數(shù)%m,n為x,y方向的網(wǎng)格數(shù)，例如（2-0）/0.01=200;%e為誤差，p為精確解symstemp;u=zeros(n+1,m+1);x=0+(0:m)*h;y=0+(0:n)*h;fori=1:n+1u(i,1)=sin(pi*y(i));u(i,m+1)=exp(1)*exp(1)*sin(pi*y(i));endfori=1:nforj=1:mf(i,j)=(pi*pi-1)*exp(x(j))*sin(pi*y(i));endendt=zeros(n-1,m-1);fork=1:kmaxfori=2:nforj=2:mtemp=h*h*f(i,j)/4+(u(i,j+1)+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;endendt(i,j)=sqrt(t(i,j));ifk>kmaxbreak;endifmax(max(t))<epbreak;endendfori=1:n+1forj=1:m+1p(i,j)=exp(x(j))*sin(pi*y(i));e(i,j)=abs(u(i,j)-exp(x(j))*sin(pi*y(i)));endend運行程序：[peuxyk]=wudianchafenfa(0.1,20,10,10000,1e-6);surf(x,y,u);xlabel('x');ylabel('y');zlabel('u');title('五點差分法解橢圓型偏微分方程');3.3橢圓型方程的五點差分格式的應(yīng)用舉例及結(jié)果應(yīng)用舉例：給定如下Lapla