May-08-2022, 10:23 PM
Hi!
I have a system of 10 nonlinear ODE with 10 variables (y[0], y[1], ... y[9]), there are 10 initial conditions for y[i] for x=0 and 4 conditions for x = 6.804.
I try to solve a boundary value problem using scipy.integrate.solve_bvp according to the following algorithm: https://docs.scipy.org/doc/scipy/referen...e_bvp.html
I get the error "_vhstack_dispatcher() takes 1 positional argument but 9 were given" though I tried to use np.dstack instead of np.hstack. Please help to fix the error.
My code is below:
I have a system of 10 nonlinear ODE with 10 variables (y[0], y[1], ... y[9]), there are 10 initial conditions for y[i] for x=0 and 4 conditions for x = 6.804.
I try to solve a boundary value problem using scipy.integrate.solve_bvp according to the following algorithm: https://docs.scipy.org/doc/scipy/referen...e_bvp.html
I get the error "_vhstack_dispatcher() takes 1 positional argument but 9 were given" though I tried to use np.dstack instead of np.hstack. Please help to fix the error.
My code is below:
import numpy as np from scipy.integrate import solve_bvp import matplotlib.pyplot as plt from scipy.constants import g def fun(x, y): return np.vstack((ptb / y[4] * y[5] / np.sqrt(y[5]**2 + y[6]**2) * function_deo(y[4], y[5], y[6], y[9], ptb, mqb, pes) - 1 / y[2]**2 + y[1]**2 / y[2], ptb / y[4] * y[6] / np.sqrt(y[5]**2 + y[6]**2) * function_deo(y[4], y[5], y[6], y[9], ptb, mqb, pes) - y[0] * y[1] / y[2]), y[0], y[1] / y[2], -mqb * function_deo(y[4], y[5], y[6], y[9], ptb, mqb, pes), y[6] * y[1] / y[2] - y[7], -y[5] * 2.0 * y[1] / y[2] + y[6] * y[0] / y[2] - y[8] * 1.0 / y[2], y[5] * (-2.0 / y[2]**3 + y[1]**2 / y[2]**2) + y[6] * (-y[0] * y[1]) / y[2] + y[8] * y[1] / y[2]**2, 0, ptb / y[4]**2 * (np.sqrt(y[5]**2 + y[6]**2)) * function_deo(y[4], y[5], y[6], y[9], ptb, mqb, pes)) #Implement evaluation of the boundary condition residuals: def bc(ya, yb): return np.array([ya[-0.03899746484, 1.05800280714, -0.30648856279, -0.605427, 1, -0.065851, 1.361805, 1.419320, -8.219418e-5, 0.6], yb[3.230021179060866e-17, 0.5519654178809713, 2.221001221001221, np.pi]]) #Define the initial mesh with 20 nodes: x = np.linspace(0.0, 6.804, 20) #To obtain solutions, different initial guesses for y are used. They are denoted by subscripts a and b. y_a = np.zeros((10, x.size)) y_b = np.zeros((10, x.size)) y_b[0] = 3.230021179060866e-17 res_a = solve_bvp(fun, bc, x, y_a) res_b = solve_bvp(fun, bc, x, y_b) x_plot = np.linspace(0, 6.804, 100) y_plot_a = res_a.sol(x_plot)[0] y_plot_b = res_b.sol(x_plot)[0] #plt.plot(x_plot, y_plot_a, label='y_a') #plt.plot(x_plot, y_plot_b, label='y_b') plt.legend() plt.xlabel("x") plt.ylabel("y") plt.show() 