SciPy:solve_bvp问题二阶差异.情商 [英] SciPy: solve_bvp Problem 2nd Order Diff. Eq

问题描述

试图解决二阶差异.eq.有两个边界条件，但我尝试的任何方法似乎都没有用，而且我找不到包含与表达式中所用术语完全/相似的术语的教程，至少对我而言，scipy文档并未真正说明如何使用清楚地是resolve_bvp.

我有一个等式:y''+ 2/r * y'= A * y + b * y ^ 3其中y是r的函数.

我将其重写如下:

y1 = y(r)

y2 = y1'

所以y2'= -2/r * y2 + y1(A + b * y1 ^ 2)

在y'(0)= 0的情况下，y(r = 10)=常数

A，b和常数是已知的.

我有以下代码，但是它似乎不起作用，并且如上所述，我对文档有些困惑，因此不胜感激！

  def fun(x，y):返回np.vstack((y [2]，-2/x * y [1] + y [0] *(A + b * y [0] * y [0])))def bc(ya，yb):返回np.array([ya [2]，yb [1] -constant])x = np.linspace(1，10，10)ya = np.zeros((3，x.size))yb = np.zeros((3，x.size))sol_1 = resolve_bvp(有趣，bc，x，ya)sol_2 = resolve_bvp(有趣，bc，x，yb) 

谢谢！

==========编辑=========================我已经计算出一个解析解，只是看我是否还能从数值上找到相同的解，我认为主要问题是解有两个单独的区域，其中r≤r.R(输入)，其中r> R(输出).这导致两个不同的解决方案(仅在各自的域中有效)，条件为y_in(R)== y_out(R)和y_in'(R)== y_out'(R).

Trying to solve a 2nd order diff. eq. with 2 boundary conditions and nothing I try seems to work and I can't find a tutorial which includes all/similar terms to what I have in my expression and, at least to me, the scipy documentation doesn't really explain how to use solve_bvp clearly.

I have the equation: y'' + 2/r * y' = A * y + b * y^3 where y is a function of r.

I've rewritten it as follows:

y1 = y(r)

y2 = y1'

so that y2' = -2/r * y2 + y1(A + b * y1^2)

With conditions that y'(0) = 0, y(r=10) = constant

A, b and the constant are known.

I have the following code but it doesn't seem to work, and as said, I'm a little confused by the documentation so any help would be appreciated!

def fun(x, y):
    return np.vstack((y[2], -2/x*y[1]+y[0]*(A+b*y[0]*y[0])))

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


x = np.linspace(1, 10, 10)
ya = np.zeros((3, x.size))
yb = np.zeros((3, x.size))

sol_1 = solve_bvp(fun, bc, x, ya)
sol_2 = solve_bvp(fun, bc, x, yb)

Thanks!

==========EDIT========================= There is an analytical solution, for which I have calculated, it was just seeing if I could also find the same solution numerically, I think the major issue is that there are two separate regions for the solutions, one where r < R (in) and one where r > R (out). This leads to two different solutions (only valid in their respective domain) with the conditions that y_in(R) == y_out(R) and y_in'(R) == y_out'(R). Full 2-part solution where Radius = 1, a=99, b = 1 and constant = 1, y(inf) = constant

From Lutz Lehmann's solution, it gets the right shape (at least for the inner region, although not on the right scales).

I'm just unsure how you'd go about coding all the equating solutions and I suppose even getting their solutions in the first place, although Lutz's answer is an amazing point in the right direction. Thanks

解决方案

Issues

• The order of the equation is 2, thus the dimension of the state vector is 2, the value is always y[0], the derivative y[1], there is no y[2], probably a remnant from the translation from Matlab.

• Also in the boundary conditions, there is no ya[2], the derivative value is ya[1], and the function value in the second one is yb[0].

• The initial solution guess has to have the same number 2 of state components.

• Why do you compute two solutions with identical data?

• Remark: It is not necessary to convert the return values into numpy types, the solver checks and converts anyway.

BVP framework with singularity treatment

The ODE is singular at r=0, so one has to treat the first segment in a special way. The mean value theorem gives

(y'(r)-y'(0))/r->y''(0)  for  r->0,

so that in that limit r->0 you get

3*y''(0) = a*y(0) + b*y(0)^3`. 

This allows to define the first arc as

y(r) = y0 + (a*y0 + b*y0^3)*r^2/6
y'(r) = (a*y0 + b*y0^3)*r/3

up to order r³ and r². So if you want precision 1e-9 in y(r), then the first segment should not be longer than 1e-3.

Instead of trying to eliminate y0 from the equations for y(h) and y'(h) to get a condition connecting ya[0] and ya[1], let the solver do also this work and add y0 as parameter to the system. Then the boundary conditions have 3 slots corresponding to the added virtual dimension for the parameter, which can be naturally filled with the equations y(h)=ya[0] and ya[1]=y'(h) and the right boundary condition.

Altogether you can define the system as

h = 1e-3;

def fun(r, y, p):
    return  y[1], -2/r*y[1]+y[0]*(a+b*y[0]*y[0]) 

def bc(ya, yb, p):
    y0, = p
    yh = y0 + y0*(a+b*y0*y0)*h*h/6;
    dyh = y0*(a+b*y0*y0)*h/3
    return ya[0]-yh, ya[1]-dyh, yb[0]-c


x = np.linspace(h, 10, 10)
ya = np.zeros((2, x.size))

sol = solve_bvp(fun, bc, x, ya, p=[1])
print(sol.message,f"y(0)={sol.p[0]}");
plt.plot(sol.x, sol.y[0]);

so that with example parameters a, b, c = -1, 0.2, 3 you get a convergent solver call with y(0)=2.236081087849196 and the resulting plot

这篇关于SciPy:solve_bvp问题二阶差异.情商的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！