六点对称格式.docVIP

利用差分法求解下列问题：方向，方向上。在观察数值解与精确解的误差。算法描述网格剖分取，方向，方向上。差分格式选用六点对称格式，即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)/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/dx^2;%****************储存三对角元素***************************%%为减小储存量，只储存三对角元素%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 %**********计算方程组右边的值************** 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)%这里的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)/A(2,n);for i=n-1:-1:1 d(i)=(d(i)-A(1,i)*d(i+1))/A(2,i);endx=d;三、结果：