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?

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

anuttarasammyak said:

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?

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 ?

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

anuttarasammyak said:

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

renormalize said:

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.

anuttarasammyak said:

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

Ax_xiom said:

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.

pasmith said:

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 said:

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

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.

anuttarasammyak said:

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 said:

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.

Ax_xiom said:

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

Ax_xiom said:

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?

After some thinking, I don't think this is possible. For an arbitrary stage i, the expression for the mass ratio will be this $$N_i = \frac{m_i + \alpha_i p + U}{s m_i + \alpha_i p + U}$$ where U is the mass of all upper stages. And we (ideally) want to find ##\frac{m_i + \alpha_i p + U}{U}## in terms of ##N, \alpha_i## and ##s##. We only have one equation and need to eliminate 2 variables (##m_i, p##) so one will be free, unless we involve equations from other stages.

So stage 2 will have to involve an equation from stage 1 (which would be the final stage, I'm denoting these in reverse order for ease). Stage 3 will have to involve an equation from stage 2, which would involve an equation from stage 1. This continues with each successive stage as the equations get more and more complex.