.微分方程数值解法微分方程数值解法课程设计报告课程设计报告班级：班级：_姓名：姓名： _学号：学号：_成绩：成绩：2017 年 6 月 21 日.摘要自然界与工程技术中的很多现象，可以归结为微分方程定解问题。其中，常微分方程求解是微分方程的重要基础内容。但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。 ，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、 RungeKutta 方法、Adams 法以及椭圆型方程、抛物型方程的有限差分方法等，通过具体的算例，结合 MATLAB 求解画图，初步给出了一般常微分方程数值解法的求解过程。同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用范围有一个直观的感受。关键词：微分方程数值解、关键词：微分方程数值解、MATLABMATLAB.目录摘摘 要要 . 2 2目目 录录 . 3 3第一章第一章 常微分方程数值解法的基本思想与原理常微分方程数值解法的基本思想与原理 . 4 41.11.1 常微分方程数值解法的基本思路常微分方程数值解法的基本思路. 4 41.21.2 用用 matlabmatlab 编写源程序编写源程序 . 4 41.31.3 常微分方程数值解法应用举例及结果常微分方程数值解法应用举例及结果. . 5 5第二章第二章 常系数扩散方程的经典差分格式的基本思想与原理常系数扩散方程的经典差分格式的基本思想与原理 . .6 62.12.1 常系数扩散方程的经典差分格式的基本思路常系数扩散方程的经典差分格式的基本思路 . 6 62.22.2 用用 matlabmatlab 编写源程序编写源程序. . 7 7.2.32.3 常系数扩散方程的经典差分格式的应用举例及结果常系数扩散方程的经典差分格式的应用举例及结果. 8 8第三章第三章 椭圆型方程的五点差分格式的基本思想与原理椭圆型方程的五点差分格式的基本思想与原理 . .10103.13.1 椭圆型方程的五点差分格式的基本思路椭圆型方程的五点差分格式的基本思路 . 10103.23.2 用用 matlabmatlab 编写源程序编写源程序 . . 10103.33.3 椭圆型方程的五点差分格式的应用举例及结果椭圆型方程的五点差分格式的应用举例及结果. . 1212第四章第四章 总结总结 . . 1212参考文献参考文献 . 1212.第一章 常微分方程数值解法的基本思想与原理1.1 常微分方程数值解法的基本思路常微分方程数值解法(numerical methods forordinary differentialequations)计算数学的一个分支.是解常微分方程各类定解问题的数值方法 .现有的解析方法只能用于求解一些特殊类型的定解问题， 实用上许多很有价值的常微分方程的解不能用初等函数来表示，常常需要求其数值解.所谓数值解，是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.1.2 用 matlab 编写源程序龙格库塔法：M 文件:function dx=Lorenz(t,x)%r=28,sigma=10,b=8/3dx=-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=;for i=1:800 x1=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,:)=x y z;.endplot3(Y(:,1),Y(:,2),Y(:,3);1.3 常微分方程数值解法的应用举例及结果应用举例：dx(t)dt a(x(t) y(t)dy(t) rx(t) y(t) x(t)z(t)dtdz(t)dt x(t)y(t)bz(t)a=10,b=8/3,0r+,当 1r24.74 时，c 和c都变成不稳定的，此时存在混沌和奇怪吸引子。运行结果：龙格库塔法：.欧拉法：.第二章 常系数扩散方程的经典差分格式的基本思想与原理2.1 常系数扩散方程的经典差分格式的基本思路用有限差分法解常系数扩散方程u2u a2tx有加权隐式差分格式n1unjujnn1n11 unununj12ujuj1j12ujj1a(1)2hh2 0其中0 1，当1时为 Crank-Nicolson 格式，当1时为向后差分格式，当 02时为向前差分格式。加权隐式格式稳定的条件是2a11，当0 ，无限制，21 2.当11。2加权隐式格式是两层隐式格式，用第n 层计算第 n+1 层节点值的时候，要解线性方程组。2.2 用 matlab 编写源程序M 文件：function M = chase(a,b,c,f)% 追赶法求解三对角矩阵方程，Ax=f% a 是对角线下边一行的元素% b 是对角线元素% c 是对角线上边一行的元素% M 是求得的结果，以列向量形式保存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 对应beta(1) = c(1)/b(1);for i = 2:(n-1) beta(i) = c(i)/( b(i)-a(i)*beta(i-1) );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);endendM 文件：function output = diffuse_equation(a0,t_max,h,tao,D,a1,b1,c1,a2,b2,c2)% 一维扩散方程的有限差分法，采用隐式六点差分格式(Crank-Nicolson)% a0: x 的最大值% t:_max: t 的最大值% h: 空间步长% tao: 时间步长% D：扩散系数% a1,b1,c1 是（x=0）边界条件的系数；a2,b2,c2 是（x=a0）边界条件的系数x = 0:h:a0;n = length(x);t = 0:tao:t_max;.k = length(t);P = tao * D/h2;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%追赶法求解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(i,:); u(i,:) = M; plot(x,u(i,:); axis(0 a0 0 t_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);%使图向变为黑色figuresurf(X,T,u),xlabel(x),ylabel(t),zlabel(u),title(有限差分解);%colormap(gray(1);%使图向变为黑色运行程序：% 一维扩散方程的有限差分法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;.%调用扩散方程子函数求解u = diffuse_equation(a0,t_max,h,tao,D,a1,b1,c1,a2,b2,c2);2.3 常系数扩散方程的经典差分格式的应用举例及结果应用举例：考虑常系数扩散方程的初边值问题u2utx2,0 x 1,t 0,u(x,0) sin,0 x 1,u(0,t) u(1,t) 0，t 0n1unjujnn1n11 unununj1 2ujuj1j1 2ujj1 (1) 022hh其中0 1，取j 10,h 比， 对不同的时间步长(1,x1 jh( j 0,1, , j),为时间步长，2为网格10h1 11计算当 0,1,时初边值问题的解u(0.4,0.4),1,2,4,8)，4 2,2并且与精确解比较，分析比较结果。运行结果：.第三章 椭圆型方程的五点差分格式的基本思想与原理3.1 椭圆型方程的五点差分格式的基本思路对 Laplace 方程的第一边值问题2u2uu 22 0,x, yDxyux, yx, y,x, yD利用 taylor 展开可得逼近它的五点差分格式的差分逼近huijij,(xi, yj)Dh-hu ui1, j 2uijui1, j2ui, j1 2uijui, j1k2 0,(xi, yj)Dh其中h,k分别为x轴和y轴步长，边界条件可以由ux, yx, y离散可得，当(x, y)D时有ijxi, yj。注意五点格式计算节点是由边界的已知节点， 计算内部节点， 计算时需要联立大型方程组，该方程组可以用迭代法求解。3.2 用 matlab 编写源程序.M 文件：function p e u x y k=wudianchafenfa(h,m,n,kmax,ep)% g-s 迭代法解五点差分法问题 %kmax 为最大迭代次数 %m,n 为 x,y 方向的网格数，例如（2-0）/0.01=200;%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); u(i,m+1)=exp(1)*exp(1)*sin(pi*y(i);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运行程序：p e u x y k=wudianchafenfa(0.1,20,10,10000,1e-6); surf(x,y,u);. xlabel(x);ylabel(y);zlabel(u); title(五点差分法解椭圆型偏微分方程);3.3 椭圆型方程的五点差分格式的应用举例及结果应用举例：给定如下 Laplace 方程（poisson 方程的特殊情况）的定值问题u 0, 0 x 17,0 y 10,u(0, y) u(17, y) u(x,0) 0,u(x,10) 100.利用椭圆型方程的五点格式，计算该问题的近似解，并且画出近似解的图形运行结果：第四章 总结通过本次课程设计，在实践操作中我学会了很多，同时也认识到了自身的不足之处。课程设计是对在课本上学到的知识的一种实际应用， 由于我对课本知识掌握的不牢固， 所以在实践操作中遇到了很多问题。 自课程设计中我深刻的认识到了理论知识的重要性， 只要有牢固的掌握理论知识，才能在实践操作中以理论知识为基础，得心应手地处理各种问题。参考文献微分方程数值解法微分方程数值解法 （第四版）（第四版） 李荣华，刘播李荣华，刘播高等教育出版社高等教育出版社.