我希望调整this method of lines based solution of a pde,以便k是空间和时间的函数,或者如果不满足阈值条件,k等于零。这不是我正在解决的实际偏微分方程,但它是一个足够好的例子
我可以很容易地使k成为t的函数(见下面的MWE),但是我不知道如何使它同时成为空间和时间的函数,因为阈值标准的复杂性使得我的ode的一个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
k函数不正确的MWE
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()
我会这样做的基础上,你现有的代码
相关问题 更多 >
编程相关推荐