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

微分方程數(shù)值解課程設(shè)計(jì)報(bào)告班級：______________姓名：_________學(xué)號：___________成績：2017年6月21日目錄一、摘要 算量較大，每前進(jìn)一步需要計(jì)算四次函數(shù)值。在用龍格庫塔方法時(shí),要注意的選擇要合適，太大，會(huì)使計(jì)算量加大，太小，較大，可能會(huì)使誤差增大。因此選擇合適的很重要。我們要在考慮精度的基礎(chǔ)上，選擇合適的。2.2算法步驟 2.2.1、四階龍格－庫塔（R-K）方法流程圖：輸入待求微分方程、求解的自變量范圍、初值以及求解范圍內(nèi)的取點(diǎn)數(shù)等。輸入待求微分方程、求解的自變量范圍、初值以及求解范圍內(nèi)的取點(diǎn)數(shù)等。確定求解范圍內(nèi)的步長k=取點(diǎn)數(shù)？否求解：求解并輸出：是結(jié)束算法2.2.2、Adams4階外插法流程圖：輸入待求微分方程、求解的自變量范圍、初值以及求解范圍內(nèi)的取點(diǎn)數(shù)等。輸入待求微分方程、求解的自變量范圍、初值以及求解范圍內(nèi)的取點(diǎn)數(shù)等。確定求解范圍內(nèi)的步長k=取點(diǎn)數(shù)？否求解：求解并輸出：是結(jié)束算法2.2.3、實(shí)例求解流程：輸入求解的自變量范圍輸入求解的自變量范圍求出待求簡單微分方程的真值解用MATLAB自帶函數(shù)ode45求解待求微分方程結(jié)束用自編函數(shù)Adams4階外插法方法求解待求微分方程開始2.3用matlab編寫源程序 Matlab程序源代碼：定義Rk4.m文件functiondy=Rk4(x,y)dy=zeros(3,1);dy(1)=10*(-y(1)+y(2));dy(2)=28*y(1)-y(2)-y(1)*y(3);dy(3)=y(1)*y(2)-8*y(3)/3;end定義Adams4.m文件function[x,y,z]=adams4(x1,y1,z1,x2,y2,z2,x3,y3,z3,h)%Adams外插法kfy=0;ksy=0;kty=0;kfz=0;ksz=0;ktz=0;kfx=10*(y3-x3);%eval_r(abx);kfy=x3*(28-z3)-y3;%eval_r(aby);kfz=x3*y3-8/3*z3;%eval_r(abz);ksx=10*(y2-x2);ksy=x2*(28-z2)-y2;;ksz=x2*y2-8/3*z2;ktx=10*(y1-x1);kty=x1*(28-z1)-y1;ktz=x1*y1-8/3*z1;x=x3+h/12*(23*kfx-16*ksx+5*ktx);y=y3+h/12*(23*kfy-16*ksy+5*kty);z=z3+h/12*(23*kfz-16*ksz+5*ktz);end定義exe11.m文件[t,y]=ode45(Rk4,[0,30],[12,2,9])suptitle('Runge-Kutta4階法')%總標(biāo)題subplot(2,2,1);plot(t,y(:,1));gridon;legend('x關(guān)于t的變化關(guān)系圖',1);xlabel('t','FontSize',14);ylabel('x','FontSize',14);subplot(2,2,2);plot(t,y(:,2));gridon;legend('y關(guān)于t的變化關(guān)系圖',1);xlabel('t','FontSize',14);ylabel('y','FontSize',14);subplot(2,2,3)plot(t,y(:,3));gridon;legend('z關(guān)于t的變化關(guān)系圖',1);xlabel('t','FontSize',14);ylabel('z','FontSize',14);subplot(2,2,4)plot3(y(:,1),y(:,2),y(:,3));gridon;legend('x,y,z的空間關(guān)系圖',1);xlabel('x','FontSize',14);ylabel('y','FontSize',14);zlabel('y','FontSize',14);view(40,60);%鎖定同樣的視圖，便于比較定義exe12.m文件a=1:1:3000;b=1:1:3000;c=1:1:3000;t=1:1:3000;x1=5;y1=5;z1=10;x2=7.83;y2=14.30;z2=12.34;x3=15.32;y3=22.87;z3=28.07;fori=1:1:3000[x,y,z]=adams4(x1,y1,z1,x2,y2,z2,x3,y3,z3,0.01);a(i)=x;b(i)=y;c(i)=z;x1=x2;y1=y2;z1=z2;x2=x3;y2=y3;z2=z3;x3=x;y3=y;z3=z;fprintf('x=%f''y=%f''z=%f\n',x,y,z);%顯示迭代值endsuptitle('Adams4階外插法')%總標(biāo)題subplot(2,2,1);plot(t,a);gridon;legend('x關(guān)于t的變化關(guān)系圖',1);xlabel('t','FontSize',14);ylabel('x','FontSize',14);subplot(2,2,2);plot(t,b);gridon;legend('y關(guān)于t的變化關(guān)系圖',1);xlabel('t','FontSize',14);ylabel('y','FontSize',14);subplot(2,2,3);plot(t,c);gridon;legend('z關(guān)于t的變化關(guān)系圖',1);xlabel('t','FontSize',14);ylabel('z','FontSize',14);subplot(2,2,4)plot3(a,b,c);gridon;legend('x,y,z的空間關(guān)系圖',1);xlabel('x','FontSize',14);ylabel('y','FontSize',14);zlabel('y','FontSize',14);view(40,60);%鎖定同樣的視圖，便于比較2.4常微分方程數(shù)值解法應(yīng)用舉例 題目： MIT的氣象學(xué)家洛侖茲（E.Lorenz）在1963年研究大氣對天氣的影響時(shí)，提出了Lorenz方程： 其中，，此方程初值問題的解存在且唯一。當(dāng)時(shí)，Lorenz方城有兩個(gè)不穩(wěn)定的不動(dòng)點(diǎn)，和，一個(gè)不穩(wěn)定的不動(dòng)點(diǎn)。當(dāng)時(shí)，和都變成不穩(wěn)定的。此時(shí)，存在混沌和一個(gè)奇怪吸引子。解：由于時(shí)，和都變成不穩(wěn)定的。此時(shí)，存在混沌和一個(gè)奇怪吸引子。所以我們?nèi)韺ｉT研究該現(xiàn)象。利用Matlab程序計(jì)算：命令窗口輸入： Runge-Kutta4階法：>>exe11Adams4階外插法：>>exe12結(jié)果輸出： 由于迭代次數(shù)的關(guān)系，結(jié)果數(shù)據(jù)量很大，故這里以圖片來展示結(jié)果。Runge-Kutta4階法結(jié)果：Adams4階外插法結(jié)果： 三、常系數(shù)擴(kuò)散方程的經(jīng)典差分格式3.1有限差分法的基本思路 用有限差分法解常系數(shù)擴(kuò)散方程有加權(quán)隱式差分格式其中，當(dāng)時(shí)為Crank-Nicolson格式，當(dāng)時(shí)為向后差分格式，當(dāng)時(shí)為向前差分格式。加權(quán)隱式格式穩(wěn)定的條件是加權(quán)隱式格式是兩層隱式格式，用第n層計(jì)算第n+1層節(jié)點(diǎn)值的時(shí)候，要解線性方程組。3.2算法步驟 實(shí)驗(yàn)步驟如下：(1)輸入，確定加權(quán)隱式格式的參數(shù)；(2)定義向量v，把初邊值條件離散，得到的值存入向量v；(3)利用差分格式由第n層計(jì)算第n+1層，建立相應(yīng)線性方程組，求解并且存入向量v；(4)計(jì)算到，輸出。3.3用matlab編寫源程序 Matlab程序源代碼：定義parabola.m文件function[u,tt,uu]=parabola(k,q,l)%q=thetal=taoh=0.1;%h:空間步長t=l*h*h;x=0.1:0.1:0.9;%x坐標(biāo)取值forn=1:9u(n)=sin(x(n)*pi);%初始值endu=u';v=zeros(9,1);%初始化v矩陣forn=1:kA=zeros(9,9);A(1,1)=1+2*q*t/h/h;fori=2:9A(i,i)=1+2*q*t/h/h;A(i-1,i)=-q*t/h/h;A(i,i-1)=-q*t/h/h;endB=zeros(9,1);B(1,1)=(1-q)*t/h/h*(u(2)-2*u(1))+u(1);B(9,1)=(1-q)*t/h/h*(u(8)-2*u(9))+u(9);fori=2:8B(i,1)=(1-q)*t/h/h*(u(i-1)-2*u(i)+u(i+1))+u(i);%加權(quán)隱式格式endv=inv(A)*B;%求解矩陣vw=u;u=v;v=w;n=n+1;endtt=k*t;uu=zeros(9,1);%初始化uu矩陣fori=1:9uu(i,1)=sin(pi*x(i))*exp(-pi*pi*tt);%精確值end定義exe20.m文件[u,tt,uu]=parabola(100,0.25,0.1);fprintf('τ=0.25，θ=0.1時(shí)0.4處的準(zhǔn)確值為：%f\n',uu(4))[u,tt,uu]=parabola(100,0.25,0.1);fprintf('τ=0.25，θ=0.1時(shí)0.4處的數(shù)值值為：%f\n',u(4))fprintf('誤差為：%f\n\n\n',abs(uu(4)-u(4)))x=0.1:0.1:0.9;subplot(2,2,1);plot(x,u,x,uu,'o');gridon;title('τ=0.25，θ=0.1');[u,tt,uu]=parabola(100,0.25,0.5);fprintf('τ=0.25，θ=0.5時(shí)0.4處的準(zhǔn)確值為：%f\n',uu(4))[u,tt,uu]=parabola(100,0.25,0.5);fprintf('τ=0.25，θ=0.5時(shí)0.4處的數(shù)值值為：%f\n',u(4))fprintf('誤差為：%f\n\n\n',abs(uu(4)-u(4)))x=0.1:0.1:0.9;subplot(2,2,2);plot(x,u,x,uu,'o');gridon;title('τ=0.25，θ=0.5');[u,tt,uu]=parabola(100,0.5,0.1);fprintf('τ=0.5，θ=0.1時(shí)0.4處的準(zhǔn)確值為：%f\n',uu(4))[u,tt,uu]=parabola(100,0.5,0.1);fprintf('τ=0.5，θ=0.1時(shí)0.4處的數(shù)值值為：%f\n',u(4))fprintf('誤差為：%f\n\n\n',abs(uu(4)-u(4)))x=0.1:0.1:0.9;subplot(2,2,3);plot(x,u,x,uu,'o');gridon;title('τ=0.5，θ=0.1');[u,tt,uu]=parabola(100,0.5,0.5);fprintf('τ=0.5，θ=0.5時(shí)0.4處的準(zhǔn)確值為：%f\n',uu(4))[u,tt,uu]=parabola(100,0.5,0.5);fprintf('τ=0.5，θ=0.5時(shí)0.4處的數(shù)值值為：%f\n',u(4))fprintf('誤差為：%f\n\n\n',abs(uu(4)-u(4)))x=0.1:0.1:0.9;subplot(2,2,4);plot(x,u,x,uu,'o');gridon;title('τ=0.5，θ=0.5');3.4有限差分法應(yīng)用舉例 題目：考慮常系數(shù)擴(kuò)散方程的初邊值問題其解析解為用加權(quán)隱式格式近似求解其中，取為時(shí)間步長，為網(wǎng)格比，對不同的時(shí)間步長()，計(jì)算當(dāng)時(shí)初邊值問題的解，并且與精確解比較，分析比較結(jié)果。用有限差分法近似求解，并且與精確解比較，分析結(jié)果。解：利用Matlab程序計(jì)算：命令窗口輸入：exe20結(jié)果輸出：τ=0.25，θ=0.1時(shí)0.4處的準(zhǔn)確值為：0.354466τ=0.25，θ=0.1時(shí)0.4處的數(shù)值值為：0.356486誤差為：0.002020τ=0.25，θ=0.5時(shí)0.4處的準(zhǔn)確值為：0.006840τ=0.25，θ=0.5時(shí)0.4處的數(shù)值值為：0.006696誤差為：0.000143τ=0.5，θ=0.1時(shí)0.4處的準(zhǔn)確值為：0.354466τ=0.5，θ=0.1時(shí)0.4處的數(shù)值值為：0.357343誤差為：0.002877τ=0.5，θ=0.5時(shí)0.4處的準(zhǔn)確值為：0.006840τ=0.5，θ=0.5時(shí)0.4處的數(shù)值值為：0.007115誤差為：0.000275誤差結(jié)論：根據(jù)圖和數(shù)據(jù)可以看出：當(dāng)，時(shí)，得到的結(jié)果誤差最小。當(dāng)改變q值和k值，可得到相應(yīng)的結(jié)果。四、橢圓型方程的五點(diǎn)差分格式4.1五點(diǎn)差分法的基本思路 對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)立大型方程組，該方程組可以用迭代法求解。4.2算法步驟 主要步驟：(1)首先取定，對求解區(qū)域劃分網(wǎng)格，按照網(wǎng)格定義矩陣v1，使得矩陣?yán)锩婷總€(gè)元素對應(yīng)求解區(qū)域中的每個(gè)節(jié)點(diǎn)；(2)由邊界條件定義v1矩陣中邊界元素的值，其余元素定義為零；(3)定義與v1同型的零矩陣v2；(4)用五點(diǎn)格式公式通過矩陣v1迭代計(jì)算矩陣v2，迭代精度為0.1；(5)畫圖。4.3用matlab編寫源程序 Matlab程序源代碼：定義fivepoint.m文件function[peuxyk]=fivepoint(h,m,n,kmax,ep)%h步長%m,n為x,y方向的網(wǎng)格數(shù)，例如（2-0）/0.01=200;%kmax為最大迭代次數(shù)%e為誤差，p為精確解symstemp;u=zeros(n+1,m+1);x=0+(0:m)*h;y=0+(0:n)*h;for(i=1:n+1)u(i,1)=sin(pi*y(i));%(0,y)u(i,m+1)=exp(1)*exp(1)*sin(pi*y(i));%(17,y)endfor(i=1:n)for(j=1:m)f(i,j)=(pi*pi-1)*exp(x(j))*sin(pi*y(i));%求區(qū)域中方程的每個(gè)節(jié)點(diǎn)endendt=zeros(n-1,m-1);for(k=1:kmax)for(i=2:n)for(j=2:m)%五點(diǎn)差分temp=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));if(k>kmax)break;%達(dá)到最大迭代次數(shù)，終止循環(huán)endif(max(max(t))<ep)break;endendfor(i=1:n+1)for(j=1:m+1)p(i,j)=exp(x(j))*sin(pi*y(i));%精確解e(i,j)=abs(u(i,j)-exp(x(j))*sin(pi*y(i)));%誤差endend4.4五點(diǎn)差分法應(yīng)用舉例 題目： 給定如下Laplace方程（poisson方程的特殊情況）的定解問題利用橢圓型方程的五點(diǎn)格式，計(jì)算該問題的近似解，并且畫出近似解的圖形。解：利用Matlab程序計(jì)算：命令窗口輸入：[peuxyk]=fivepoint(0.025,80,40,10000,1e-11);k=3952;figuresurf(x,y,u);xlswrite('D:\define.xls',e,'data1');%保存數(shù)據(jù)xlabel('x');ylabel('y');zlabel('u');title('五點(diǎn)差分法解橢圓型偏微分方程');figuresurf(x,y,e);xlswrite('D:\wucha.xls',e,'data2');%保存數(shù)據(jù)title('步長為0.025時(shí)的誤差');figuresurf(x,y,p);xlswrite('D:\jinshi.xls',e,'data32');%保存數(shù)據(jù)title('精確解曲面圖');結(jié)果輸出： 精確解近似解誤差2.02E-052.02E-052.02E-054.04E-054.04E-054.04E-056.02E-056.02E-056.02E-057.97E-057.97E-057.97E-059.87E-059.87E-059.87E-050.0001170.0001170.0001170.0001350.0001350.0001350.0001520.0001520.0001520.0001680.0001680.0001680.0001820.0001820.0001820.0001960.0001960.0001960.0002090.0002090.0002090.000220.000220.000220.000230.000230.000230.0002380.0002380.0002380.0002450.0002450.0