微分方程数值解实验

微分方程数值解课程设计报告班级：______________姓名：_________学号：___________成绩：2017年6月21日目录一、摘要 算量较大，每前进一步需要计算四次函数值。在用龙格库塔方法时,要注意的选择要合适，太大，会使计算量加大，太小，较大，可能会使误差增大。因此选择合适的很重要。我们要在考虑精度的基础上，选择合适的。2.2算法步骤 2.2.1、四阶龙格－库塔（R-K）方法流程图：输入待求微分方程、求解的自变量范围、初值以及求解范围内的取点数等。输入待求微分方程、求解的自变量范围、初值以及求解范围内的取点数等。确定求解范围内的步长k=取点数？否求解：求解并输出：是结束算法2.2.2、Adams4阶外插法流程图：输入待求微分方程、求解的自变量范围、初值以及求解范围内的取点数等。输入待求微分方程、求解的自变量范围、初值以及求解范围内的取点数等。确定求解范围内的步长k=取点数？否求解：求解并输出：是结束算法2.2.3、实例求解流程：输入求解的自变量范围输入求解的自变量范围求出待求简单微分方程的真值解用MATLAB自带函数ode45求解待求微分方程结束用自编函数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阶法')%总标题subplot(2,2,1);plot(t,y(:,1));gridon;legend('x关于t的变化关系图',1);xlabel('t','FontSize',14);ylabel('x','FontSize',14);subplot(2,2,2);plot(t,y(:,2));gridon;legend('y关于t的变化关系图',1);xlabel('t','FontSize',14);ylabel('y','FontSize',14);subplot(2,2,3)plot(t,y(:,3));gridon;legend('z关于t的变化关系图',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的空间关系图',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阶外插法')%总标题subplot(2,2,1);plot(t,a);gridon;legend('x关于t的变化关系图',1);xlabel('t','FontSize',14);ylabel('x','FontSize',14);subplot(2,2,2);plot(t,b);gridon;legend('y关于t的变化关系图',1);xlabel('t','FontSize',14);ylabel('y','FontSize',14);subplot(2,2,3);plot(t,c);gridon;legend('z关于t的变化关系图',1);xlabel('t','FontSize',14);ylabel('z','FontSize',14);subplot(2,2,4)plot3(a,b,c);gridon;legend('x,y,z的空间关系图',1);xlabel('x','FontSize',14);ylabel('y','FontSize',14);zlabel('y','FontSize',14);view(40,60);%锁定同样的视图，便于比较2.4常微分方程数值解法应用举例 题目： MIT的气象学家洛仑兹（E.Lorenz）在1963年研究大气对天气的影响时，提出了Lorenz方程： 其中，，此方程初值问题的解存在且唯一。当时，Lorenz方城有两个不稳定的不动点，和，一个不稳定的不动点。当时，和都变成不稳定的。此时，存在混沌和一个奇怪吸引子。解：由于时，和都变成不稳定的。此时，存在混沌和一个奇怪吸引子。所以我们取来专门研究该现象。利用Matlab程序计算：命令窗口输入： Runge-Kutta4阶法：>>exe11Adams4阶外插法：>>exe12结果输出： 由于迭代次数的关系，结果数据量很大，故这里以图片来展示结果。Runge-Kutta4阶法结果：Adams4阶外插法结果： 三、常系数扩散方程的经典差分格式3.1有限差分法的基本思路 用有限差分法解常系数扩散方程有加权隐式差分格式其中，当时为Crank-Nicolson格式，当时为向后差分格式，当时为向前差分格式。加权隐式格式稳定的条件是加权隐式格式是两层隐式格式，用第n层计算第n+1层节点值的时候，要解线性方程组。3.2算法步骤 实验步骤如下：(1)输入，确定加权隐式格式的参数；(2)定义向量v，把初边值条件离散，得到的值存入向量v；(3)利用差分格式由第n层计算第n+1层，建立相应线性方程组，求解并且存入向量v；(4)计算到，输出。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坐标取值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);%加权隐式格式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时0.4处的准确值为：%f\n',uu(4))[u,tt,uu]=parabola(100,0.25,0.1);fprintf('τ=0.25，θ=0.1时0.4处的数值值为：%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时0.4处的准确值为：%f\n',uu(4))[u,tt,uu]=parabola(100,0.25,0.5);fprintf('τ=0.25，θ=0.5时0.4处的数值值为：%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时0.4处的准确值为：%f\n',uu(4))[u,tt,uu]=parabola(100,0.5,0.1);fprintf('τ=0.5，θ=0.1时0.4处的数值值为：%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时0.4处的准确值为：%f\n',uu(4))[u,tt,uu]=parabola(100,0.5,0.5);fprintf('τ=0.5，θ=0.5时0.4处的数值值为：%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有限差分法应用举例 题目：考虑常系数扩散方程的初边值问题其解析解为用加权隐式格式近似求解其中，取为时间步长，为网格比，对不同的时间步长()，计算当时初边值问题的解，并且与精确解比较，分析比较结果。用有限差分法近似求解，并且与精确解比较，分析结果。解：利用Matlab程序计算：命令窗口输入：exe20结果输出：τ=0.25，θ=0.1时0.4处的准确值为：0.354466τ=0.25，θ=0.1时0.4处的数值值为：0.356486误差为：0.002020τ=0.25，θ=0.5时0.4处的准确值为：0.006840τ=0.25，θ=0.5时0.4处的数值值为：0.006696误差为：0.000143τ=0.5，θ=0.1时0.4处的准确值为：0.354466τ=0.5，θ=0.1时0.4处的数值值为：0.357343误差为：0.002877τ=0.5，θ=0.5时0.4处的准确值为：0.006840τ=0.5，θ=0.5时0.4处的数值值为：0.007115误差为：0.000275误差结论：根据图和数据可以看出：当，时，得到的结果误差最小。当改变q值和k值，可得到相应的结果。四、椭圆型方程的五点差分格式4.1五点差分法的基本思路 对Laplace方程的第一边值问题利用taylor展开可得逼近它的五点差分格式的差分逼近其中分别为轴和轴步长，边界条件可以由离散可得，当时有。注意五点格式计算节点是由边界的已知节点，计算内部节点，计算时需要联立大型方程组，该方程组可以用迭代法求解。4.2算法步骤 主要步骤：(1)首先取定，对求解区域划分网格，按照网格定义矩阵v1，使得矩阵里面每个元素对应求解区域中的每个节点；(2)由边界条件定义v1矩阵中边界元素的值，其余元素定义为零；(3)定义与v1同型的零矩阵v2；(4)用五点格式公式通过矩阵v1迭代计算矩阵v2，迭代精度为0.1；(5)画图。4.3用matlab编写源程序 Matlab程序源代码：定义fivepoint.m文件function[peuxyk]=fivepoint(h,m,n,kmax,ep)%h步长%m,n为x,y方向的网格数，例如（2-0）/0.01=200;%kmax为最大迭代次数%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));%求区域中方程的每个节点endendt=zeros(n-1,m-1);for(k=1:kmax)for(i=2:n)for(j=2:m)%五点差分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;%达到最大迭代次数，终止循环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五点差分法应用举例 题目： 给定如下Laplace方程（poisson方程的特殊情况）的定解问题利用椭圆型方程的五点格式，计算该问题的近似解，并且画出近似解的图形。解：利用Matlab程序计算：命令窗口输入：[peuxyk]=fivepoint(0.025,80,40,10000,1e-11);k=3952;figuresurf(x,y,u);xlswrite('D:\define.xls',e,'data1');%保存数据xlabel('x');ylabel('y');zlabel('u');title('五点差分法解椭圆型偏微分方程');figuresurf(x,y,e);xlswrite('D:\wucha.xls',e,'data2');%保存数据title('步长为0.025时的误差');figuresurf(x,y,p);xlswrite('D:\jinshi.xls',e,'data32');%保存数据title('精确解曲面图');结果输出： 精确解近似解误差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