Finite Difference implicit solver (Crank-Nicolson) compare with analytical solution

Question

Khairul Hasan 2019-10-5

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/483750-finite-difference-implicit-solver-crank-nicolson-compare-with-analytical-solution

编辑： darova 2021-3-26

why my results from analytical and imilicit method doesn't match?

clc
close all
clf
L = 100;
delt = 10;
delx = 10;
D = 10;
n = (L/delx)+1;
lamda = (D*delt)/delx^2;
k = (0:delt:500);
T= zeros(length(k),n);
% Initial & Boundary Condition
T(:,1) = 16;
T(:,n) = 11;
T(1,2:n-1) = 16;
% Finite Difference implicit solver (Crank-Nicolson)
% Setup tri-diagonal matrix entries for the internal nodes
a = zeros(1,n);
b = zeros(1,n);
c = zeros(1,n);
d = zeros(1,n);
a(1)= 0; b(1)= 1; c(1)= 0; d(1)= 16;
a(n)= 0; b(n)= 1; c(n)= 0; d(n)= 11;
for q = 2:1:n-1
    a(q)= -(lamda);
    b(q)= 2*(1+ lamda);
    c(q)= -(lamda);
    d(q)= (2*T(1,q)) + (2*((T(1,q-1)- 2*T(1,q)+ T(1,q+1))/delx^2)*delt);
end 
for i = 2: length (k)
    T(i,:)= Tridiag (a,b,c,d); % Call tri-diagonal function 
    for q = 2:1:n-1
        d(q)= (2*T(i,q)) + (2*((T(i,q-1)- 2*T(i,q)+ T(i,q+1))/delx^2)*delt);% this allow to move to next time step
    end
end
% Analytical Solution
T1 = 16;
T2 = 11;
D = 10;
t=[10 20 40 80 160 250 500];
x=[0:10:100];
T_ana= zeros(length(x),length(t));
for i= 1:length(x)
    for j = 1 : length (t)
        constant = 0;
        for q = 1:50
            term1 = ((T2 - T1)*cos(q * pi))/q;
            term2 = sin((q*pi*x(i))/L);
            term3 = exp(-(D*(q)^2*(pi)^2*t(j))/(L)^2);
            constant = constant + (term1*term2*term3);
        end
        tfinal = T1+((T2 - T1)*x(i))/L+((2/pi)*(constant));
        T_ana(i,j)= tfinal;
    end
end