Matlab有限差分法的Python实现

Question

鉴于我的老师创建的Matlab代码：

function [] = explicitWave(T,L,N,J)
% Explicit method for the wave eq.
% T: Length time-interval
% L: Length x-interval
% N: Number of time-intervals
% J: Number of x-intervals
k=T/N;
h=L/J;
r=(k*k)/(h*h);
k/h

x=linspace(0,L,J+1); % number of points = number of intervals + 1

uOldOld=f(x); % solution two time-steps backwards. Initial condition
disp(uOldOld)
uOld=zeros(1,length(x)); % solution at previuos time-step
uNext=zeros(1,length(x));

% First time-step
for j=2:J
    uOld(j)=(1-r)*f(x(j))+r/2*(f(x(j+1))+f(x(j-1)))+k*g(x(j));
end

% Remaining time-steps
for n=0:N-1
    for j=2:J
        uNext(j)=2*(1-r)*uOld(j)+r*(uOld(j+1)+uOld(j-1))-uOldOld(j);
    end
    uOldOld=uOld;
    uOld=uNext;
end
 plot(x,uNext,'r')
end

我试图通过使用以下代码在Python中实现这一点：

import numpy as np
import matplotlib.pyplot as plt

def explicit_wave(f, g, T, L, N, J):
    """

    :param T: Length of Time Interval
    :param L: Length of X-interval
    :param N: Number of time intervals
    :param J: Number of X-intervals
    :return:
    """
    k = T/N
    h = L/J
    r = (k**2) / (h**2)

    x = np.linspace(0, L, J+1)
    Uoldold = f(x)
    Uold = np.zeros(len(x))
    Unext = np.zeros(len(x))

    for j in range(1, J):
        Uold[j] = (1-r)*f(x[j]) + (r/2)*(f(x[j+1]) + f(x[j-1])) + k*g(x[j])

    for n in range(N-1):
        for j in range(1, J):
            Unext[j] = 2*(1-r) * Uold[j]+r*(Uold[j+1]+Uold[j-1]) - Uoldold[j]

        Uoldold = Uold
        Uold = Unext
    
    
    plt.plot(x, Unext)
    plt.show()

    return Unext, x

但是，当我使用相同的输入运行代码时，在打印它们时会得到不同的结果。我的意见：

g = lambda x: -np.sin(2*np.pi*x)
f = lambda x: 2*np.sin(np.pi*x)
T = 8.0
L = 1.0
J = 60
N = 480

Python绘图结果与精确结果进行比较。x-es表示实际解决方案，红线表示函数：

Matlab绘图结果，x-es表示精确解，红线表示函数：

你能看到我在翻译这段代码时可能犯的任何明显错误吗

如果有人需要确切的解决方案：

exact = lambda x,t: 2*np.sin(np.pi*x)*np.cos(np.pi*t) - (1/(2*np.pi))*np.sin(2*np.pi*x)*np.sin(2*np.pi*t)
 

Answer 1

我通过调试发现了错误。这里的主要问题是代码：

Uoldold = Uold
Uold = Unext

因此，在Python中，当您将一个新变量定义为与旧变量相等时，它们将成为彼此的引用（即相互依赖）。让我以列表为例说明这一点：

a = [1,2,3,4]
b = a
b[1] = 10
print(a)
>> [1, 10, 3, 4]

所以这里的解决方案是使用.copy() 因此：

Uoldold = Uold.copy()
Uold = Unext.copy()