%一維的熱傳傳導(dǎo)方程的

%初邊值問題。

%北太天元數(shù)值計算軟件做心臟射頻消融術(shù)有關(guān)的數(shù)值模擬

%激光頭處的溫度作為時間的函數(shù)設(shè)置為如下的函數(shù)

% 這個函數(shù)和方波的Fourier變換有關(guān)

Tc = 50;

N = 200;

L = 1;

NN = 1:N;

NNinv = 1.0./NN;

fh = @(t) Tc/2 + ...

???Tc/pi*( ...

????(NNinv.*(1-cos(NN*pi))) * ( sin(NN*pi*t/L) )' ...

???) ;

t=0:0.1:5;

n = length(t);

y = zeros(1,n);

for i=1:n

???y(i) = fh(t(i));

end

plot(t, y, 'r*');

xl = 0;

xr = 1;

tbeg = 0;

dt = 0.0001

tend = 1e-2;

t = tbeg:dt:tend;

nt = length(t);

dx = 0.01;

x = xl:dx:xr;

nx = length(x);

T = zeros(nt,nx);

T(1,:) = 37; % 初始溫度是37攝氏度

i=1;

T(i,1) = fh(t(i)); %左邊界的溫度

T(i,end) = 37; %右邊界的溫度

K = 0.001*dt/dx/dx; % 熱傳導(dǎo)系數(shù)設(shè)置為 0.001

for i=2:nt

???T(i,2:end-1) = T(i-1,2:end-1) + ...

??????K * (?...

???????????????T(i-1,1:end-2) ...

????????????-2*T(i-1,2:end-1) ...

????????????+T(i-1,3:end) ...

?????????);?

T(i,1) = fh(t(i)); %左邊界的溫度

T(i,end) = 37; %右邊界的溫度

end

load_plugin("time")

pause(2)

clf; close all;

hold on

for i=1:nt

???plot(x(1:20),T(i,1:20))

???pause(0.1)

end

hold off