matlab编写抛物线方程的最简差分格式

% function chap3_parabolic_1D
% test the finite difference method for 1D parabolic equation
% u_t = a*u_xx + f(x)，0<x<L，0<t<T,
% u(x,0) = phi(x)，0<x<L,
% u(0, t) = 0, u(0,t) = 0,





L =1;
T= 1;
a =1;
N = 500;% number of subintervals of the partition of time
J = 10;% number of subintervals of the partition of space
dt = T/N;
dx = L/J;
x =(0:dx:L)';
r= a*dt/dx^2


f = @f1;%y = zeros(size(x));
init = @init_val_1;%sin(pi*x);
u_0 = init(x);%9*1
u_0 = u_0(2:end-1);
f_n = f(x);
f_n= f_n(2:end-1);
f_n = f(x);
f_n = f_n(2: end-1);
%  % forward scheme
% u_fw_T = forward(x,r,dx,dt,J,N,u_0,f_n);
% u_fw_T = [0 ;u_fw_T;0];
% backward scheme
% u_bw_T = backward(x,r,dx,dt,J,N,u_0,f_n);
% u_bw_T = [0;u_bw_T;0];
% % crank-nicolson scheme
% u_cn_T = crank_nicolson(x, r,dx,dt,J,N,u_0,f_n);
% u_cn_T = [0;u_cn_T;0];
% % richardson scheme
u_ri_T = richardson (x, r, dx, dt,J, N, u_0, f_n);
u_ri_T = [0;u_ri_T;0];

u_ex=exact_sol(x, N*dt) ; 
%y = exp(-pi^2*t).*sin(pi*x);
    
% figure(1)
% plot (x, u_fw_T, 'g',x, u_ex,'r')
% title(' forward');
% figure(2)
% plot(x, u_bw_T, 'g', x, u_ex,'r')
% title(' backward') ;
% figure(3)
% plot(x, u_cn_T, 'g', x, u_ex, 'r')
% title(' Crank-Nicolson') ;
figure(4)
plot(x,u_ri_T, 'g', x, u_ex, 'r')
title('Richardson') ;
% end    
    
function u_n1 = forward(x, r, dx, dt,J, N, u_n,f_n)
d_0 = (1-2*r)*ones(J-1,1);
d_1 = r*ones(J-1,1);
d_m1 = d_1;
C = spdiags([d_m1 d_0 d_1],[-1 0 1],J-1,J-1);
for n = 1:N
    u_n1 = C*u_n + dt*f_n;
    u_n = u_n1;
    figure(1)
    plot(x,[0;u_n;0],'g');
    title('forward') ;
%pause;
end
end  
        
function u_n1 = backward(x,r,dx,dt,J,N,u_n,f_n)
d_0 = (1+2*r)*ones(J-1,1);
d_1 = -r*ones(J-1,1);
d_ml = d_1;
C = spdiags([d_ml d_0 d_1],[-1 0 1],J-1,J-1);
for n = 1:N
    b_n = u_n + dt*f_n;
    u_n1 = C\b_n;
    u_n = u_n1;
    figure(2)
    plot(x,[0;u_n;0],'g');
    title('backward') ;
%pause;
end
end

function u_n1 = crank_nicolson(x,r, dx,dt,J,N,u_n,f_n)
d_0 = (1+r)*ones(J-1,1);
d_1 = -0.5*r*ones(J-1,1);
d_m1 = d_1;
Cl = spdiags([d_m1 d_0 d_1],[-1 0 1],J-1,J-1);
d_0 = (1-r)*ones(J-1,1);
d_1 = 0.5*r*ones(J-1,1) ;
d_m1 = d_1;
Cr = spdiags([d_m1 d_0 d_1],[-1 0 1],J-1,J-1);
for n = 1:N
b_n = Cr*u_n + dt*f_n;
u_n1 = Cl\b_n;
u_n = u_n1;
figure(3)
plot(x,[0;u_n;0],'g');
title('Crank-Nicolson');
%pause;
end
end

function u_n1 =richardson(x, r, dx, dt,J,N,u_n,f_n)
d_0 = (1-2*r)*ones(J-1,1);
d_1 = r*ones(J-1,1);
d_m1 = d_1;
C = spdiags([d_m1 d_0 d_1],[-1 0 1],J-1,J-1);
u_n1 = C*u_n + dt*f_n;
u_nm1 = u_n;
u_n = u_n1;%向前差分先算一部

d_0 = -4*r*ones(J-1,1);%122
d_1 = 2*r*ones(J-1,1);
d_m1 = d_1;
C = spdiags([d_m1 d_0 d_1],[-1 0 1],J-1,J-1);
for n = 1:N
    u_n1 = C*u_n + u_nm1 + dt*f_n;
    u_nml = u_n;
    u_n = u_n1;
    figure(4)
    plot(x,[0 ;u_n;0], 'g');
    title('Richardson');
%pause;
end
end

function y = f1(x)%134
y = zeros(size(x));
end

function y = init_val_1(x)
y = sin(pi*x);
end

function y = exact_sol(x, t)
y = exp(-pi^2*t).*sin(pi*x);
end

u_n=u_0;
 d_0 = (1-2*r)*ones(J-1,1);
d_1 = r*ones(J-1,1);
d_m1 = d_1;
C = spdiags([d_m1 d_0 d_1],[-1 0 1],J-1,J-1);
for n = 1:N
    u_n1 = C*u_n + dt*f_n;
    u_n = u_n1;
    figure(1)
    plot(x,[0;u_n;0],'g');
    title('forward') ;
end