Finite Differential Method for the Heat Equation, CN Method

앞서 다뤘던 conditionally stable 한 FTCS Method 와는 다르게,

2019/02/26 - [Programming/MATLAB] - Finite Difference Method for the Heat Equation, FTCS Method

CN Method, Crank - Nicolson Method 는 Heat Equation을 위한 Finite Differential Method 이면서 stable 하다.

CN Method 는,

j 와 j+1 의 시간변화와, 공간에서 중앙의 2차 미분값 차이의 평균을 계산한다.

식으로 표현하자면,

가 된다.

, 로 치환을 하고 다시 정리를 해보자.

왼쪽은 시간이 j+1 일 때, 오른쪽은 시간이 j 일 때로 보기좋게 나누어 졌다.

이제,     라 하고, CN Method 를 Matrix Notation 으로 바꾸어보면,

이 때,

의 approximation 은  을 이용해 계산해준다. 이 때,

따라서 , 는  이다.

이를 이용해서 MATLAB 으로 CN Method 를 구현해보자.

행렬 U를 Output 으로 하는 HeatCN function 을 만들자.

U = heatCN(alpha,L,G0,GL,F,h,k,tf)

Output:

U : N+2, M+1의 온도를 store 하는 matrix

Input:

alpha: 열 전도율

L: 길이

G0,GL: 시간을 변수로 하는 Boundary Condition 의 벡터화 된 함수

F: 길이를 변수로 하는 초기 온도 값을 output으로 하는 벡터화 된 함수 

h: mesh의 공간 dimension을 지정

k: mesh의 시간 dimension을 지정

tf: final time value, 종료 시간 값

MATLAB:

function U = heatCN(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);
    
    % 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 CN method
    p = alpha * k / 2 / h^2;
    
    A = eye(n)*(1+2*p)+diag(ones(n-1,1)*(-p),1)+diag(ones(n-1,1)*(-p),-1);
    B = eye(n)*(1-2*p)+diag(ones(n-1,1)*(p),1)+diag(ones(n-1,1)*(p),-1);
    gj = zeros(n,1);
    
    for j = 2:m
        gj(1) = p * (G0((j-1)*k) + G0(j*k));
        gj(end) = p * (GL((j-1)*k) + GL(j*k));
        U(:,j) = A\(B*U(:,j-1) + gj);
    end

end

아래 스크립트를 실행해보면

alpha = 0.01; L = 1; tf = 12; h = 0.02; k = 0.01;
F = @(x) 125*sin(pi*x) + 50*sin(3*pi*x);
G0 = @(t) zeros(size(t));
GL = @(t) zeros(size(t));
U = heatCN(alpha,L,G0,GL,F,h,k,tf);
image(U,CDataMapping,scaled);
colormap jet; colorbar; axis off;

다음과 같은 figure가 나온다