在 MATLAB 上使用 Crank-Nicolson 方法的热方程

Question

我正在尝试在这个方程的 matlab 中实现曲柄尼科尔森方法：

du/dt-d²u/dx²=f(x,t)
u(0,t)=u(L,t)=0
u(x,0)=u0(x)

与 :

- f(x,t)=20*exp(-50(x-1/2)²) if t<1/2; elso f(x,t)=0
- (x,t) belong to [0,L] x R+

边界条件为：

- U0(x)=0
- L = 1
- T = 1

这是我的数学思维：

of the form A*Un+1=B*Un+ht/2*Fn

我的问题是我得到的结果不一致而且我找不到我的错误。我的图表由正峰和负峰组成，这与等式

完全错误

这是我的代码：

%Parameters discretization in space according to the variable x
Nx=50; %Number of intervals, we will have Nx+1 discretization points
L=1; %Wire length
hx = 1/Nx ; %  Step of discretization in space
Xx = hx*[0:Nx]' ;    % Vector of discretized space

%Parameters discretization in space according to the variable t
Nt=200; %Number of intervals, we will have Ny+1 discretization points
T=1; %Time for which the heat propagation will be simulated
ht = 1/Nt ; %  Step of discretization in time
Xt = ht*[0:Nt]' ;    % Vector of discretized time

F=zeros(Nx+1,Nt+1); %Creation of the matrix that will contain the values of the function F(x,t)

for i=1:(Nt/2-1)
    F(:,i)=20*exp(-50*([0:hx:L]-1/2).^2); %Insertion in the matrix of the function F(x,t)=20*exp(-50(x-1/2)²) if t<1/2 and 0 otherwise
end

U=zeros(Nx+1,Nt+1); %Creation of the matrix that will contain the solutions of the equation

A=(1+2*alpha)*diag(ones(1,Nx+1))-alpha*diag(ones(1,Nx),1)-alpha*diag(ones(1,Nx),-1);
A(1,:)=0; %Zeroing of the first line, to enter the boundary conditions
A(end,:)=0; %Zeroing the last line to enter the boundary conditions
A(1,1)=1; %Boundary condition
A(end,end)=1; %Boundary condition

B=(1-2*alpha)*diag(ones(1,Nx+1))+alpha*diag(ones(1,Nx),1)+alpha*diag(ones(1,Nx),-1); %We write in the matrix B the terms which are repeated on the diagonals
B(1,:)=0; %Zeroing of the first line, to enter the boundary conditions
B(end,:)=0; %Zeroing the last line to enter the boundary conditions
B(1,1)=1; %Boundary condition
B(end,end)=1; %Boundary condition


for i = 1:(Nt)
    U(:,i+1)=((B*U(:,i)+(ht/2)*F(:,i)+(ht/2)*F(:,i+1))\A); %Iterative resolution of the system, we advance by one time step at each loop
end

surf(Xt,Xx,U)
xlabel("t");
ylabel("x");
title("Temperature distribution in the wire as a function of time")
shading interp

result

Answer 1

线性系统A*x=b的解在Matlab中表示为

x = A\b;

“分母位置b左边的A”，即A^(-1)*b（必要时使用pseudo-inverses）。

应用此更正，

for i = 1:(Nt)
    U(:,i+1)=A\(B*U(:,i)+(ht/2)*F(:,i)+(ht/2)*F(:,i+1)); 
end

解决方案的结果

在 MATLAB 上使用 Crank-Nicolson 方法的热方程

Heat equation with the Crank-Nicolson method on MATLAB

matlab

differential-equations