网格中的填充点 - 正向欧拉算法 - 错误输出

Question

我会非常简短地尝试向那些没有数学经验的人解释我在做什么，这真的很简单。

我们正在尝试填充一个网格，如下：

我们找到橙色点，U(j,n+1)，使用它下面的三个连续点，U(j-1,n), U(j,n), U(j,n+1)

其中U在整个底行的值是给定的，并且是周期性的。所以理论上我们可以填满整个网格。

橙点计算公式为：

U(j,n+1) = U(j,n) + (delta_t / (2 * delta_x)) * (U(j+1,n) - U(j-1,n))

我们可以很容易地把它写成一个线性方程组如下：

现在我们只要重复这个乘以这个矩阵的过程（遍历时间变量）就可以了。这是一种用数值方法逼近偏微分方程解的简单方法。

我编写了执行此操作的代码，然后将我的最后一行与微分方程的已知解进行比较。

这是代码

import math
import numpy

def f(x):
    return math.cos(2 * math.pi * x)


def solution(x, t):
    return math.cos(2 * math.pi * (x + t))


# setting everything up
N = 16

Lambda = 10 ** (-20)
Delta_x = 1/(N+1)
Delta_t = Lambda * Delta_x * Delta_x
t_f = 5
v_0 = numpy.zeros((N, 1))

# Filling first row, initial condition was given
for i in range(N):
    v_0[i, 0] = f(i * Delta_x)

# Create coefficient matrix
M = numpy.zeros((N, N))
for i in range(N):
    M[i, i - 1] = -Delta_t / (2 * Delta_x)
    M[i, i] = 1
    M[i, (i + 1) % N] = Delta_t / (2 * Delta_x)

# start iterating through time
v_i = v_0
for i in range(math.floor(t_f / Delta_t) - 1):
    v_i = numpy.dot(M, v_i)

v_final = v_i
if (Delta_t * math.ceil(t_f / Delta_t) != t_f): #we don't reach t_f exactly using Delta_t
    v_final = (1/2) * (v_i + numpy.dot(M, v_i))

u = numpy.zeros(v_final.shape)
for i in range(N):
    u[i, 0] = solution(i * Delta_x, t_f)

for x in range(v_final.shape[0]):
    print (v_final[x], u[x])

从理论上讲，我应该能够找到足够小的 lambda，使得 v_final 和已知解 u 非常相似。

但我不能。不管我把 lambda 做得多小，我把网格做得多好，我似乎都会收敛到一些不正确的地方。他们不近。

我这辈子都弄不明白这个问题。有谁知道哪里出了问题？

Answer 1

您应该 Delta_x = 1.0/N，因为您将区间划分为 N 个单元格。

你在从 u[0] 到 u[N] 的网格上得到 N+1 个点，但是根据边界条件 u[N]=u[0]，你也只使用长度为 [= =13=] 保存所有节点值。

根据你给定的公式，你有 gamma = dt/(2*dx)，因此反向计算应该是 dt = gamma*2*dx 或者在你的变量名中

Delta_t = Lambda * 2 * Delta_x

或者您针对的是 O(dt, dx²) 方法的错误，因此 dt = c*dx^2 是有意义的，但没有像 c=1e-20 这样荒谬的因素，如果您希望时间离散化误差相对于 space 离散化误差较小，则 c=0.1 或 c=0.01 就足够了。

import numpy as np

def f(x):
    return np.cos(2 * np.pi * x)

def solution(x, t):
    return f(x + t)


# setting everything up
N_x = 16

Lambda = 1e-2
Delta_x = 1./N_x
Delta_t = Lambda * Delta_x * Delta_x
t_f = 5
N_t = int(t_f/Delta_t+0.5); t_f = N_t*Delta_t

# Filling first row, initial condition was given
x = np.arange(0,N_x,1) * Delta_x
v_0 = f(x)

# Create coefficient matrix
M = np.zeros((N_x, N_x))
for i in range(N_x):
    M[i, i - 1] = -Delta_t / (2 * Delta_x)
    M[i, i] = 1
    M[i, (i + 1) % N_x] = Delta_t / (2 * Delta_x)

# start iterating through time
v_i = v_0[:]
for i in range(N_t):
    v_i = np.dot(M, v_i)

v_final = v_i

u = solution(x, t_f)

for vx, ux in zip(v_final, u):
    print (vx, ux)

Euler 方法也不是最精确的方法，对于 N_x=16，预期误差在 exp(L*t_f)*dx^2 = e^5/N_x^2=0.58 范围内，其中 L=1 被视为近似 Lipschitz 常数。现在，如果您增加到 N_x=50，此误差估计会减少到 0.06，这在结果中也可见。

x 离散化问题的 t 精确解是 cos(2*pi*(x+c*t)) 其中 c=sin(2*pi*dx)/(2*pi*dx)。如果您与该公式进行比较，误差应该很小 O(dt).

网格中的填充点 - 正向欧拉算法 - 错误输出

Filling points in a grid - Forward Euler algorithm - wrong output

python

math

numerical-methods

differential-equations