Пространственно-временные параметры в pde - решение с использованием метода линий в Python - PullRequest
Новая
       20

Пространственно-временные параметры в pde - решение с использованием метода линий в Python

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

Я хочу адаптировать этот метод решения на основе линий pde так, чтобы k было функцией как пространства, так и времени, ИЛИ равно нулю, если пороговый критерий не выполняется. Это не фактический pde, который я решаю, но это достаточно хороший пример.

Я могу достаточно легко сделать k-функцию t (см. Ниже MWE), однако я не могу понять, как сделать из нее функциюпространства и времени со сложностью порогового критерия, который отправляет один член RHS моей оды в ноль в любых узлах, в которых порог не превышен.

В качестве иллюстрации рассмотрим 

kSpace = X
kTime = 0.001*t if t < 2 and k = 0 otherwise.

if kSpace*kTime > 0.002:
    return k = kSpace*kTime
else:
    return 0

MWE с неправильной функцией k 

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

plt.interactive(False) # so I can see plots

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

def k(t):
    if t< 2:
        return 0.001*t
    else:
        return 0.1


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(t)*(u[i+1] - 2*u[i] + u[i-1]) / h**2 - 1.0
    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 голосов
/ 01 октября 2019

Я бы сделал это так, основываясь на вашем существующем коде. 

def k(x, t):
    return value_that_depends_on_x_and_t

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(X[i], t)*(u[i+1] - 2*u[i] + u[i-1]) / h**2 - 1.0
    return dudt
...