I have a system of 3 differential equations (will be obvious from the code I believe) with 3 boundary conditions. I managed to solve it in MATLAB with a loop to change the initial guess bit by bit without terminating the program if it is about to return an error. However, on scipy
's solve_bvp
, I always get some answer, although it is wrong. So I kept changing my guesses (which kept changing the answer) and am giving pretty close numbers to what I have from the actual solution and it's still not working. Is there some other problem with the code perhaps, due to which it's not working? I just edited their documentation's code.
import numpy as np
def fun(x, y):
return np.vstack((3.769911184e12*np.exp(-19846/y[1])*(1-y[0]), 0.2056315191*(y[2]-y[1])+6.511664773e14*np.exp(-19846/y[1])*(1-y[0]), 1.696460033*(y[2]-y[1])))
def bc(ya, yb):
return np.array([ya[0], ya[1]-673, yb[2]-200])
x = np.linspace(0, 1, 5)
#y = np.ones((3, x.size))
y = np.array([[1, 1, 1, 1, 1], [670, 670, 670, 670, 670], [670, 670, 670, 670, 670] ])
from scipy.integrate import solve_bvp
sol = solve_bvp(fun, bc, x, y)
The actual solution is given below in the figure.
MATLAB Solution to the BVP