微分方程数值解实验报告.docx
- 文档编号:12203685
- 上传时间:2023-04-17
- 格式:DOCX
- 页数:17
- 大小:563.58KB
微分方程数值解实验报告.docx
《微分方程数值解实验报告.docx》由会员分享,可在线阅读,更多相关《微分方程数值解实验报告.docx(17页珍藏版)》请在冰豆网上搜索。
微分方程数值解实验报告
微分方程数值解法
课程设计报告
班级:
姓名:
学号:
成绩:
2017年6月21日
摘要
自然界与工程技术中的很多现象,可以归结为微分方程定解问题。
其中,常微分方程求解是微分方程的重要基础内容。
但是,对于许多的微分方程,往往很难得到甚至不存在精确的解析表达式,这时候,数值解提供了一个很好的解决思路。
,针对于此,本文对常微分方程数值解法进行了简单研究,主要讨论了一些常用的数值解法,如欧拉法、改进的欧拉法、Runge—Kutta方法、Adams法以及椭圆型方程、抛物型方程的有限差分方法等,通过具体的算例,结合MATLAB求解画图,初步给出了一般常微分方程数值解法的求解过程。
同时,通过对各种方法的误差分析,让大家对各种方法的特点和适用范围有一个直观的感受。
关键词:
微分方程数值解、MATLAB
目录
摘要2
目录3
第一章常微分方程数值解法的基本思想与原理4
1.1常微分方程数值解法的基本思路4
1.2用matlab编写源程序4
1.3常微分方程数值解法应用举例及结果5
第二章常系数扩散方程的经典差分格式的基本思想与原理6
2.1常系数扩散方程的经典差分格式的基本思路6
2.2用matlab编写源程序7
2.3常系数扩散方程的经典差分格式的应用举例及结果8
第三章椭圆型方程的五点差分格式的基本思想与原理10
3.1椭圆型方程的五点差分格式的基本思路10
3.2用matlab编写源程序10
3.3椭圆型方程的五点差分格式的应用举例及结果12
第四章总结12
参考文献12
第一章常微分方程数值解法的基本思想与原理
1.1常微分方程数值解法的基本思路
常微分方程数值解法(numericalmethodsforordinarydifferentialequations)计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题,实用上许多很有价值的常微分方程的解不能用初等函数来表示,常常需要求其数值解.所谓数值解,是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.
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:
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,:
)=[xyz];
end
plot3(Y(:
1),Y(:
2),Y(:
3));
1.3常微分方程数值解法的应用举例及结果
应用举例:
a=10,b=8/3,0 , ,r-1)和 (- ,- ,r-1),一个稳定的不动点0=(0,0,0),当r>24.74时,c和 都变成不稳定的,此时存在混沌和奇怪吸引子。 运行结果: 龙格库塔法: 欧拉法: 第二章常系数扩散方程的经典差分格式的基本思想与原理 2.1常系数扩散方程的经典差分格式的基本思路 用有限差分法解常系数扩散方程 有加权隐式差分格式 其中 ,当 时为Crank-Nicolson格式,当 时为向后差分格式,当 时为向前差分格式。 加权隐式格式稳定的条件是 ,当 ,无限制, 当 。 加权隐式格式是两层隐式格式,用第n层计算第n+1层节点值的时候,要解线性方程组。 2.2用matlab编写源程序 M文件: functionM=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); fori=(n-1): (-1): 1 a(i+1)=a(i); end %将a矩阵和n对应 beta (1)=c (1)/b (1); fori=2: (n-1) beta(i)=c(i)/(b(i)-a(i)*beta(i-1)); end y (1)=f (1)/b (1); fori=2: n y(i)=(f(i)-a(i)*y(i-1))/(b(i)-a(i)*beta(i-1)); end M(n)=y(n); fori=(n-1): (-1): 1 M(i)=y(i)-beta(i)*M(i+1); end end M文件: functionoutput=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/h^2; P1=1/P+1; P2=1/P-1; u=zeros(k,n);%初始条件 u(1,: )=exp(x);%求A矩阵的对角元素d d=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矩阵的对角元素下面一行元素e e=-ones(1,n-1); e(1,n-1)=-b2;%求A矩阵的对角元素上面一行元素f f=-ones(1,n-1); f(1,1)=-b1; R=zeros(k,n);%求R %追赶法求解 fori=2: k R(i,1)=(b1*P2-h*a1)*u(i-1,1)+b1*u(i-1,2)+2*h*c1; forj=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([0a00t_max]);pause(0.1) end output=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));%使图向变为黑色 figure surf(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常系数扩散方程的经典差分格式的应用举例及结果 应用举例: 考虑常系数扩散方程的初边值问题 其中 ,取 为时间步长, 为网格比,对不同的时间步长 ,计算当 时初边值问题的解u(0.4,0.4),并且与精确解比较,分析比较结果。 运行结果: 第三章椭圆型方程的五点差分格式的基本思想与原理 3.1椭圆型方程的五点差分格式的基本思路 对Laplace方程的第一边值问题 利用taylor展开可得逼近它的五点差分格式的差分逼近 其中 分别为 轴和 轴步长,边界条件可以由 离散可得,当 时有 。 注意五点格式计算节点是由边界的已知节点,计算内部节点,计算时需要联立大型方程组,该方程组可以用迭代法求解。 3.2用matlab编写源程序 M文件: function[peuxyk]=wudianchafenfa(h,m,n,kmax,ep) %g-s迭代法解五点差分法问题%kmax为最大迭代次数 %m,n为x,y方向的网格数,例如(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+1 u(i,1)=sin(pi*y(i)); u(i,m+1)=exp (1)*exp (1)*sin(pi*y(i)); end fori=1: n forj=1: m f(i,j)=(pi*pi-1)*exp(x(j))*sin(pi*y(i)); end end t=zeros(n-1,m-1); fork=1: kmax fori=2: n forj=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)); ifk>kmax break; end ifmax(max(t)) break; end end fori=1: n+1 forj=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))); end end 运行程序: [peuxyk]=wudianchafenfa(0.1,20,10,10000,1e-6); surf(x,y,u); xlabel('x');ylabel('y');zlabel('u'); title('五点差分法解椭圆型偏微分方程'); 3.3椭圆型方程的五点差分格式的应用举例及结果 应用举例: 给定如下Laplace方程(poisson方程的特殊情况)的定值问题 利用椭圆型方程的五点格式,计算该问题的近似解,并且画出近似解的图形 运行结果: 第四章总结 通过本次课程设计,在实践操作中我学会了很多,同时也认识到了自身的不足之处。 课程设计是对在课本上学到的知识的一种实际应用,由于我对课本知识掌握的不牢固,所以在实践操作中遇到了很多问题。 自课程设计中我深刻的认识到了理论知识的重要性,只要有牢固的掌握理论知识,才能在实践操作中以理论知识为基础,得心应手地处理各种问题。 参考文献 《微分方程数值解法》(第四版)李荣华,刘播高等教育出版社
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 数值 实验 报告
![提示](https://static.bdocx.com/images/bang_tan.gif)