初边值问题matlab代码,给定初边值条件的椭圆形方程数值解法

function fan=calculate(step)

%% 设置沿t方向与x方向的选取采样点的步长(自由设置)

h=step;%x方向

r=0.1;

k=r*h.^2;%t方向

%% 变化范围、采样点数目、采样点的计算

T=0.001;

X=1;

n=T/k;%t采样数目

m=1/h;%x方向采样数目 = 区间长度 / 采样步长

t=0:k:T;% 采样点的坐标存储在数组 t[]、x[]中

x=0:h:X;%初值点、边界上的点也存在这两个数组中

tt=t(2:n+1); %#okxx=x(2:m);%数值差分预测区，内部点

t0=sin(pi*(0:h:1))';%初值点，边值点均为0，不需要存储

%% 以精确解并作出图像

exact=zeros(1,m-1);%记录精确解

for i=1:1:m-1

exact(i)=exp(1).^(-pi^2*T)*sin(pi*xx(i));

end

figure;

plot(xx,exact,'g');%仅画出t=T时的图像，而且仅画出内部需要近似的点

hold on

%% 使用古典隐式格式求数值解、追赶法求解方式是逐层求解

A=eye(m-1);%每一层矩阵A相同，事先求出

for i=1:1:m-1

A(i,i)=1+2*r;

if i~=1&&i~=m-1

A(i,i+1)=-r;

A(i,i-1)=-r;

end

end

A(1,2)=-r;%补充两点

A(m-1,m-2)=-r;

lastpoint=t0(2:m);%上一层节点

for j=1:1:n

newpoint=A\l