加权六点格式求解抛物型偏微分方程（Matlab实现）

加权六点格式求解抛物型偏微分方程（Matlab实现）

clear
clc
at=[0,1];
ax=[0,1];
t1=0.01;
x1=0.01;
theta=0.5;
th=at(1):t1:at(2);
xh=ax(1):x1:ax(2);
n1=length(xh);
n2=length(th);
u=zeros(n2,n1);
syms x t;
bot=20;%时间边界
boxr=0;%空间左边界
boxl=100;%空间右边界
%下面求边界值
bt=zeros(1,n1);
for i=1:n1
   bt(i)=subs(bot,x,xh(i));
   u(1,i)=bt(i);
end
bxr=zeros(1,n2);
bxl=zeros(1,n2);
for i=1:n2
   bxr(i)=subs(boxr,t,th(i));
   bxl(i)=subs(boxl,t,th(i));
   u(i,1)=bxr(i);
   u(i,n1)=bxl(i);
end
p=t1/(x1)^2;
A=zeros(n1-2,n1-2);
for i=1:n1-2
   A(i,i)=1+2*p*theta; 
end
for i=2:n1-2
   A(i,i-1)=-p*theta; 
end
for i=1:n1-3
   A(i,i+1)=-p*theta; 
end
c=zeros(n1-2,1);
for i=2:n2
   b=zeros(n1-2,1);
   for j=2:n1-1
      b(j-1)=p*(1-theta)*u(i-1,j-1)+(1-2*p*(1-theta))*u(i-1,j)+p*(1-theta)*u(i-1,j+1); 
   end
   b(1)=b(1)+p*theta*u(i,1);
   b(n1-2)=b(n1-2)+p*theta*u(i,n1);
   aa=A\b;
   for j=1:n1-2
    u(i,j+1)=aa(j); 
   end
  % plot(xh,u(i,:));
  % hold on
end
uu=zeros(n2,n1);
for i=1:n2
    for j=1:n1
        uu(i,j)=exp(-pi^2*th(i))*sin(pi*xh(j));
    end
end
%mesh(xh,th,u);%数值解图
%mesh(xh,th,uu);%真实值图
u2=uu-u;
%mesh(xh,th,u2);%误差图;