计算科学 - 跳蛙法求解振弦波动方程误差过大 - 吾爱随笔录

跳蛙法求解振弦波动方程误差过大

计算科学 Python 计算物理学误差估计

2021-12-20 14:36:32

我一直在试图弄清楚过去两天我做错了什么。我不知道我是否真的做错了什么，或者在通常的越级问题中错误是否应该如此之大。我试着用一个精确的数值来检查我的数值解。具体来说，字符串的初始条件是和具有精确解 pde 是

y (x, 0) = x (1 - x)

$y(x,0) = x(1-x)$

\frac{\partial y}{\partial t} |_{t = 0} = 0

$\frac{\partial y}{\partial t} \vert_{t=0} = 0$

y (x, t) = \sum_{k = 1 k odd}^{\infty} \frac{8}{k^{3} π^{3}} \sin k π x \cos v k π t

$y(x,t) = \sum_{k=1\\ \textrm{k odd}}^{\infty} \frac{8}{k^3 \pi^3} \sin{}k \pi x \cos{v k \pi t}$

v^{2} y_{x x} = y_{t t}

$v^2 y_{xx} = y_{tt}$

第一个初始条件给我所有，第二个初始条件给我其中 C 是 Courant 数找到进一步的时间步 $y_{i,0}$ $i$

y_{i, 1} = y_{i, 0} + \frac{C^{2}}{2} (y_{i + 1, 0} + y_{i - 1, 0} - 2 y_{i, 0})

$y_{i,1} = y_{i,0} + \frac{C^2}{2}(y_{i+1, 0} + y_{i-1, 0} - 2y_{i,0})$

v Δ t / Δ x

$v\Delta t / \Delta x$

y_{i, j + 1} = 2 y_{i, j} - y_{i, j - 1} + C^{2} (y_{i + 1, j} + y_{i - 1, j} - 2 y_{i, j})

$y_{i,j+1} = 2y_{i,j} - y_{i,j-1} + C^2(y_{i+1,j} + y_{i-1,j} - 2y_{i,j})$

我确保满足 Courant 条件。波的幅度为 0.25 个单位，而仅在 316 个时间步之后就有 0.1 个单位的误差。大约是 0.316 秒。我不认为越级方法具有如此弱的可预测性（是吗？）。如果我增加时间步长，那么直到误差达到 0.1 的时间步数只会以相同的比例减少。我今天尝试将其与另一个精确的解决方案相匹配，但错误更大。我确信我的方法有问题，但我似乎无法弄清楚是什么。我非常有信心我的重复关系是正确的。如果误差真的应该这么大，那么对于这类问题的误差是否有任何分析公式可以用来正确选择我的参数？

这是我在 python 中的代码（使用更快的算法进行编辑，因为我只需要知道最后两个时间步的值即可计算下一个时间步的值）：

import numpy as np
from math import *

# Exact Soln
def y_exact(x,t):
    if t==0:
        return x*(1-x) 
    else:
        sum = 0
        for i in range(1,n):
            if i%2 == 1:
                sum += 8*sin(i*pi*x)*cos(v*i*pi*t)/((i**3)*(pi**3))
        return sum


n = 200
v = 1
tolerance = 0.1
L = 1
T = 1
x = np.linspace(0,L,10)
t = np.linspace(0,T,1000)
dx = x[1] - x[0]
dt = t[1] - t[0]
c = dx/dt
C = v/c
C2 = C*C
nx = len(x)
y = np.zeros((nx,3))
max_error = 0

# Initial condition 1
for i in range(nx):
    y[i,0] = y_exact(x[i],0)

# Initial condition 2
for i in range(nx):
    if i == 0 or i == (nx - 1) : # Boundary Condition
        y[i,1] = 0
    else:
        y[i,1] = y[i,0] + 0.5*C2*(y[i+1,0] + y[i-1,0] - 2*y[i,0])

for j in range(2,len(t)):
    if max_error >= tolerance :
        print(f"Error exceeded {tolerance} at time step number {j-1}.")
        break
    for i in range(nx):
        if i==0  or i == (nx -1):
            y[i,2] = 0
        else:
            y[i,2] = 2*y[i,1] - y[i,0] + C2*(y[i+1,1] + y[i-1,1] - 2*y[i,1])
            error = abs(y[i,2] - y_exact(x[i],dt*j))
            if error >= max_error :
                max_error = error
    for k in range(nx):
        y[i,0] = y[i,1]
        y[i,1] = y[i,2]

print(f"Maximum value of error was {max_error}")
```

1个回答

我使用另一种变体来精确解波动方程（与级数和相同），并通过将循环转换为 numpy 数组操作来缩短计算时间。使用此更改的代码，我没有观察到您得到的数字错误

v = 1
def f(x): 
    '''2-periodic odd-symmetric continuation of y(x,0)'''
    x = (x%2+1)%2-1; 
    return x*(1-abs(x))
def y_exact(x,t): return 0.5*(f(x+v*t)+f(x-v*t))

tolerance = 0.1
L = 1
T = 10
x = np.linspace(0,L,10)
t = np.linspace(0,T,1000)
dx = x[1] - x[0]
dt = t[1] - t[0]
c = dx/dt
C = v/c
C2 = C*C
nx = len(x)
y = np.zeros((nx,3))
max_error = 0

# Initial condition 1
y[:,0] = y_exact(x,0)

# Initial condition 2
y[[0,-1],1] = 0
y[1:-1,1] = y[1:-1,0] + 0.5*C2*(y[2:,0] + y[:-2,0] - 2*y[1:-1,0])

for j in range(2,len(t)):
    if max_error >= tolerance :
        print(f"Error exceeded {tolerance} at time step number {j-1}.")
        break
    y[[1,-1],2] = 0
    y[1:-1,2] = 2*y[1:-1,1] - y[1:-1,0] + C2*(y[2:,1] + y[:-2,1] - 2*y[1:-1,1])
    error = max(abs( y[:,2] - y_exact(x,dt*j) ))
    max_error = max(max_error,error)

    y[:,0] = y[:,1]
    y[:,1] = y[:,2]

print(f"Maximum value of error was {max_error}")

这运行到最后，最大报告误差为0.005687，远低于容差。人们必须更仔细地检查替换消除了编码错误的位置。

找到了。在最后一步，环绕，你有

    for k in range(nx):
        y[i,0] = y[i,1]
        y[i,1] = y[i,2]

迭代变量是k，使用的数组索引是i。

其它你可能感兴趣的问题

上一篇如何求解刚度矩阵？下一篇如何计算罗盘算子给出的无噪声图像的梯度？