Решающая система дес с использованием метода линий - пространственно дискретная - PullRequest
Новая
       14

Решающая система дес с использованием метода линий - пространственно дискретная

0 голосов
/ 30 сентября 2019

Я пытаюсь адаптировать этот пример, используя метод линий, поэтому решите pde для системы дес. Это не та система, которую я пытаюсь решить, просто пример.

Как мне включить вторую де? Я пытался изменить odefunc, чтобы он был 

def odefunc(u,t):
    dudt = np.zeros(X.shape)
    dvdt = np.zeros(X.shape)

    dudt[0] = 0 # constant at boundary condition
    dudt[-1] = 0

    dvdt[0] = 0
    dvdt[-1] = 0


    # for the internal nodes
    for i in range (1, N-1):
        dudt[i] = k* (u[i+1] - 2*u[i] + u[i-1]) / h**2
        dvdt[i] = u[i]
    return [dudt,dvdt]

на основе нескольких примеров решений систем од, которые я нашел - но, поскольку у меня уже был массив для dudt, я подозреваю, что это может быть моей проблемой. Я также не знаю, как должны выглядеть мои начальные условия, поэтому они имеют правильные размеры и т. Д.

Я получаю ошибку The array return by func must be one-dimensional, but got ndim=2. в строке sol = odeint(odefunc, init, tspan)

Пример содин pde 

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

plt.interactive(False) 

N = 5 # number of points to discretise
L = 1.0 # length of the rod
X = np.linspace(0,L,N)
h = L/ (N  - 1)

k = 0.02

def odefunc(u,t):
    dudt = np.zeros(X.shape)

    dudt[0] = 0 # constant at boundary condition
    dudt[-1] = 0


    # for the internal nodes
    for i in range (1, N-1):
        dudt[i] = k* (u[i+1] - 2*u[i] + u[i-1]) / h**2
    return dudt

init = 150.0 * np.ones(X.shape) # initial temperature
init[0] = 100.0 # boundary condition
init[-1] = 200.0 # boundary condition

tspan = np.linspace(0.0, 5.0, 100)
sol = odeint(odefunc, init, tspan)

for i in range(0, len(tspan), 5):
    plt.plot(X,sol[i], label = 't={0:1.2f}'.format(tspan[i]))

# legend outside the figure
plt.legend(loc='center left', bbox_to_anchor=(1,0.5))
plt.xlabel('X position')
plt.ylabel('Temperature')

# adjust figure edges so the legend is in the figure
plt.subplots_adjust(top=0.89, right = 0.77)
plt.show()

1 Ответ

0 голосов
/ 30 сентября 2019

Из документации odeint функции ode определяет y как вектор, но он определен как список векторов.

Решает проблему начальных значений для жестких или не жестких систем од первого порядка:

dy/dt = func(y, t, ...) [or func(t, y, ...)]

, где y может быть вектором.

Вы можете просто использовать функцию ravel, чтобы создать вектор из списка векторов. Вот пример из вашего кода: 

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
plt.interactive(False) 
N = 5 # number of points to discretise
L = 1.0 # length of the rod
X = np.linspace(0,L,N)
h = L/ (N  - 1)
k = 0.02
def odefunc2(y,t):

    u = y[:N]
    v = y[N:]

    dudt = np.zeros_like(u)
    dvdt = np.zeros_like(v)    

    dudt[0] = 0 # constant at boundary condition
    dudt[-1] = 0

    dvdt[0] = 0
    dvdt[-1] = 0
    # for the internal nodes
    for i in range (1, N-1):
        dudt[i] = k* (u[i+1] - 2*u[i] + u[i-1]) / h**2
        dvdt[i] = u[i]
    return np.ravel([dudt, dvdt])

initu = 150.0 * np.ones(X.shape) # initial temperature
initu[0] = 100.0 # boundary condition
initu[-1] = 200.0 # boundary condition

initv = np.zeros(X.shape) # initial v

init = np.ravel([initu, initv])

tspan = np.linspace(0.0, 5.0, 100)
sol = odeint(odefunc2, init, tspan)

fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2)

for i in range(0, len(tspan)):
    ax1.plot(X, sol[i, :N], label = 't={0:1.2f}'.format(tspan[i]))
    ax2.plot(X, sol[i, N:], label = 't={0:1.2f}'.format(tspan[i]))

plt.show()
...