scipy bvp solver for four differential equation systems

Question

I am trying to solve this systems but I get error.

I have to definition y3d=0 because y3'=0 in the equation systems. but when I did this, program cant solve. if I say y3d=y[3] then program run,

equation system that ı have to solve is like this:

dy1/dx=y2

dy2/dx=-y3*y1

dy3/dx=0

dy4/d=y1**2+y2**2 and boundary condition y1(0)=y1(1)=0 and y4(0)= 0 y4(1)=1

can scipy handle this?

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

def eqts(x,y):
    

    y1d=y[1]
    
    y2d=-y[2]*y[0]
    
    y3d=0
    
       
    y4d=y[0]**2+y[1]**2
    
    return np.vstack((y1d,y2d,y3d,y4d))


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

x = np.linspace(0,1,10)
y= np.zeros((4,x.size))
y[2,:]=1



sol=solve_bvp(eqts,bc,x,y)

Unfortunately I get the following error message ;


ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 1, the array at index 0 has size 10 and the array at index 2 has size 1

The BC return vector should be `[ya[0],yb[0],ya[3],yb[3]-1]` per the given BC in the given formulas. — Lutz Lehmann, Dec 29 '21 at 21:07
yes you are right. I adjusted boundary condition but I get same error. How should I define the y3d? — Determinant, Dec 29 '21 at 22:14
The initial `y` needs to contain a function table for every component of the state, thus start with `y=np.zeros([4,len(x)])`. Then the line after that should also make sense. You could put the eigenvalue `y3` into the parameter that is called `p` in the `solve_bvp` interface. — Lutz Lehmann, Dec 29 '21 at 22:51
thank you, y=np.zeros((4,x.size)). If I put the y3 into the p parameter, how can I find the different y3 value? actually I have to find the y3 with this method that use four differential equation system with y1,y2,y3 and y4. because I will need to use this method later for more complex problem. — Determinant, Dec 29 '21 at 23:13

score 1 · Answer 1 · answered Dec 29 '21 at 21:02

1

Well, first of all, in your script, your boundary conditions are overdetermined. Nowhere it is said that y3(0) = 0 or y3(1) = 0. Actually, it is not: y3(t) is a constant but it is not zero. If you impose such condition y3(t) = 0, things will not work at all. On top of that, this system looks non-linear (quadratic) but actually is a linear system. You can solve it explicitly without python. If I am not mistaken, the only way you can have a solution is when y3 > 0, which gives you

y1(t)  =  B * sin(k*pi*t)
y2(t)  =  k*pi*B * cos(k*pi*t)
y3(t)  =  k^2*pi^2
y4(t)  =  t  +  (k^2pi^2 - 1) B^2 * sin(2*k*pi*t) / (4*k*pi)

where B = sqrt( 2*pi*k / (k^2*pi^2 + 1) )
and k is an arbitrary non-zero integer

or at least something along those lines.

answered Dec 29 '21 at 21:02

Futurologist

1,874
2
7
9

1

The BC are contradictory as at the same time `yb[3]=0` and `yb[3]=1` is demanded, that is, two different values for `y4(1)`. Note the index shift, python arrays start with a zero index. – Lutz Lehmann Dec 29 '21 at 21:09
@LutzLehmann Yes, you are right, I guess I did not pay enough attention to the OP. – Futurologist Dec 29 '21 at 21:21
yes, I can solve analytical but I need to get y3 values and y function's graph that provide the equation with python. Because ı want to see if bvp solver can solve the differential system that I know the analytical solution. If ıt solve then I will use the solver at more complex problem. – Determinant Dec 29 '21 at 22:27

score 1 · Answer 2 · answered Mar 03 '22 at 11:55

1

Determinant

The main issue is here

    y1d=y[1]
    
    y2d=-y[2]*y[0]
    
    y3d=0
    
       
    y4d=y[0]**2+y[1]**2

In your implementation all y1d, y2d, and y4d are vector, but y3d is scalar! You may use y3d = np..zeros_like(y1d) solve_bvp requires to return a rhs with same size of y see https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html?highlight=solve_bvp#scipy.integrate.solve_bvp

answered Mar 03 '22 at 11:55

atronoush

21
2

One could say in more detail that `solve_bvp` calls its ODE function in multi-point evaluation or "vectorized" mode. Which means that the argument `t` is an numpy array of the time grid and `y` contains in each component the list of all the corresponding component values over the grid. Thus the origin of the vectors in `y[0]` etc. One could also shortly set `y3d = 0*y[2]`. – Lutz Lehmann Mar 03 '22 at 12:32

Lutz Lehmann · Answer 3 · 2021-12-30T17:27:04.703

In scipy's solve_bvp there is the possibility to treat constant components as parameters that are to be fitted. Thus one can define the system in dimension 3 as

def eqn(t,x,p): return x[1], -p[0]*x[0],x[0]**2+x[1]**2

def bc(x0,x1,p): return x0[0], x1[0], x0[2], x1[2]-1

t_init = np.linspace(0,1,6)
x_init = np.zeros([3,len(t_init)])
x_init[0,1]=1; # break the symmetry that gives a singular Jacobian

Another variant to get a general initial guess would be to fill it with random noise.

To get different solutions the initial conditions need to be different, one point of influence is the frequency square. Setting it close to the expected values gives indeed different solutions.

for w0 in np.arange(1,10)*3:
    res = solve_bvp(eqn,bc,t_init,x_init,p=[w0**2])
    print(res.message,res.p[0])

This gives as output

The algorithm converged to the desired accuracy. 9.869869348810965
The algorithm converged to the desired accuracy. 39.47947747757308
The algorithm converged to the desired accuracy. 88.82805487260174
The algorithm converged to the desired accuracy. 88.8280548726001
The algorithm converged to the desired accuracy. 157.91618155379464
The algorithm converged to the desired accuracy. 157.91618155379678
The algorithm converged to the desired accuracy. 355.31352482772894
The algorithm converged to the desired accuracy. 355.3135248277307
The algorithm converged to the desired accuracy. 157.91618163528014

As one can see, in the higher frequencies the given initial frequency gets balanced against the lower frequency behavior of the initial functions. This is a general problem, if not forced to stay orthogonal to the lower frequency solutions, the solver tends towards the smoother lower frequencies.

Adding plot commands plt.plot(res.x, res.y[0]) etc. shows the expected sinusoidal solutions.

thank you so much. is the bvpsolver's algotithm based on collocation medthod? — Determinant, Dec 30 '21 at 13:55
Yes, afaik. Combined with adaptive mesh refinement. But not the other way, no coarsening of the interval subdivision to increase the local error density to the desired error level if the solution changes drastically during the solver iterations, making a formerly "curved" segment "flat". — Lutz Lehmann, Dec 30 '21 at 14:04
thank you, now I plotted the y function for first p0 but I dont get sinus function. but p0 values looks right. — Determinant, Dec 30 '21 at 15:50
I see perfectly the sine function. Note that for the basis frequency you get only half a period. — Lutz Lehmann, Dec 30 '21 at 17:27
yes but the blue line is not sine function, I didnt understand, first eigenfunction's analytical soluiton is sin(pi*x) — Determinant, Dec 30 '21 at 17:43
Yes it is, it is `sin(pi*t)`, times a constant, over the interval `[0,1]`, half a period. — Lutz Lehmann, Dec 30 '21 at 17:45

scipy bvp solver for four differential equation systems

3 Answers3