CODEWITHSUNDEEP
pythonscipydifferential-equationsboundary
Determinant
Determinant

Reputation: 11

scipy bvp solver for four differential equation systems

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

Upvotes: 1

Views: 809

Answers (3)

atronoush
atronoush

Reputation: 31

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

Upvotes: 1

Lutz Lehmann
Lutz Lehmann

Reputation: 26040

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.

enter image description here

Upvotes: 0

Futurologist
Futurologist

Reputation: 1914

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.

Upvotes: 1

Related Questions