clear all
close all
clc
%for example-11，用隐式差分格式求解抛物型偏微分方程

N=1000;
h=1/N;
x=0:h:1;    %空间步长
y=@(x) exp(-20*(x-0.5).^2)-exp(-20*(x-1.5).^2)-exp(-20*(x+0.5).^2); %初始条件

delta_t=0.0001; %时间步长
t=0:delta_t:0.1;%时间上的迭代次数，决定了最终得到哪个时刻的解
y_n=zeros(N+1,length(t));
y_n(1,1)=0;  %边界条件
y_n(N+1,1)=0; %边界条件
%yy=zeros(2:N,1); 

%精确解
y_r=@(xx,t)(1/sqrt(1+80*t).*(exp(-20*(xx-0.5).^2/(1+80*t))...
    -exp(-20*(xx-1.5).^2/(1+80*t))-exp(-20*(xx+0.5).^2/(1+80*t))));

%================数值解（隐式差分格式）==========================
rr=delta_t/h^2

A=zeros(N-1,N-1); %给系数矩阵赋值
A(1,1)=1+2*rr; A(1,2)=-rr; %将系数矩阵第一行和最后一行单独赋值
A(N-1,N-1)=1+2*rr; A(N-1,N-2)=-rr;
for i=2:N-2       %对系数矩阵中中间行赋值
    A(i,i)=1+2*rr;
    A(i,i-1)=-rr;
    A(i,i+1)=-rr;
end
B(1:N-1,1)=y(h:h:(N-1)*h);

AA=A^-1; %矩阵求逆
yy=zeros(N-1,1);
for ii=1:length(t)-1
%给方程右边的列矩阵赋值
yy=AA*B;   %采用矩阵乘法获得解的列矩阵
B=yy;      
%获得新的方程右边的列矩阵的值，对于此特殊问题，S=0，边值为0，B与yy相等，但不具有普遍性
%y_n(N+1,i)=yy;
end

y_n(2:N,length(t))=yy(1:N-1,1);
%注意解析结果的最后一个值不为0
error=y_n(:,length(t))-y_r(x,delta_t*(length(t)-1))';
%注意：y_(:,length(t))以列存储，y_r(x,delta_t*N_t)以行数组存储，所以需要转置运算
figure(1),plot(x,y_n(:,length(t)),'-r',x,y_r(x,delta_t*(length(t)-1)),'b',x,y(x),'-k')
legend('数值结果','解析结果','初始状态')
figure(2),plot(x,error)
legend('数值结果与解析结果间的误差')