六点对称格式_第1页
六点对称格式_第2页
六点对称格式_第3页
六点对称格式_第4页
全文预览已结束

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

1、利用差分法求解下列问题:方向,方向上。在观察数值解与精确解的误差。一、 算法描述1. 网格剖分取,方向,方向上。2. 差分格式选用六点对称格式,即Crank-Nicolson格式,是将向前差分格式和向后差分格式做算术平均得到的一种差分格式。,得到矩阵形式:3. 初边值处理直接给出。在本问题中,边值都为0,所以程序简化了很多,编程时和都直接取为0。二、 程序%*清理缓存*clear all;clc;dx=0.1; %空间步长dt=0.01; %时间步长t0=0;tn=0.25; %时间起点;时间终点x0=0;xn=1; %空间起点;空间终点N=(xn-x0)/dx; %空间网格数T=(tn-t0

2、)/dt; %时间网格数for i=1:N+1 x(i)=(i-1)*dx;end%*边界条件*%u(0,t)=0,u(1,t)=0;%可知在边界上的值始终为0,后面的时间每层迭代都直接取为0.%*初始条件*u=sin(pi*x)+x.*(1-x);r=dt/dx2;%*储存三对角元素*%为减小储存量,只储存三对角元素%A=zeros(3,N-1);for i=1:N-1 A(1,i)=-r/2;A(1,N-1)=0; A(2,i)=1+r; A(3,1)=0;A(3,i)=-r/2;endA=mylu(A); %LU分解a=r/2;b=1-r;c=r/2;for j=1:T %*计算方程组右

3、边的值* d(1)=b*u(2)+c*u(3)+0.02; for i=2:N-2 d(i)=a*u(i)+b*u(i+1)+c*u(i+2)+0.02; end d(N-1)=a*u(N-1)+b*u(N)+0.02; % d=myzhuigan(A,d); %追赶法求解 u=0 d 0; %两个0是给出的边界条件endplot(x,u,o) %绘制数值解图形holdu_xt=exp(-pi*pi*T*dt).*sin(pi.*x)+x.*(1-x);plot(x,u_xt,r) %绘制精确解图形 e=u_xt-u %误差%*over*子程序1:function LU=mylu(A)%这里的

4、LU分解是为了分解只储存三对角元素的矩阵% for example% A=2 2 2 2 0;2 1 1 1 1;0 -1 -1 -1 -1;% result:% LU =% 2.0000 2.0000 2.0000 2.0000 0% 2.0000 2.0000 2.0000 2.0000 2.0000% 0 -0.5000 -0.5000 -0.5000 -0.5000%n=size(A,2); %矩阵的列数for i=2:n A(3,i)=A(3,i)/A(2,i-1); A(2,i)=A(2,i)-A(3,i)*A(1,i-1);endLU=A;子程序2:function x=myzhuigan(A,d)n=length(d);%*追的过程*for i=2:n d(i)=d(i)-A(3,i)*d(i-1);end%*赶的过程*d(n)=d(n)

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论