Best way to numerically solve this system of equations?

[Ax_xiom]

TL;DR

I have a system of equations that I need solving numerically

So I'm working on a project where I am trying to work out what the ideal sizes of rocket stages are and I am using Excel to allow the user to interact with this quickly. The method I am using is derived from this video, and I end up with this system of equations: \begin{align}
1 - c_{1}\lambda + \frac{sN_{1}}{1 - sN_{1}} = 0 \\
1 - c_{2}\lambda + \frac{sN_{2}}{1 - sN_{2}} = 0 \\
c_{1}\ln(N_{1}) + c_{2}\ln(N_{2}) - V_{f} = 0 \end{align}
Where I am trying to solve for $N_{1}$ and $N_{2}$ The number of equations like (1) and (2) can be as many as 5, and the number of $c_{n}\ln(N_{n})$ terms in (3) is also variable. As far as I know, these equations cannot be solved analytically, so they must be solved numerically. Here are the options that I've considered:

1. Newton's method with the entire system
   • Pros:
     • Easiest to think about conceptually
     • Converges quickly from a reasonable starting point
   • Cons:
     • Hard to implement a variable amount of stages
     • More work to set up
     • Have to consider the order of the columns and rows of the Jacobian
2. Solve for $\lambda$ and $N_{n}$ after
   • It is possible to show that $N_{n} = \frac{1-\frac{1}{c_{n}\lambda}}{s}$. Using this, it is possible to substitute these into equation (3) and numerically solve for $\lambda$ before solving for $N_{n}$ using the aforementioned expression.
   • Pros:
     • Less computation
     • Easier to adjust for different amounts of stages
   • Cons:
     • The resulting expression is extremely badly behaved, diverging to $-\infty$ when using Newton-Raphson unless given an extremely close guess

So what method should I use? Have I overlooked any methods that I may be better off using?

 

[anuttarasammyak]

I would try let computer draw three graphs (1)(2)(3) and observe crossing points, if not ask AI the solutions direct.

 

[Ax_xiom]

I would try let computer draw three graphs (1)(2)(3) and observe crossing points.

A few problems:

1. I want to implement up to 5 stages, which would result in 6 equations and 6 variables that would need to be graphed which would require a 6D space to visualise
2. I'm not sure how easy that would be to use or how accurate it would be
3. Can Excel even draw 3D graphs?

 

[anuttarasammyak]

There are six kind of alphabets in your equations. 3 equations for 6 unknowns seems desperate. Which are numbers given for numerical calculation and which are parameters for which we would like to get solutions ?

 

[renormalize]

Which are numbers given for numerical calculation and which are parameters for which we would like to get solutions ?

From skimming through the video referenced by the OP, the known rocket parameters are:

• $c_n$ is the exhaust velocity of the $n^\text{th}-$stage rocket motor.
• $s_n$ is the $n^\text{th}-$stage structural factor (fraction of the total stage-mass that is not fuel).
• $V_f$ is the desired final velocity of the last stage.

Unknowns to be solved for:

• $N_n$ are mass ratios of the various stages. These values ultimately determine the mass of fuel in each stage that's required to achieve the final velocity.
• $\lambda$ is a Lagrange-multiplier parameter.

But the OP is essentially reinventing the wheel: there's a nice online calculator that already does what they want: https://space.geometrian.com/calcs/opt-multi-stage.php.

 

[anuttarasammyak]

Following #5, we get equation for $\lambda$
$$c_1 \ln(1-\frac{1}{c_1 \lambda})+c_2 \ln(1-\frac{1}{c_2 \lambda})=V_f+(c_1+c_2)\ln s$$
if not mistaken. We can solve it to get solution $\lambda$. For an example case of
$$c_1=c_2=s=V_f=1$$ ,though I do know these values are physical or not, we get
$$\lambda=-\frac{1}{\sqrt{e}-1}$$

 

[Ax_xiom]

There are six kind of alphabets in your equations. 3 equations for 6 unknowns seems desperate. Which are numbers given for numerical calculation and which are parameters for which we would like to get solutions ?

There's 3 equations in my post, but that's the version for 2 stages. I want to be able to expand it up to 5 stages where the system of equations will be this: $$\begin{align}

1 - c_{1}\lambda + \frac{sN_{1}}{1 - sN_{1}} = 0 \\

1 - c_{2}\lambda + \frac{sN_{2}}{1 - sN_{2}} = 0 \\

1 - c_{3}\lambda + \frac{sN_{3}}{1 - sN_{3}} = 0 \\

1 - c_{4}\lambda + \frac{sN_{4}}{1 - sN_{4}} = 0 \\

1 - c_{5}\lambda + \frac{sN_{5}}{1 - sN_{5}} = 0 \\

c_{1}\ln(N_{1}) + c_{2}\ln(N_{2}) + c_{3}\ln(N_{3}) + c_{4}\ln(N_{4}) + c_{5}\ln(N_{5}) - V_{f} = 0 \end{align}$$

But the OP is essentially reinventing the wheel: there's a nice online calculator that already does what they want: https://space.geometrian.com/calcs/opt-multi-stage.php.

