Finite Difference Method for the Heat Equation, FTCS Method

U 행렬을 Output으로 하는 function heatFTCS 를 구현하자

Euler's Method, Midpoint Method 와 같은 유한차분법을 Heat Equation 에도 적용이 가능하다.

이 MATLAB Function 은 Forward Time Central Space (FTCS) Finite Method 인데 원리는 다음과 같다.

위 행렬 U 에서 는 i (distance) metre 에서 j (time) seconds 일 때의 Temperature 를 나타낸다.

FTCS Method 는 다음과 같다.

한마디로,

같은 t 의 인접한 u 3세트를 가지고 다음 한 개의 항의 Temperature 를 유도하는 방식이다.

% alpha: thermal conductivity, the coefficient α in the heat equation
% L:length L
% G0 and GL: (vectorized) function handles with 1 parameter defining 
%            the time-dependent boundary conditions G0(t) and GL(t)
% F: (vectorized) function handle with 1 parameter defining the 
%    initial temperature distribution F(x)
% h: mesh size in the spatial dimension
% k: mesh size in the time dimension
% tf: final time value

function U = heatFTCS(alpha,L,G0,GL,F,h,k,tf)
    N = round(L/h)-1;
    M = round(tf/k);
    U = zeros(N+2,M+1);
    [n,m] = size(U);
    r = alpha*k/h^2;
    
    % Put x in first column
    for i = 1:n-1
        U(i+1,1) = i * h; 
    end
    % Put initial temperature distribution by F(x) on first column
    U(:,1) = F(U(:,1));

    % put t in first and last rows
    for j = 1:m-1
        U(1,j+1) = j * k;
        U(end,j+1) = j * k;
    end
    
    %Put time-dependent boundary condition on first and last rows
    U(1,:) = G0(U(1,:));
    U(end,:) = GL(U(end,:));
    
    % Apply FTCS Method
    for j = 2:m
        for i = 2:n-1
            U(i,j) = r*U(i+1,j-1) + (1-2*r)*U(i,j-1) + r*U(i-1,j-1);
        end
    end
end

Input Value:

alpha: Heat Equation 에서 Thermal conductivity의 값

L: Heat Equation 에서의 길이, Length

G0, GL: t 를 variable로 하는 x=0, x=L 에서의 Boundary Condition

F: x 를 variable 로 하는 t=0 에서의 Boundary Condition

h: U의 Mesh Size에서 Rows 를 결정짓는 값. 공간 (x) 를 표현한다

k: U의 Mesh Size에서 Columns 를 결정짓는 값. 시간 (t) 를 표현한다

tf: 최종시간 값, final t (seconds)

G0, GL, F 는 함수로 표현되며 모두 벡터화(Vectorized) 되었다.

*벡터화(Vectorized) 란?

만약, 

F = @(x) 0; 

처럼 정의된 함수 F 가 있을 때 input 이 어느 size 이어도 관계없이 0 을 return 한다.

이 때 이 함수는 vectorized 되지 않았다.

하지만

F = @(x) zeros(size(x));

처럼 정의되었다면, input의 size에 따라서 return 값의 size도 같이 달라질 것이다.

이런 함수를 Vectorized 되었다고 한다.

Matrix U 를 만들기 위해서는 L/h 와 tf/k 가 정수이어야 한다.

h와 k는 각각 L과 tf의 stepsize 이고, L/h와 tf/k 를 roundup 해준다.

N = round(L/h)-1;

M = round(tf/k);

그 후, N+2, M+1 사이즈의 영행렬을 만들어준다.

U = zeros(N+2,M+1);

U의 첫번째 열에 0부터 L까지 h의 stepsize를 가진 값들을 채워넣는다 (distance)

for i = 1:n-1 U(i+1,1) = i * h; end

그 후, F 함수를 이용해 첫번째 열에 Boundary Condition 을 apply 시킨다.

U(:,1) = F(U(:,1));

U의 첫번째 행과 마지막 행에 0부터 tf까지 k의 stepsize를 가진 값들을 채워넣는다 (time)

for j = 1:m-1 U(1,j+1) = j * k;

U(end,j+1) = j * k;

end

그 후, G0, GL 함수를 이용해 첫번째와 마지막 행에 Boundary Condition 을 apply 시킨다.

U(1,:) = G0(U(1,:)); U(end,:) = GL(U(end,:));

마지막으로 FTCS Method 를 Apply 시켜준다. 

for j = 2:m

for i = 2:n-1

U(i,j) = r*U(i+1,j-1) + (1-2*r)*U(i,j-1) + r*U(i-1,j-1);

end

end

*여기서 주의할점, FTCS 는 U행렬에서 previous column의 값으로 계산되는 method 이므로,

 column 을 먼저 계산해야지, row 를 먼저 계산하게 되면 틀린 값을 return한다.

다음 script를 실행했을 때의 결과,

alpha = 0.1; L = 1; tf = 50;
F = @(x) zeros(size(x)); G0 = @(t) sin(1.2*t); GL = @(t) -2*sin(0.5*t);
h = 0.01; k = 0.0005;
U = heatFTCS(alpha,L,G0,GL,F,h,k,tf);
image(U,'CDataMapping','scaled');
colormap jet; colorbar; axis off;

result: