Python numerical solution for a nonlinear second order ODE with two boundary conditions

2.6k Views Asked by At

I want to solve numerical the next equation, in Python

$$u''(x) = \left( a - \Big(b\big(u(x)^{2}\big)\Big) \right) \big(u'(x)\big)^{3}$$

it is a nonlinear second order $ODE$ with two $B.C$. $(u(x=0)=m_1, du/dx(x=N)=m_2)$ the x-axis have$ N$ gridpoints ($a, b$ 1D $N$-th arrays, not a number)

Do someone can help me to solve this equation numerically in Python?

Theoretically if I use the central difference method for $u'', u' (j=1,N)$ the next step will be to use the Newton iteration method for find a solution for each $j$. But the calculation of its Jacobian for this Nth system is something that I want to avoid.

Please, any suggestion ? from someone that have solve a similar problem...

1

There are 1 best solutions below

0
On

I'm not sure why you want to avoid computing the Jacobian. Granted, it is a bit messy, but it will probably give you the best method.

As you say, after central differences you get a nonlinear system of equations. You can try solving it with scipy.optimize.roots (which estimates the Jacobian itself).

A completely different method for solving your problem is shooting. The idea is that you solve the ODE as an initial value problem from $x=0$ to $x=N$ with initial conditions $u(0) = m_1$ and $du/dx(0) = \alpha$ for some value of $\alpha$, compare the resulting value of $du/dx(N)$ with the desired value $m_2$, and vary $\alpha$ in order to hit $m_2$. This is (in my opinion) simpler to program, but depending on the stability properties of the ODE it may not work if $N$ is large, or it may only work if you integrate the ODE from $x=N$ to $x=0$.