scipy 在边值问题的定义中如何正确使用强迫项

slmsl1lt  于 2023-08-05  发布在  其他
关注(0)|答案(1)|浏览(120)

我给出了一个我试图实现的代码的简化示例

from scipy.integrate import solve_bvp
import numpy as np
p = [1, 2]
a = 0.5
num = 100
x = np.linspace(0, 1, num)

var = np.ones((1, int(num/2)))
for i in range(int(num/2)):
    var[0, i] = i**2

def fun(x, y):
    force = np.zeros_like(x)
    force[x <= 1-a] = -((np.polyval(p, y[0]))[x <= 1-a])*np.interp(x, x, var)

    return np.vstack((y[1], y[2], y[3], force))

def bc(ya, yb):
    return np.array([ya[0], ya[1], yb[0]-1, yb[2]])

y_guess = 0*np.ones((4, x.size))
res = solve_bvp(fun,bc,x,y_guess)

字符串
我面临以下错误

ValueError: object too deep for desired array


我还想澄清的是,我已经将var矩阵构造为x^2,并希望使用np.interp函数来用于我试图求解的实际四阶边值问题的强制项。我的强迫项是不连续的,由多值函数给出,x在(0,0.5)中,否则为零。删除interp函数可以使代码正常工作,但我无法找到包含此术语的正确方法。

tzcvj98z

tzcvj98z1#

为了首先解决错误,interp需要三个参数:xxpfpxpfp是你知道的x和y值,它们应该是一对一对应的1D数组。目前,var的形状为(1,50),这不是1D,因此出现了“太深”错误。您可以修复此问题并将其更改为var.ravel(),但随后会出现错误,即形状不匹配,因为x.shape = (100,)var.ravel().shape = (50,)
无论哪种方式,我都会将您的代码更改为具有单独的函数来处理f(y)var(x)。现在,f(y)实际上是不正确的,因为它实际上是x和y的函数,即f=f(x,y)(since f(x>a,y) = 0).为了简单起见,我也将var(x)称为g(x)vars是一个Python关键字,所以我避免使用var,因为它很接近)。您希望ODE调用这两个函数,因此可以将它们设为参数,并使用functools.partial处理只需要x和y的函数的solve_bvp
对于f(x,y),我使用np.polynomial.Polynomial创建多项式并调用它(注意:这类采用与np.polyval相反顺序的系数,因此颠倒)。然后我将x > a的点更新为0。
对于g(x),我只是使用了您提到的x**2表单,但您可以根据需要更改此表单。但是,因为g(x)只是x的函数,并且x是常数(您设置这些值),所以g(x)(或x的任何函数)是常数,实际上不需要函数。
下面是我将使用的solve_bvp代码,尽管sol.successFalse,这可能是这个问题被编造出来的结果,没有任何有效的解决方案。

import numpy as np
from scipy.integrate import solve_bvp
from functools import partial

p = [1., 2.]
poly = np.polynomial.Polynomial(p[::-1])
a = 0.5
num = 100
x = np.linspace(0, 1, num)

def g(x):
    return x**2

def f(x, y):
    res = poly(y)
    res[x > a] = 0.
    return res

def ode_full(x, y, f, g):
    yp1 = y[1]
    yp2 = y[2]
    yp3 = y[3]
    yp4 = -f(x, y[0])*g(x)
    
    return np.vstack([yp1, yp2, yp3, yp4])

def bc(ya, yb):
    return np.array([ya[0], ya[1], yb[0]-1, yb[2]])

y_guess = np.ones((4, x.size))

ode = partial(ode_full, f=f, g=g)
sol = solve_bvp(ode, bc, x, y_guess)
print(sol.success)

字符串

相关问题