微分方程数值解实验.doc

微分方程数值解课程设计报告 班级：_ 姓名： _ 学号： _ 成绩： 2017年 6月 21 日 目 录一、摘 要1二、常微分方程数值解22.1 4阶Runge-Kutta法和Adams4阶外插法的基本思路22.2 算法流程图22.3 用matlab编写源程序22.4 常微分方程数值解法应用举例4三、 常系数扩散方程的经典差分格式63.1 有限差分法的基本思路63.2 算法流程图73.3 用matlab编写源程序73.4 有限差分法应用举例8四、 椭圆型方程的五点差分格式104.1 五点差分法的基本思路104.2 算法流程图114.3 用matlab编写源程序114.4 五点差分法应用举例12五、 自我总结16六、参考文献16一、摘 要自然界与工程技术中的很多现象，可以归结为微分方程定解问题。其中，常微分方程求解是微分方程的重要基础内容。但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、RungeKutta方法、Adams法以及椭圆型方程、抛物型方程的有限差分方法等，通过具体的算例，结合MATLAB求解画图，初步给出了一般常微分方程数值解法的求解过程。同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用范围有一个直观的感受。关键词：微分方程数值解、MATLAB 二、常微分方程数值解2.1 基本思路常微分方程数值解法(numerical methods forordinary differential equations)计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题，实用上许多很有价值的常微分方程的解不能用初等函数来表示，常常需要求其数值解.所谓数值解，是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.常微分方程初值问题的数值解法是求方程（1）的解在点列上的近似值，这里是到的步长，一般略去下标记为。 （1） 经典的方法是一个四阶的方法，它的计算公式是： （2）方法的优点是：单步法、精度高，计算过程便于改变步长，缺点是计算量较大，每前进一步需要计算四次函数值。在用龙格库塔方法时,要注意的选择要合适，太大，会使计算量加大，太小，较大，可能会使误差增大。因此选择合适的很重要。我们要在考虑精度的基础上，选择合适的。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文件-function dy = 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);grid on;legend(x关于t 的变化关系图,1);xlabel(t,FontSize,14);ylabel(x,FontSize,14);subplot(2,2,2);plot(t,y(:,2);grid on;legend(y关于t 的变化关系图,1);xlabel(t,FontSize,14);ylabel(y,FontSize,14);subplot(2,2,3)plot(t,y(:,3);grid on;legend(z关于t 的变化关系图,1);xlabel(t,FontSize,14);ylabel(z,FontSize,14);subplot(2,2,4)plot3(y(:,1),y(:,2),y(:,3);grid on;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;for i=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=%fy=%fz=%fn,x,y,z);%显示迭代值end suptitle(Adams4阶外插法) %总标题subplot(2,2,1);plot(t,a);grid on; legend(x关于t 的变化关系图,1);xlabel(t,FontSize,14);ylabel(x,FontSize,14);subplot(2,2,2);plot(t,b);grid on; legend(y关于t 的变化关系图,1);xlabel(t,FontSize,14);ylabel(y,FontSize,14);subplot(2,2,3);plot(t,c);grid on; legend(z关于t 的变化关系图,1);xlabel(t,FontSize,14);ylabel(z,FontSize,14);subplot(2,2,4)plot3(a,b,c);grid on; 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=theta l=taoh=0.1; %h :空间步长t=l*h*h; x=0.1:0.1:0.9; %x坐标取值for n=1:9 u(n)=sin(x(n)*pi); %初始值endu=u; v=zeros(9,1); %初始化v矩阵for n=1:k A=zeros(9,9); A(1,1)=1+2*q*t/h/h; for i=2:9 A(i,i)=1+2*q*t/h/h; A(i-1,i)=-q*t/h/h; A(i,i-1)=-q*t/h/h; end B=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); for i=2:8 B(i,1)=(1-q)*t/h/h*(u(i-1)-2*u(i)+u(i+1)+u(i); %加权隐式格式 end v=inv(A)*B; %求解矩阵v w=u; u=v; v=w; n=n+1; endtt=k*t; uu=zeros(9,1); %初始化uu矩阵for i=1:9 uu(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处的准确值为：%fn,uu(4)u,tt,uu=parabola(100,0.25,0.1);fprintf(=0.25，=0.1 时 0.4处的数值值为：%fn,u(4)fprintf(误差为：%fnnn,abs(uu(4)-u(4)x=0.1:0.1:0.9;subplot(2,2,1);plot(x,u,x,uu,o);grid on;title(=0.25，=0.1);u,tt,uu=parabola(100,0.25,0.5);fprintf(=0.25，=0.5 时 0.4处的准确值为：%fn,uu(4)u,tt,uu=parabola(100,0.25,0.5);fprintf(=0.25，=0.5 时 0.4处的数值值为：%fn,u(4)fprintf(误差为：%fnnn,abs(uu(4)-u(4)x=0.1:0.1:0.9;subplot(2,2,2);plot(x,u,x,uu,o);grid on;title(=0.25，=0.5);u,tt,uu=parabola(100,0.5,0.1);fprintf(=0.5，=0.1 时 0.4处的准确值为：%fn,uu(4)u,tt,uu=parabola(100,0.5,0.1);fprintf(=0.5，=0.1 时 0.4处的数值值为：%fn,u(4)fprintf(误差为：%fnnn,abs(uu(4)-u(4)x=0.1:0.1:0.9;subplot(2,2,3);plot(x,u,x,uu,o);grid on;title(=0.5，=0.1);u,tt,uu=parabola(100,0.5,0.5);fprintf(=0.5，=0.5 时 0.4处的准确值为：%fn,uu(4)u,tt,uu=parabola(100,0.5,0.5);fprintf(=0.5，=0.5 时 0.4处的数值值为：%fn,u(4)fprintf(误差为：%fnnn,abs(uu(4)-u(4)x=0.1:0.1:0.9;subplot(2,2,4);plot(x,u,x,uu,o);grid on;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 p e u x y k= fivepoint (h,m,n,kmax,ep)%h步长%m,n为x,y方向的网格数，例如（2-0）/0.01=200;%kmax为最大迭代次数%e为误差，p为精确解syms temp;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; end end t(i,j)=sqrt(t(i,j); if(kkmax) break;%达到最大迭代次数，终止循环 end if(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);%误差 endend-4.4 五点差分法应用举例题目：给定如下Laplace方程（poisson方程的特殊情况）的定解问题利用椭圆型方程的五点格式，计算该问题的近似解，并且画出近似解的图形。解： 利用Matlab程序计算：命令窗口输入：p e u x y k=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