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

 

[Ax_xiom]

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.$$




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.

So not related to this specific topic but related to this project. So turns out the equations I was working with are likely wrong. So the method I'm using starts like this: $$
\Delta V = c_1 \ln(\frac{m_1 + m_2 + (\alpha_1 + \alpha_2 + 1)p}{s m_1 + m_2 + (\alpha_1 + \alpha_2 + 1)p})
+ c_2 \ln(\frac{m_2 + (\alpha_2 + 1)p}{s m_2 + (\alpha_2 + 1)p}) $$

Where $m_i$ is the mass of fuel tanks in the ith stage of the rocket, $p$ is the mass of the payload, $\alpha_i$ is the mass of any ancillary mass (engines, aerodynamic parts etc.) of the ith stage, divided by the payload mass, $s$ is the fraction of the fuel tank that is dry mass and $\Delta V$ is the target Delta-V. The method continues with stating: $$
\begin{align}
N_1 = \frac{m_1 + m_2 + (\alpha_1 + \alpha_2 + 1)p}{s m_1 + m_2 + (\alpha_1 + \alpha_2 + 1)p}\nonumber \\
N_2 = \frac{m_2 + (\alpha_2 + 1)p}{s m_2 + (\alpha_2 + 1)p} \nonumber
\end{align} $$
And having us try to find: $$
\begin{align}
\frac{m_1 + m_2 + (\alpha_1 + \alpha_2 + 1)p}{m_2 + (\alpha_2 + 1)p}\nonumber \\
\frac{m_2 + (\alpha_2 + 1)p}{p} \nonumber
\end{align} $$
In terms of $N_1, N_2$ and our other dimensionless variables. I will be calling these the "stage ratios" as these are the ratio of mass between each stage. Unfortunately, it seems that doing this creates equations that get more and more complex with each stage. For the second stage, the stage ratio is this: $$\frac{(1-s)N_2}{1 - s N_2}(1+\alpha_2) $$
For the first stage, the stage ratio is this: $$ \frac{(1-s)N_1}{1 - s N_1}(1+\alpha_1(\frac{m_2+(\alpha_2+1)p}{p})^{-1}) $$
Which involves the stage ratio for stage 2. If we had 3 stages and wanted to find the stage ratio for stage 1, this would be the expression: $$\frac{(1-s)N_1}{1 - s N_1}(1+\alpha_1(\frac{m_3+(\alpha_3+1)p}{p}\frac{m_2+m_3+(\alpha_2+\alpha_3+1)p}{m_3+(\alpha_3+1)p})^{-1}) $$ Which involves the stage ratios for stage 2 and 3. This will get increasingly complex for more and more stages. So is there a better way to approach this problem or should I settle for only allowing the user to specify the structural ratio?

 

[anuttarasammyak]

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:

We should solve the equation of $\lambda$
$$\sum_{j=1}^n c_j\ln(1-\frac{1}{c_j \lambda}) - \ln s \sum_{j=1}^n c_j - V_f=0$$
N's are given by solution $\lambda$ as
$$ N_j=\frac{1}{s}(1-\frac{1}{c_j \lambda}) $$
We can draw the graph of LHS, $y=f(\lambda)$ by computer to observe whereabout y=0. If we can expect $$ \frac{1}{c_j \lambda} << 1 $$ for all j, $\lambda_0$ satisfying
$$-\frac{n}{\lambda_0}- \ln s \sum_{j=1}^n c_j - V_f=0$$
is a good initial point for numerical methods.
[EDIT] I corrected careless mistakes.

 

[Ax_xiom]

We should solve the equation of $\lambda$
$$\sum_{j=1}^n \ln(1-\frac{1}{c_j \lambda}) - \ln s \sum_{j=1}^n c_j - V_f=0$$
N's are given by solution $\lambda$ as
$$ N_j=\frac{1}{s}(1-\frac{1}{c_j \lambda}) $$
We can draw the graph of LHS, $y=f(\lambda)$ by computer to observe whereabout y=0. If we can expect $$ \frac{1}{c_j \lambda} << 1 $$ for all j, $\lambda_0$ satisfying
$$-\frac{1}{\lambda_0}\sum_{j=1}^n \frac{1}{c_j } - \ln s \sum_{j=1}^n c_j - V_f=0$$
is a good initial point for numerical methods.

This is nice, but @pasmith pointed out a way to make the Jacobian upper triangular, which makes the computation of $\Delta X_n$ much easier to scale and makes everything much more convenient. The bigger problem I am currently facing is that the system of equations I initially came up with was wrong and the correct versions being extremely complex

 

[anuttarasammyak]

Thanks for the comment. I should appreciate it if you share values of setting parameters : $c_j$'s , s and $V_f$. I would like to go to physics from mathematics.

 

[Ax_xiom]

Thanks for the comment. I should appreciate it if you share values of setting parameters : $c_j$'s , s and $V_f$. I would like to go to physics from mathematics.

the $c_i$'s and $V_f$ will be variable but the $c_i$'s will typically be around $2400-3000ms^{-1}$ and $V_f$ will be anywhere from $2600 - 9000ms^{-1}$

 

[anuttarasammyak]

Thanks. How about s ?

Say s=0.5, please find attached my calculation in EXCEL sheet where I try inputting $K=\frac{1}{\lambda}$ so that f(K) becomes around zero.

Is N>1 OK ? Does these results have a sense of physics?
Some graphs of y=f(x)

For s=0.1

 

[Ax_xiom]

Thanks. How about s ?

Say s=0.5, please find attached my calculation in EXCEL sheet where I try inputting $K=\frac{1}{\lambda}$ so that f(K) becomes around zero.

View attachment 372959

Is N>1 OK ? Does these results have a sense of physics?

s will be around 0.1. These results make sense but this is a solved issue currently.

Which involves the stage ratios for stage 2 and 3. This will get increasingly complex for more and more stages. So is there a better way to approach this problem or should I settle for only allowing the user to specify the structural ratio?

This is the issue I am currently having