# How to solve second order coupled non linear differential equations

For a project I am doing, I have to solve the following system of differential equations numerically using my own code:

$$x^2K'' = KH^2 + K(K^2-1)$$ and, $$x^2H'' = 2K^2H + \alpha H(H^2-x^2)$$ Here, $$K$$ and $$H$$ are the dependent variables, and $$\alpha$$ is a parameter which I have to vary from 0 to infinity. An analytical solution exists only for $$\alpha=0$$ case. The boundary conditions are, $$K(\epsilon) = 1 \\ H(\epsilon) = 0$$ and, $$K(\bar{\epsilon}) = 0 \\ H(\bar{\epsilon}) = \bar{\epsilon}$$ Here, $$\epsilon$$ is an arbitrary small number, (say 0.0001), and $$\bar{\epsilon}$$ is an arbitrarily large number, (say 400).

In short, boundary condition is that $$K$$ goes to 0 while $$H$$ goes to $$x$$.

What I have tried:

For the analytical case, I have tried RK4, RKF(Fehlberg), RK4 based Shooting, and scipy's solve_ivp. However, all suffer from the same problem, they diverge from analytical solution at around 10 (see attached graphs, first is of the two functions compared to analytical, and second is of derivatives).

I have also tried iterative finite difference approximation since I have initial guess, it suffers from same issue.

The only inbuilt function that works is solve_bvp. Can anyone suggest me two or three numerical methods for this system, for arbitrary $$\alpha$$, to solve as bvp and point at some good resources for those? What I have is an initial good guess for the solutions, which is the analytical solution for $$\alpha=0$$ case.

EDIT: This might be helpful. The solution for analytical case is,

$$K = \frac{x}{\sinh(x)}$$ and $$H = \frac{x}{\tanh(x)} - 1$$ Edit2: Also, I even tried introducing new variables, $$k$$ and $$h$$, which: Trial1: Take out a power of x, i.e. $$K=xk$$ and likewise for $$H$$. Trial2: Take out asymptotic behavior, i.e. $$K = e^{-x}k$$ and $$H=xh$$

I recast the equations, still the same result.

UPDATE: So I tried a couple more methods and implemented solve_bvp's residual control method as well. Based on suggestion of @superbee I implemented the Finite Difference Method with Newton Ralphson technique and it works!! Briefly, I generated a system of non-linear algebraic equations by giving an initial guess and finite difference approximation. I found the Jacobian, and found an updated guess as:

$$y^{(i+1)} = y^{(i)} - J^{-1} F$$

• Could you update the plots? I can barely see the numbers and labels. Apr 25 at 19:56
• I would suspect that you can analytically show that $K(t)>0, H(t)>0$ if you have $K(0)>0, H(0)>0$. But numerically, this may not be the case and then you get a runaway effect. You might want to enforce this constraint in a different way. Apr 26 at 0:23
• Separately, what happens if you decrease the step size? Apr 26 at 0:23
• @nicoguaro Updated, I hope it is clear now. Also, red and green are completely overlaid. Apr 26 at 4:02
• All the typical ODE integrators guarantee is that the solution converges as $\Delta t\to 0$. They do not guarantee that structural properties such as positivity of solutions are satisfied. In your case, you want that $H(0)=0$, which is a structural property. If that is important to you, you need to find ways to enforce this condition -- because without this, all you will get is that the solution converges to this value, but probably not satisfy it exactly. Apr 26 at 16:35

You could use a discretization method such as the finite element or finite difference method with a linearization technique such as Picard or Newton method to solve this problem. A similar question is this.

The main idea of these solution techniques is as follows:

• Pick an initial guess for the solution.
• Linearize your equation and write an updated solution in terms of a previous solution.
• Solve a sequence of linear problems until you achieve some convergence criterion.

I suggest that you check the following reference where this is explained step-by-setp.

• I guess here we are talking about doing a sequence of steps in x, using a nonlinear solver for each step. Using an implicit ODE integrator such as LSODE would do essentially the same. Apr 28 at 14:20
• @MaximUmansky, I don't think so. You linearize and then discretize your differential equation and solve the resulting linear system. Apr 28 at 15:41
• @nicoguaro Thank you for this, I am currently reading he reference you attached, hopefully it will work. Apr 29 at 3:46
• @nicoguaro Yes, you can solve it this way, for the whole profile at once, as a nonlinear problem. I would also add there artificial time evolution to make it easier to converge, but maybe for a 1D problem it is not needed. Apr 29 at 4:48

Are you sure that you have the correct boundary conditions? I think the equation you state for the analytical case, $$H=x/\tanh(x)-1$$ is incompatible with the boundary condition $$H(\bar{\epsilon})=\bar{\epsilon}$$ where $$\bar{\epsilon}$$ is large. Maybe this is why you are getting incorrect behaviour for large $$x$$?

I tried solving your equation for the $$\alpha=0$$ case with the fixed BC $$H(\bar{\epsilon})=\bar{\epsilon}-1$$ and was able to find good agreement to the analytic solution using a Newton-Gauss-Seidel relaxation method.

• Yes, $\bar{\epsilon} - 1$ is a better estimate. However, while using the inbuilt functions, as well as my written code, I gave exact analytical conditions. It still fails, at around x = 12. I do not think it is a stable IVP. That's great, could you please see if your method works for other cases? Say $\alpha = 0.1$. Also, what is a good resource to learn this method? Apr 30 at 4:50
• My method seems to work for $\alpha=0.1$, though I didn't attempt to verify my solution in any way (I don't think there is an analytic solution for $\alpha\neq0$ to compare against). Try looking up Newton iterative methods (the Gauss-Seidel part is just a small optimization). Here is a report that shows how to solve with Newton's method and finite differences: lakeheadu.ca/sites/default/files/uploads/77/docs/… Apr 30 at 6:16
• There is no analytical solution, however there is an integral of a function $\xi(K, H, K', H')$ which has the value 1.10-1.13. I can share the form of integral if you want to try it. Regardless, thanks a lot, I will see that report and try to implement it. May 1 at 3:50
• I am still stuck over implemeneting this. The crux is I know how to use it to construct a matrix for a 2nd order DE, but can't seem to modify it for coupled DEs. Could you guide regarding the structure of the matrix? May 10 at 4:25
• Yes, now it works. However, a word of caution: I tried shooting and it failed. Because shooting converts BVP to IVP and the IVP seems unstable. May 18 at 3:19