heat1d.m

应用数值方法解算一维热传导问题

function [x,t,u,Frames] = heat1d
% Heat diffusion in a wire

% Setup parimeter
L = 1.0;
T = 150;
K = 0.001;
n = 150;
dt = T/n;
l =20;
dx = L/l;

u = zeros(l+1,n+1);

% Create mesh
t = (0:n)*dt;
x = (0:l)*dx;

% Set i.c. and b.c.s
u(:,1) = u0(x)';

u(1,:) = bcLeft;

u(l+1,:) = bcRight;

alpha = K*dt/dx^2;

for k=1:n
    for i=2:l
        u(i,k+1) = alpha*u(i-1,k) + (1 - 2*alpha)*u(i,k) + alpha*u(i+1,k)+ dt*f;
    end
end

for k=1:n+1
    plot(x,u(:,k));
    axis([0 1 0 10]);
    Frames(k) = getframe;
end

% Some functions

function f = f
%f = 0.1;
f = 0.0;

function u0 = u0(x)
%u0 = 0.0;
for i = 1:length(x)
   if x(i) < 0.5
       u0(i) = 5.0;
  else
       u0(i) = 0.0;
   end
end

function bcLeft = bcLeft
bcLeft = 0.0;

function bcRight = bcRight
bcRight = 0.0;