проверка сходимости с помощью интегратора python - PullRequest
Новая
       35

проверка сходимости с помощью интегратора python

0 голосов
/ 09 января 2020

Я пытаюсь интегрировать разницу между моим численным и точным решением уравнения теплопроводности, хотя я не уверен, что было бы наилучшим способом решения этой проблемы. Существует ли специальный c интегратор, который позволил бы мне это сделать? Я надеюсь интегрировать его в $ x $.

Пока у меня есть следующий код: 

import numpy as np 
from matplotlib import pyplot 
from mpl_toolkits.mplot3d import Axes3D
from scipy import linalg
import matplotlib.pyplot as plt
import math


def initial_data(x):
    y = np.zeros_like(x)
    for i in range(len(x)):
        if (x[i] < 0.25):
            y[i] = 0.0
        elif (x[i] < 0.5):
            y[i] = 4.0 * (x[i] - 0.25)
        elif (x[i] < 0.75):
            y[i] = 4.0 * (0.75 - x[i])
        else:
            y[i] = 0.0

    return y

def heat_exact(x, t, kmax = 150):
    """Exact solution from separation of variables"""

    yexact = np.zeros_like(x)
    for k in range(1, kmax):
        d = -8.0* 
      (np.sin(k*np.pi/4.0)-2.0*np.sin(k*np.pi/2.0)+np.sin(3.0*k*np.pi/4.0))/((np.pi*k)**2)
         yexact += d*np.exp(-(k*np.pi)**2*t)*np.sin(k*np.pi*x)

     return yexact

def ftcs_heat(y, ynew, s):
    ynew[1:-1] = (1 - 2.0 * s) * y[1:-1] + s * (y[2:] + y[:-2])
    # Trivial boundary conditions
    ynew[0] = 0.0
    ynew[-1] = 0.0

Nx = 198
h = 1.0 / (Nx + 1.0)
t_end = 0.25
s = 1.0 / 3.0 # s = delta / h^2
delta = s * h**2
Nt = int(t_end / delta)+1

x = np.linspace(0.0, 1.0, Nx+2)
y = initial_data(x)
ynew = np.zeros_like(y)
for n in range(Nt):
    ftcs_heat(y, ynew, s)
    y = ynew


fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111)
x_exact = np.linspace(0.0, 1.0, 200)
ax.plot(x, y, 'kx', label = 'FTCS')
ax.plot(x_exact, heat_exact(x_exact, t_end), 'b-', label='Exact solution')
ax.legend()
plt.show()


diff = (y - heat_exact(x_exact,t_end))**2 # squared difference between numerical and exact solutions 

x1 = np.trapz(diff, x=x)  #(works!) 
import scipy.integrate as integrate

x1 = integrate.RK45(diff,diff[0],0,t_end) #(preferred but does not work)

Я хочу интегрировать переменную diff (квадратная разница). Любые предложения приветствуются, спасибо.

Редактировать: Я хотел бы использовать метод RK45, однако я не уверен, каким должен быть мой y0, я попытался x1 = integrate.RK45(diff,diff[0],0,t_end), но получаю следующую ошибку вывода: 

  raise ValueError("`y0` must be 1-dimensional.")

  ValueError: `y0` must be 1-dimensional.

1 Ответ

0 голосов
/ 09 января 2020

Под интеграцией вы подразумеваете, что хотите найти область между y и heat_exact? Или вы хотите знать, совпадают ли они в пределах указанного предела c? Последний может быть найден с numpy.isclose. В первом случае вы можете использовать несколько встроенных функций интеграции numpy.

Например: 

np.trapz(diff, x=x)

О, разве последняя строка не должна быть diff = (y - heat_exact(x_exact,t_end))**2? Моя интеграция этого diff дала 8.32E-12, что выглядит правильно, судя по графикам, которые вы мне дали.

Проверьте также scipy.integrate

...