Matlab代码的Python实现-有限差分法

编程问答 2022-06-03

问题描述

给出我老师创建的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,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,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,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代表实际的解决方案，红线是函数：

Actual solution compared to Finite Difference Method

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

Exact solution compared to FDM Matlab

翻译此代码时，您可能会看到任何明显的错误吗？

万一有人需要确切的解决方案：

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)

解决方法

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

Uoldold = Uold
Uold = Unext

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

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

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

Uoldold = Uold.copy()
Uold = Unext.copy()
matlabnumerical-methodsnumpypython