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

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了网格中的填充点 - 前向欧拉算法 - 输出错误相关的知识,希望对你有一定的参考价值。

我将简要地试着向那些对数学经验不足的人解释我在做什么,这真的很简单。

我们正在尝试填充网格,如下所示:

enter image description here

我们找到橙色点,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))

我们可以很容易地将它写成线性方程组,如下所示:

enter image description here

现在我们只需要重复这个乘以这个矩阵(迭代时间变量)的过程。这是一种数值逼近偏微分方程解的简单方法。

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

这是代码

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,我如何制作网格,我似乎都会收敛到不正确的东西。他们并不亲密。

我不能为我的生活找出问题所在。有谁知道什么可能是错的?

答案

当你将间隔分成Delta_x = 1.0/N细胞时,你应该有N

N+1u[0],您可以在网格上获得u[N]点,但根据边界条件u[N]=u[0],您还只使用长度为N的数组来保存所有节点值。

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

Delta_t = Lambda * 2 * Delta_x

或者你的目标是O(dt, dx²)的方法的错误,以便有dt = c*dx^2有意义,但没有像c=1e-20这样的荒谬因素,如果你想时间离散化误差小对空间离散化误差,c=0.1c=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)

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


t离散问题的x精确解是cos(2*pi*(x+c*t)),其中c=sin(2*pi*dx)/(2*pi*dx)。如果您与该公式进行比较,则错误的大小应该非常小O(dt)

以上是关于网格中的填充点 - 前向欧拉算法 - 输出错误的主要内容,如果未能解决你的问题,请参考以下文章

Pandas:使用日期时间索引进行分组前向填充

我如何使用python将终端输出到网格中的tkinter框架?

地理网格搜索算法

当窗口/分区使用前向填充时,向 pyspark sql 中的 last() 函数添加条件

HDU - 6311 Cover (欧拉路径)

前向填充多列可重用功能代码