MATLAB代码
clear
a = 1;
dx = 0.01;%dx要大
x = 0:dx:1;
dt = 0.00005;%dt要小
t = 0:dt:1;
u = zeros(length(x),length(t));
% 初始条件,x的函数
u(:,1) = sin(pi*x);
% 边界条件,t的函数
m1 = 1 + 0.0*sin(t);
m2 = 1 - 0.2*sin(10*t);
% 构造A矩阵
A = -2*eye(length(x)) + diag(ones(1,length(x)-1),1) + diag(ones(1,length(x)-1),-1);
for n = 1:length(t) - 1
u(:,n+1) = u(:,n) + a^2*dt/dx^2*A*u(:,n);
u(1,n+1) = m1(n+1);
u(end,n+1) = m2(n+1);
% 温度随时间传导的动画
plot(x,u(:,n+1))
axis([x(1) x(end) 0 1])
getframe;
end
% [T,X] = meshgrid(t,x);
% surf(X,T,u);
% shading interp