求解∂u/∂t=1/16 ∆u的六点对称格式，Du Fort-Frankel 格式，ADI格式，LOD格式matlab代码

以下是求解∂u/∂t=1/16 ∆u的六点对称格式、Du Fort-Frankel 格式、ADI格式、LOD格式的matlab代码：

六点对称格式：

L = 1; % length of the domain
N = 100; % number of grid points
h = L/(N-1); % grid spacing
dt = 0.01; % time step
T = 10; % final time
D = 1/16; % diffusion coefficient

x = linspace(0,L,N)';
u0 = exp(-50*(x-0.5).^2); % initial condition

% 6-point stencil coefficients
a = -D*dt/h^2;
b = 1 + 4*a;
c = -a;

% set up matrix A
A = sparse(N,N);
A(1,1) = 1;
A(N,N) = 1;
for i = 2:N-1
    A(i,i-1:i+1) = [a,c,b];
end

% time stepping
t = 0;
u = u0;
while t < T
    t = t + dt;
    u = A*u;
end

plot(x,u);

Du Fort-Frankel 格式：

L = 1; % length of the domain
N = 100; % number of grid points
h = L/(N-1); % grid spacing
dt = 0.01; % time step
T = 10; % final time
D = 1/16; % diffusion coefficient

x = linspace(0,L,N)';
u0 = exp(-50*(x-0.5).^2); % initial condition

% set up matrix A
a = -D*dt/h^2;
b = 1 + 2*a;
c = a;
A = sparse(N,N);
A(1,1) = 1;
A(N,N) = 1;
for i = 2:N-1
    A(i,i-1:i+1) = [a,c,b];
end

% Du Fort-Frankel scheme
B = inv(A)*(2*eye(N)-A);
C = 2*eye(N)+A;

% time stepping
t = 0;
u = u0;
while t < T
    t = t + dt;
    u = B*u;
    u = C*u;
end

plot(x,u);

ADI格式：

L = 1; % length of the domain
N = 100; % number of grid points
h = L/(N-1); % grid spacing
dt = 0.01; % time step
T = 10; % final time
D = 1/16; % diffusion coefficient

x = linspace(0,L,N)';
u0 = exp(-50*(x-0.5).^2); % initial condition

% set up matrices A1, A2, B1, B2
a = -D*dt/(2*h^2);
b = 1 + 2*a;
A1 = sparse(N,N);
A1(1,1) = 1;
A1(N,N) = 1;
for i = 2:N-1
    A1(i,i-1:i+1) = [a,a,b];
end
B1 = sparse(N,N);
B1(1,1) = 1;
B1(N,N) = 1;
for i = 2:N-1
    B1(i,i-1:i+1) = [-a,-a,b];
end
a = -D*dt/h^2;
b = 1 + 2*a;
A2 = sparse(N,N);
A2(1,1) = 1;
A2(N,N) = 1;
for i = 2:N-1
    A2(i,i-1:i+1) = [0,a,b];
end
B2 = sparse(N,N);
B2(1,1) = 1;
B2(N,N) = 1;
for i = 2:N-1
    B2(i,i-1:i+1) = [0,-a,b];
end

% time stepping
t = 0;
u = u0;
for n = 1:T/dt
    u = A1\(B1*u);
    u = A2\(B2*u);
end

plot(x,u);

LOD格式：

L = 1; % length of the domain
N = 100; % number of grid points
h = L/(N-1); % grid spacing
dt = 0.01; % time step
T = 10; % final time
D = 1/16; % diffusion coefficient

x = linspace(0,L,N)';
u0 = exp(-50*(x-0.5).^2); % initial condition

% set up matrices A1, A2, B1, B2
a = -D*dt/h^2;
b = 1 + 2*a;
A1 = sparse(N,N);
A1(1,1) = 1;
A1(N,N) = 1;
for i = 2:N-1
    A1(i,i-1:i+1) = [a,a,b];
end
B1 = sparse(N,N);
B1(1,1) = 1;
B1(N,N) = 1;
for i = 2:N-1
    B1(i,i-1:i+1) = [-a,-a,b];
end
a = -D*dt/(2*h^2);
b = 1 + 4*a;
c = -a;
A2 = sparse(N,N);
A2(1,1) = 1;
A2(N,N) = 1;
for i = 2:N-1
    A2(i,i-1:i+1) = [a,c,b];
end
B2 = sparse(N,N);
B2(1,1) = 1;
B2(N,N) = 1;
for i = 2:N-1
    B2(i,i-1:i+1) = [-a,-c,b];
end

% time stepping
t = 0;
u = u0;
for n = 1:T/dt
    u = A1\(B1*u);
    u = A2\(B2*u);
end

plot(x,u);