# numerical solution for differential equation

I have these 3 equations and i want to solve them with numerical methods. so I am using scipy library but I don't know how to solve 3 equations together. R, g, sigma and density are constants.

\begin{align} &\frac{d\Phi}{ds} = -\frac{\sin\Phi}{x} + \frac{2}{R} + \frac{\Delta\rho gz}{\sigma}\, ,\\ &\frac{dx}{dx} = \cos\Phi\, ,\\ &\frac{dz}{ds} = \sin\Phi\, \end{align}

with $$\Phi(0) = x(0) = z(0) = 0$$.

From what I can tell this is an initial value ODE problem, so you should be able to use solve_ivp. This handles vector equations just fine. I'm not sure exactly how stiff your problem is, so using the defaults for most of the settings should be fine, and you can adjust these if you find the numerical solution isn't satisfactory.

There is one slightly problematic case you'll have to handle, and that's at $$s=0$$ you have to handle $$-\frac{\sin{\phi}}{x}$$ using l'Hopital's rule, and deal with anytime that x could get near 0 somehow later in the simulation.

Ok, so first, you have a problem: when $$s=0$$ both $$x(0) = 0$$ and $$\sin\Phi(0) = \sin(0) = 0$$. To figure out how to deal with this, try to extend the function $$\frac{\sin\Phi(s)}{x(s)}$$ continuously for $$s = 0$$. That means calculate the limit $$L = \lim_{s \to 0} \, \frac{\sin\Phi(s)}{x(s)}$$ The limit exists, because roughly speaking, the numerator is bounded. Notice, that by initial conditions, $$\Phi(0) = 0$$ and $$x(0) = 0$$ and are continously differentiable functions, so \begin{align} L = & \lim_{s \to 0} \, \frac{\sin\Phi(s)}{x(s)} = \lim_{s \to 0} \left(\, \frac{\Phi(s)}{x(s)}\, \, \frac{\sin\Phi(s)}{\Phi(s)}\right) \\ = & \left( \, \lim_{s \to 0} \, \frac{\Phi(s)}{x(s)}\, \right) \left( \, \lim_{s \to 0}\frac{\sin\Phi(s)}{\Phi(s)}\right) = \left( \, \lim_{s \to 0} \, \frac{\Phi(s)}{x(s)}\, \right) 1 \\ = & \lim_{s \to 0} \, \frac{\Phi(s)}{x(s)} \end{align} By l'Hospital's rule \begin{align} L = & \lim_{s \to 0} \, \frac{\Phi(s)}{x(s)} = \lim_{s \to 0} \, \frac{\, \frac{d\Phi}{ds}(s) \,}{\, \frac{dx}{ds}(s)\,} = \lim_{s \to 0} \, \frac{\, \frac{d\Phi}{ds}(s) \,}{\, \cos\Phi(s)\,}\\ = & \, \frac{\, \lim_{s \to 0} \frac{d\Phi}{ds}(s) \,}{\, \lim_{s \to 0} \cos\Phi(s)\,} = \frac{\, \lim_{s \to 0} \frac{d\Phi}{ds}(s) \,}{\, \cos\Phi(0)\,} = \frac{\, \lim_{s \to 0} \frac{d\Phi}{ds}(s) \,}{\, \cos(0)\,} = \lim_{s \to 0} \frac{d\Phi}{ds}(s)\\ =& \lim_{s \to 0} \left( - \frac{\sin\Phi(s)}{x(s)} + \frac{2}{R} + \frac{\Delta \rho g \, z(s)}{\sigma}\right)\\ =& \lim_{s \to 0} \left( - \frac{\sin\Phi(s)}{x(s)} \right) + \lim_{s \to 0} \left( \frac{2}{R} \right) + \lim_{s \to 0} \left( \frac{\Delta \rho g \, z(s)}{\sigma}\right) \\ =& \, - \, \lim_{s \to 0} \left( \frac{\sin\Phi(s)}{x(s)} \right) + \frac{2}{R} + \frac{\Delta \rho g \, z(0)}{\sigma} =\, - \, L + \frac{2}{R} + 0 \\ =& \, - \, L + \frac{2}{R} \end{align} By solving the latter equation for $$L$$, one gets $$L = \frac{1}{R} = \lim_{s \to 0} \, \frac{\sin\Phi(s)}{x(s)}$$

Based on this fact, I wrote some python code and hoped for the best, the best being that $$x(s)$$ doesn't become zero after the initial step. Of course, I do not know the constants, so I manufactured some random ones and assumed that $$\Delta$$ is also some kind of a constant and not the Laplacian applied to $$\rho$$.

import numpy as np
from scipy.integrate import solve_ivp

# y = np.array([Phi, x, z])
def F(t, y):
cos_Phi = np.cos(y)
sin_Phi = np.sin(y)
if abs(y) + abs(y) < 0.00000001:
dPhi = 1/R + k*y
else:
dPhi = - sin_Phi/y + 2/R + k*y
return np.array([dPhi, cos_Phi, sin_Phi])

R = 7
Delta = 1
rho = 0.7
g = 9.8
sigma = 2
k = Delta*rho*g/sigma

s_0 = 0
s_1 = 10
y0 = np.array([0, 0, 0])
t_span = np.linspace(s_0, s_1, num=1000)

sol = solve_ivp(F, [s_0, s_1], y0, method='Radau', t_eval=t_span)

fig, ax = plt.subplots(3,1)
ax.plot(sol.t, sol.y[0,:])
ax.plot(sol.t, sol.y[1,:])
ax.plot(sol.t, sol.y[2,:])
fig2, ax2 = plt.subplots(1,1)
ax2.set_aspect( 1 )
ax2.plot(sol.y[1,:], sol.y[2,:])
plt.show()


Here are the plots of $$\Phi, \, x, \, z$$ and this is the plot of the trajectory $$\big(x(s), \, z(s) \big)$$ 