基于MATLAB的曲柄尼科尔森热方程

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

结果

Answer 1

用Matlab将线性系统A*x=b的解表示为

x = A\b;

"A在b的左分母位置“，即A^(-1)*b (必要时使用伪逆)。

运用这一修正，

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

解决方案的结果

​

​