suanli4_3_1.m

《偏微分方程数值解法》孙志忠版中算例的一些MATLAB程序

% compace 差分格式
%算例4.3.1
%%%%%%%网格剖分%%%%%%%%%
h = 1/10;  				%%空间步长
a = 1;
dt = 1/100;			   %时间步长
s = a * dt/h;          %步长比
T = 1;	  			   %  时间长度
n = ceil(T / dt);	   % 时间网格
X0 = 1;				   % 空间长度
m = ceil(X0 / h);	   %空间网格
u = zeros(n + 1, m + 1);	%%数解初始化u(tk,xi)
% %初始值及边界值
x=0:h:X0;
phi=exp(x);
u(1,:)=phi;         %u的第一行表示初始状态
t=0:dt:T;
u(:,1)=exp(t');         %左界
u(:,m+1)=exp(1+t');
%%%%%%%%%%%%%%%
f=zeros(n+1,m+1);
for k=1:n+1
    for i=1:m+1
       f(k,i)=0;
   end
end
%%%%%%%%%%%%%%%
psi=exp(x);
for i = 2:m
    u(2,i) = dt * psi(i) + s^2/2 * ( phi(i-1) + phi(i+1) ) ...
              + ( 1 - s^2 ) * phi(i) + dt^2/2 * f(1,i);
end
%%%%%%%%%%%%%%%
A=diag( (5/6+s^2 )*ones( m-1,1 ) )+diag( (1/12-1/2*s^2) * ones (m-2,1),1 )+diag( (1/12-1/2*s^2) * ones (m-2,1),-1 );
B=diag( -(5/6+s^2 )*ones( m-1,1 ) )+diag( (1/2*s^2-1/12) * ones (m-2,1),1 )+diag( (1/2*s^2-1/12) * ones (m-2,1),-1 );
C=diag(5/3*ones(m-1,1))+diag(1/6*ones(m-2,1),1)+diag(1/6*ones(m-2,1),-1);
for k=2:n
    for i=2:m
        F(i-1)=1/12*dt^2*( f(i-1)+10*f(i)+f(i+1) );
    end
    F(1)=F(1) + (1/2*s^2-1/12) * ( u(k+1,1) + u(k-1,1) )+1/6*u(k,1);
    F(m-1)=F(m-1) + (1/2*s^2-1/12) * ( u(k+1,m+1) + u(k-1,m+1) )+1/6*u(k,m+1);
    y = A \ ( B * u(k-1,2:m )' + F' + C* u( k,2:m )' );
    for i=1:m-1
        u(k+1,i+1)=y(i);
    end
end
%%%%%%%%%%%%%%%%给出精确解
ana = zeros(n + 1, m + 1);
for k = 1 : n+1
	for i = 1 : m+1
         ana( k, i ) = exp ( (i-1) *h + (k-1)*dt);
   end;
end;
%%%%%%%%%%%%%%%给出在T=1时的最大误差
max( abs ( u (n+1,:) -ana (n+1,:) ) )
%图示化
[tt,xx]= meshgrid(t,x);
xx=xx';   tt=tt';
mesh(tt,xx,ana);
title('analytical solution');
xlabel(sprintf('steps in time domain\nfrom 0 to 1'));
ylabel(sprintf('steps in x-domain\nfrom 0 to 1'));
zlabel('U(x,y)');
pause
mesh(tt,xx,u);
title('numerical solution');
xlabel(sprintf('steps in time domain\nfrom 0 to 1'));
ylabel(sprintf('steps in x-domain\nfrom 0 to 1'));
zlabel('u(xi,tk)');
pause
mesh(tt,xx,u-ana);
title('diffeence beween analical and numeical soluion');
xlabel(sprintf('steps in time- domain\nfrom 0 to 1'));
ylabel(sprintf('steps in x domain\nfrom 0 to 1'));
zlabel('u-ana');
pause
plot(x,u(n+1,:),'-or',x,ana(n+1,:))
title('显格式数值解与精确解的比较')
legend('uh(x,1)','u(x,1)')