That's nice, but instead of specifying the structural ratio, I want to be able to specify the dry mass, non-fuel tank mass of the stage. It just turns out that doing it that way results in equations equivalent to the video.

if not mistaken. We can solve it to get solution λ. For an example case of

The value of s is the structural ratio of the fuel tanks. If we set it as 1, that means the fuel tank is just dry mass so those values aren't physical at all

 

[pasmith]

TL;DR: I have a system of equations that I need solving numerically

So I'm working on a project where I am trying to work out what the ideal sizes of rocket stages are and I am using Excel to allow the user to interact with this quickly. The method I am using is derived from this video, and I end up with this system of equations: \begin{align}
1 - c_{1}\lambda + \frac{sN_{1}}{1 - sN_{1}} = 0 \\
1 - c_{2}\lambda + \frac{sN_{2}}{1 - sN_{2}} = 0 \\
c_{1}\ln(N_{1}) + c_{2}\ln(N_{2}) - V_{f} = 0 \end{align}

It may be easier to express the first two as $$
(1 - sN_i)(1 - c_i \lambda) + sN_i = 1 - c_i\lambda + sc_iN_i \lambda = 0.$$

Where I am trying to solve for $N_{1}$ and $N_{2}$ The number of equations like (1) and (2) can be as many as 5, and the number of $c_{n}\ln(N_{n})$ terms in (3) is also variable. As far as I know, these equations cannot be solved analytically, so they must be solved numerically. Here are the options that I've considered:

Newton's method is actually fairly straightforward here; we can even solve the resulting linear system at design time.

If we order the variables as $(N_1, \dots, N_5, \lambda)$ and the equations as $$
f_i = \begin{cases} (1 - c_i \lambda) + sc_iN_i\lambda & i = 1, \dots, 5 \\
\sum_{i=1}^5 c_i \ln N_i - V_f & i= 6 \end{cases}$$ then the Jacobian is $$
\begin{pmatrix} sc_1 \lambda & 0 & 0 & 0 & 0 & c_1(sN_1 - 1) \\
0 & sc_2 \lambda & 0 & 0 & 0 & c_2(sN_2 - 1) \\
0 & 0 & sc_3 \lambda & 0 & 0 & c_3(sN_3 - 1) \\
0 & 0 & 0 & sc_4 \lambda & 0 & c_4(sN_4 - 1) \\
0 & 0 & 0 & 0 & sc_5 \lambda & c_5(sN_5 - 1) \\
c_1/N_1 & c_2/N_2 & c_3/N_3 & c_4/N_4 & c_5/N_5 & 0 \end{pmatrix}$$ This is almost upper triangular, and all we need do to make it so is to multiply the last row by $s\lambda$ and subtract each other row $i$ divided by $N_i$ to get $$
\begin{pmatrix}
sc_1 \lambda & 0 & 0 & 0 & 0 & c_1(sN_1 - 1) \\
0 & sc_2 \lambda & 0 & 0 & 0 & c_2(sN_2 - 1) \\
0 & 0 & sc_3 \lambda & 0 & 0 & c_3(sN_3 - 1) \\
0 & 0 & 0 & sc_4 \lambda & 0 & c_4(sN_4 - 1) \\
0 & 0 & 0 & 0 & sc_5 \lambda & c_5(sN_5 - 1) \\
0 & 0 & 0 & 0 & 0 & - \sum_{i=1}^5 c_i \left(s- \frac 1{N_i} \right) \end{pmatrix}
\begin{pmatrix} \Delta N_1 \\ \Delta N_2 \\ \Delta N_3 \\ \Delta N_4 \\ \Delta N_5 \\ \Delta \lambda \end{pmatrix} = - \begin{pmatrix} f_1 \\ f_2 \\ f_3 \\ f_4 \\ f_5 \\ s \lambda f_6 - \sum_{i=1}^5 \frac{f_i}{N_i} \end{pmatrix}$$ which is easily solved by back substitution for the increments. The cases for fewer than 5 stages can be deduced from this; they all follow the same pattern.

 

[Ax_xiom]

Newton's method is actually fairly straightforward here; we can even solve the resulting linear system at design time.

Do you mean it's possible to solve the linear system without using Excel's matrix operations each time?

If that is correct, I believe these would be the solutions to the system: $$

\begin{align}
\Delta\lambda = \frac{s\lambda f_6 - \sum_{i=1}^5 \frac{f_i}{N_i}}{\sum_{i=1}^5 c_i (s - \frac{1}{N_i})} \nonumber \\
\Delta N_i = \frac{c_i(1-s N_i)\Delta\lambda - f_i}{s c_i \lambda} \nonumber
\end{align}$$

If these are correct, I think this would make everything much easier

This is almost upper triangular, and all we need do to make it so is to multiply the last row by sλ and subtract each other row i divided by Ni to get

In general, what kind of operations are you allowed to do on matrices of linear systems to make them easier to solve?

Edit: I believe it's just the same types of operations you can do on linear systems of equations (like adding and subtracting them, and scaling them by a parameter