六点对称格式

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)