Using Newton-Raphson method to solve the hydrostatic equation

I'm trying to use newton-raphson method for nonlinear systems of equations as described in 'Numerical recipies' book in chapter 9.6 to solve the hydrostatic equation for a polytropic star.

For each iteration, I change the radius vector as described in the book like x→x+Δx(Δ can be calculated for each iteration) but I can't ignore the fact that it can do shell crossing - meaning the i shell moves to a higher shell then i+1i+1, this cause a negative volume. I know there is a method to avoid shell crossing by multiplying delta with a factor but I don't know how to calculate it. There is no description in the book about shell crossing.

let me be more specific:

I want to solve $F_i(x_1,....x_N)=0$ set of equations,( when i goes from 1 to N). the method is to expand $F$ to a taylor series: $\vec{F}(\vec{x}+\vec{\delta x})=\vec{F}+J\vec{\delta x} +O(x^2)$ , where $J$ is the partial deriviatives of $F$. by setting $\vec{F}(\vec{x}+\vec{\delta x})=0$ we can get N sets of equations : $J \cdot \vec{\delta x}=-F$ that can be solved by LU decomposition.

the method starts from initial point $x_{old}$. for each step we calculate $\delta x$ and add it to $x_{old}$ , meaning : $x_{new}=x_{old}+\delta x$

if the newton step $\delta x$ minimizes $f=0.5\cdot F\cdot F$ we can say this step is acceptible. By using line searches and backtracking we can guaranty to minimize $f$. Be aware that every solution to $F_i(x_1,....x_N)=0$ is a minimum of $f$ but not the oposite, so we must check if $f$ goes to zero.

My set of equations are: $F_i=\frac{dP_i/dm_i}{Gm_i/4\pi r_i^4}+1=0$

where $\delta P_i=P_{i-1}-P_i$ , $\delta m_i =0.5\cdot (m_{i+1}-m_{i-1})$ , and from the polytropic equetion: $P_i=K(\frac{m_i-m_{i-1}}{4\pi / 3 (r_i^3-r_{i-1}^3)})^\gamma$

Where do i stuck? in the first iteration. and it's suppose to be bad. the problem is shell crossing. $\delta x$ changes $r$ so much that it changes the order of the cells. meaning, $r_i$ the radius of the star (from zero) is suddenly bigger then $r_{i+1}$. I know there is a method to correct this by multiplying $\delta$ by some factor. How to calculate it? Is there any other methods to correct it?

Do anyone know how to help?

• If you share a complete set of your equations then maybe? – zhk May 6 '17 at 3:51
• I added my equations now – Noam Chai May 6 '17 at 6:37
• Is this a discretization of a BVP? Then I don't understand why your $r_i$ are not fixed. The way you've written it down, you have $2N$ variables $(m_i,r_i)$ and just $N$ equations. – Kirill May 6 '17 at 15:37