Заполнение точек в сетке - алгоритм Форвард-Эйлера - неверный вывод - PullRequest
Новая
       120

Заполнение точек в сетке - алгоритм Форвард-Эйлера - неверный вывод

0 голосов
/ 20 октября 2018

Я очень кратко попытаюсь объяснить, что я делаю, тем, кто менее опытен в математике, это действительно довольно просто.

Мы пытаемся заполнить сетку следующим образом:

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])

Теоретически, я должен быть в состоянии найти лямбду достаточно маленькой, чтобы v_final и известное решение, u, были очень похожи.

Но я не могу.Независимо от того, насколько маленьким я делаю лямбду, каким точным я делаю сетку, мне кажется, что что-то не так.Они не близко.

Я не могу на всю жизнь понять проблему.У кого-нибудь есть идеи, что может быть не так?

1 Ответ

0 голосов
/ 20 октября 2018

Вы должны иметь Delta_x = 1.0/N, так как вы делите интервал на N клетки.

Вы получаете N+1 точек на сетке от u[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.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)

Метод Эйлера также не самый точный метод, ожидаемая ошибка находится в диапазоне exp(L*t_f)*dx^2 = e^5/N_x^2=0.58 для N_x=16, где L=1 былопринимается за приблизительную константу Липшица.Теперь, если вы увеличите до N_x=50, эта оценка ошибки уменьшится до 0.06, что также видно в результатах.

Точное решение x дискретной задачи *1036* равно cos(2*pi*(x+c*t))где c=sin(2*pi*dx)/(2*pi*dx).Если сравнивать с этой формулой, ошибки должны быть очень маленькими: O(dt).

...