Consider a frictionless train that uses a constant propulsion force to accelerate but relies solely on the gravitational pull of the Earth and aerodynamic drag for deceleration. Using the simplifying assumption shown in Table 1, we compute the required propulsion force and the travel time for the train between two cities 1000km apart through a straight tunnel.

Train mass | m = 100,000 kg |

Distance travelled | d = 1000,000 m |

Earth Radius | Re = 6371,000 m |

Gravitational constant | g = 10 m^{2}/s |

Drag force | F(v) = 0.5 v^{2} N |

Propulsion force | F_{p} N |

Gravitational force | F_{g} = m_{t} g cos(θ) N |

Figure 1 shows the forces acting on the train during its trip from City A to City B along the earth circular path.

The motion for the train is governed by the 2^{nd} order equation:

where `θ(x)`, can be approximated by the following function based on the fact that `d/r << 1`

The train starts from rest so the initial values are:

$$x\left(0\right)=0$$ $$\dot{x}\left(0\right)=0$$Upon arrival to city B, the train would have travelled distance `d` and come to a stop so the final conditions can be stated as:

The unknown variables for this problem are `F _{p}` and

We treat the final conditions as constraints on the equation of motion then solve for the unknown variables `F _{p}` and

The first step is to convert the 2nd order equation into two first order equations as follows:

$$\frac{dx}{dt}=v$$ $$\frac{dv}{dt}=\frac{1}{m}\left({F}_{p}-{F}_{d}\left(\dot{x}\right)+{F}_{g}\left(\theta \left(x\right)\right)\right)$$
Next, we define the model formulas using the named variables
as shown in Table 2. The green, `B2:B4`, and yellow, `B11:B12`, ranges represent the system variables `(t,x,v)`, and the differential equations which are input arguments to the
initial value problem solver function, `IVSOLVE`. The blue, `D10:D11` range represents the design parameters
`F _{p}` and

A | B | C | D | |

1 | Variables with initial values | Parameters | ||
---|---|---|---|---|

2 | t | m | 100000 | |

3 | x | 0 | d | 1000000 |

4 | v | 0 | g | 10 |

5 | Re | 6371000 | ||

6 | Forces | |||

7 | Fg | =m*g*COS(theta) | theta | =IF(x<d/2,ATAN(Re/(d/2-x)),IF(x>d/2,ATAN(Re/(d/2-x))+PI(),PI()/2)) |

8 | Fd | =0.5*v^2 | ||

9 | Unknowns with initial guess | |||

10 | Equations | tf | 2000 | |

11 | dx/dt | =v | Fp | 1000 |

12 | dv/dt | =(Fp+Fd+Fg)/m |

We simulate the system with the guess value for a sufficient long time of 4500 seconds by executing the array formula
`=IVSOLVE(B11:B12, (t,x,v), {0,4500})`

in array `J3:L40` shown in Table 3 and plotted in Figure 2.

J | K | L | |

3 | t | x | v |

4 | 0 | 0 | 0 |

5 | 125 | 6115.463 | 96.66716 |

6 | 250 | 23618.42 | 180.3757 |

7 | 375 | 50345.74 | 243.5815 |

8 | 500 | 83614.12 | 285.3167 |

9 | 625 | 120912.9 | 308.8178 |

10 | 750 | 160244.8 | 318.6075 |

11 | 875 | 200165.6 | 318.7776 |

12 | 1000 | 239672.8 | 312.4788 |

13 | 1125 | 278105.9 | 301.8894 |

14 | 1250 | 315027.7 | 288.4672 |

15 | 1375 | 350145.5 | 273.1784 |

16 | 1500 | 383269.8 | 256.6382 |

17 | 1625 | 414271.2 | 239.2515 |

18 | 1750 | 443058.6 | 221.296 |

19 | 1875 | 469577.2 | 202.9315 |

20 | 2000 | 493781.3 | 184.2855 |

21 | 2125 | 515639.8 | 165.4439 |

22 | 2250 | 535135 | 146.4526 |

23 | 2375 | 552249.8 | 127.3546 |

24 | 2500 | 566971.2 | 108.1826 |

25 | 2625 | 579292.2 | 88.95507 |

26 | 2750 | 589207.6 | 69.68509 |

27 | 2875 | 596712.8 | 50.38592 |

28 | 3000 | 601804.6 | 31.06775 |

29 | 3125 | 604479.3 | 11.73845 |

30 | 3250 | 604737.5 | -7.59494 |

31 | 3375 | 602579.9 | -26.9254 |

32 | 3500 | 598007 | -46.2459 |

33 | 3625 | 591021.5 | -65.5461 |

34 | 3750 | 581622.5 | -84.8244 |

35 | 3875 | 569814.4 | -104.068 |

36 | 4000 | 555605.5 | -123.255 |

37 | 4125 | 539004.1 | -142.368 |

38 | 4250 | 520019.3 | -161.382 |

39 | 4375 | 498664.3 | -180.264 |

40 | 4500 | 474961 | -198.956 |

In Table 4, we define two constraints corresponding to the final conditions
$x\left(\mathrm{{t}_{f}}\right)=d$
and
$\dot{x}\left({t}_{f}\right)=0$
based on the solution obtained in Step 1.
Constraint `C10` computes the difference between the interpolated value for `x(t _{f})`
and the distance

C | |

10 | =ODEVAL(J4:J40, x, "INTERP", tf) - d |

11 | =ODEVAL(J4:J40, v, "INTERP", tf) |

What is `ODEVAL`? `ODEVAL` is one of the criterion functions that allow us to define dynamic constraints on the solution array.
Here we are simply computing values `x(t _{f})` and

Using the solver `NLSOLVE`, we solve for the unknowns `t _{f}` and

`=NLSOLVE(C10:C11, (Fp, tf))`

in array A | B | |

10 | Fp | 62648.94 |

11 | tf | 3863.575 |

12 | SSERROR | 1E-16 |

13 | ITRN | 6 |

14 | TIME (s) | 0.459 |