Ошибочная реализация scipy.integrate.RK45 - PullRequest
Новая
       41

Ошибочная реализация scipy.integrate.RK45

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

дальше моего вопроса секунд а go, у меня есть еще один вопрос при попытке реализовать решатель scipy.integrate.RK45 для решения интеграла, который вычисляет разницу между моим численным и точным решением.

Код работает до последних 3 строк фактического кода. Я получаю следующее сообщение об ошибке: raise ValueError("`y0` must be 1-dimensional.") ValueError: `y0` must be 1-dimensional

Я не уверен, каким должен быть мой ввод для y0. Любые предложения приветствуются, вероятно, что-то простое.

Пожалуйста, найдите код ниже: 

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
x = np.trapz(diff, x=x)
import scipy.integrate as integrate
x = integrate.RK45(diff,diff[0],0,t_end)
...