Introduction

Having very briefly reviewed the fundamentals of mechanics in the previous article , let's move on to practice, for even the brief theory that has been considered will suffice.



So, the task:

The stone is thrown vertically, with no initial speed from a height of h = 100 m. Neglecting air resistance, determine the law of stone movement as a function of the height of the stone above the Earth’s surface from time to time. Acceleration of free fall to be equal to 10 m / s ²

Simple task? Yes, an elementary, having an analytical solution that is easy to write for a little bit literate schoolchild. But this simple task will serve as a very illustrative example.

1. Task formalization

What starts any simulation. Under the formalization understand the receipt of mathematical expressions describing the process in progress. At the same time, assumptions are formulated: the list of simplifications of the model due to factors whose influence can be neglected.
')
For this problem, the assumption is applicable, according to which a stone can be considered a material point. A single force is applied to this point - the force of gravity, so we use the assumption of the absence of air resistance.

The next assumption is to assume that the Earth is flat, since the height from which we throw a stone is negligible in comparison with the size of the planet, and therefore the curvature of its surface can be neglected. Then gravity can be considered constant, directed perpendicular to the surface of the Earth and equal in magnitude

$$

$$P = m \, g$$

where g is the acceleration of free fall at the surface of the Earth. Now it's time to make the equations of motion for the stone. Remember these equations

$$

$$\ begin {cases} & m \, \ ddot x = \ sum F_ {kx} \\ & m \, \ ddot y = \ sum F_ {ky} \\ & m \, \ ddot z = \ sum F_ {kz} \ end {cases}$$

Their left part does not interest us so far, but in the right one there are sums of projections of forces applied to a point on the coordinate axes. Let the x and y axes be located on the surface, and the z axis is directed upward perpendicular to it. There is only one force, its projections on the x and y axes are zero, and on the z axis, the projection is negative, since the force is directed against the axis direction, i.e.

$$

$$\ begin {cases} & m \, \ ddot x = 0 \\ & m \, \ ddot y = 0 \\ & m \, \ ddot z = -m \, g \ end {cases}$$

The mass of the stone is obviously not zero, so you can safely divide both sides of the resulting equations into it.

$$

$$\ begin {cases} & \ ddot x = 0 \\ & \ ddot y = 0 \\ & \ ddot z = -g \ end {cases}$$

I will not bore, proving that the movement of the stone will occur strictly vertically, although this must be done from a formal point of view. Zero on the right side of the first two equations does not at all mean that it is impossible to move along these axes - we recall Newton's first law. I will dwell on this in the next article in more detail, but for now it’s fair to set the one-dimensionality of motion by writing the final differential equation

$$

$$\ ddot z = -g$$

The fact that we got not a lot is not enough - a mathematical model of the process of what is happening in the problem. Pathetic, yes?

Not. Analyzing this equation we conclude, for example, that the mass of a stone does not affect the law of its movement, because there is no mass in this equation. You see, even without solving the equations, we have already formally proved the validity of the experience with a feather and a piece of lead in a vacuum, which they like to show in school (and some have repeated it on the Moon ).

It’s easy to get an analytical solution, I’m not even going to bother, it’s

$$

$$z (t) = h - \ frac {g \, t ^ 2} {2}$$

But how to solve this numerically? And what is it all about - “numerically”?

2. Numerical integration of a first order differential equation

What is the first order? I said last time that the equations of motion are of the second order. Everything is correct, but most methods for solving diffur on a computer are only able to solve first-order equations. There are methods of direct integration of second-order equations (for example, the Verlet method), but not about them now.

First, this equation is of a type that allows for a decrease in the order. The right side does not depend on the unknown function (there is no z), so we recall that

$$

$$\ ddot z = \ dot v_z$$

the projection of the acceleration on the z axis is equal to the first derivative of the velocity projection on the same axis z. Well, great, then

$$

$$\ dot v_z = -g$$

here is the first order equation. Not always this number passes (I will not now about the form of Cauchy!), But in this case everything is in order. We will search for not the coordinate but the speed of the point. What's next? So what is next

$$

$$\ frac {dv_z} {dt} = -g$$

for the derivative, we know, is the ratio of an infinitely small increment of a function (velocity) in to the argument that caused it an infinitesimal increment of argument (time). Take a very small amount of time. $$$\ Delta t$so small that can be considered

$$

$$\ frac {dv_z} {dt} \ approx \ frac {\ Delta v_z} {\ Delta t}$$

What is the result? And that's what

$$

$$\ Delta v_z = -g \, \ Delta t$$

We got a speed increment. Negative increment. As it is, a stone falling down will accelerate the same! Yes, there will be. Its speed, its velocity vector, will be directed downwards. So the projection of this vector on the z axis will be negative. That's right, we get an absolute increase in the projection of the downward vector. We know the initial value of the speed is zero, which means

$$

$$v_z ^ {(0)} = 0$$

Using the fact that we can calculate the increment of speed, we calculate what the speed will be say in 0.1 seconds

$$

$$v_z ^ {(1)} = v_z ^ {(0)} - g \, \ Delta t = 0 - 10 \ cdot 0.1 = -1.0$$

and after another 0.1 seconds

$$

$$v_z ^ {(2)} = v_z ^ {(1)} - g \, \ Delta t = -1.0 - 10 \ cdot 0.1 = -2.0$$

and after another 0.1 seconds

$$

$$v_z ^ {(3)} = v_z ^ {(2)} - g \, \ Delta t = -2.0 - 10 \ cdot 0.1 = -3.0$$

Hmm, so we can continue for quite a long time, but we limit ourselves to a period of one second.

Time with	Speed, m / s
0.0	0.0
0.1	-1.0
0.2	-2.0
0.3	-3.0
0.4	-4.0
0.5	-5.0
0.6	-6.0
0.7	-7.0
0.8	-8.0
0.9	-9.0
1.0	-10.0

That is, using the formula

$$

$$v_z ^ {(i + 1)} = v_z ^ {(i)} - g \, \ Delta t$$

we got the point velocity versus time. And all you need to do is take the value of speed at the current time, and add to it the increment that the speed will receive at the new time, which is from the current $$$\ Delta t$seconds The time increment here is called the integration step . And we calculate the increment as the value of the derivative of the desired function at the current time multiplied by the step. Simply? Yes, of course. And the formula that I wrote has a name — an explicit Euler method for numerically solving differential equations . This is the so-called recurrent formula, where the new value of the calculated value depends on its previous value.

And what about the height of the point above the Earth? Yes, all the same, see

$$

$$\ frac {dz} {dt} = v_z$$

for the projection of speed is a derivative of the corresponding coordinate. We apply the formula of the Euler method for this equation, because the speed we already know

$$

$$z ^ {(i + 1)} = z ^ {(i)} + v_z ^ {(i)} \, \ Delta t$$

and according to this formula we will add one more column to our table

$$

$$\ begin {align} & z ^ {(1)} = z ^ {(0)} + v_z ^ {(0)} \, \ Delta t = 100 + 0.0 \ cdot 0.1 = 100 \\ & z ^ {(2 )} = z ^ {(1)} + v_z ^ {(1)} \, \ Delta t = 100 - 1.0 \ cdot 0.1 = 99.9 \\ & z ^ {(3)} = z ^ {(2)} + v_z ^ {(2)} \, \ Delta t = 99.9 - 2.0 \ cdot 0.1 = 99.7 \\ \ end {align}$$

and so on

Time with	Speed, m / s	Height, m
0.0	0.0	100.0
0.1	-1.0	100.0
0.2	-2.0	99.9
0.3	-3.0	99.7
0.4	-4.0	99.4
0.5	-5.0	99.0
0.6	-6.0	98.5
0.7	-7.0	97.9
0.8	-8.0	97.2
0.9	-9.0	96.4
1.0	-10.0	95.5

Hmm, well, firstly, it is noticeable that the height varies with us unevenly, since the speed changes with time. Now our derivative itself depends on time. But already in the first step, we notice something is wrong - the speed is already there, but the height is still 100 meters. How so?

This happened because at each step we set the derivative (speed) constant. The method does not provide information about what happens with the solution between steps. Accordingly, the error accumulates, we compare the solution obtained with the exact

$$

$$z (t) = 100 - 5 \, t ^ 2$$

Time with	Speed, m / s	Height, m	Exact solution, m
0.0	0.0	100.0	100.00
0.1	-1.0	100.0	99.95
0.2	-2.0	99.9	99.80
0.3	-3.0	99.7	99.55
0.4	-4.0	99.4	99.20
0.5	-5.0	99.0	98.75
0.6	-6.0	98.5	98.20
0.7	-7.0	97.9	97.55
0.8	-8.0	97.2	96.80
0.9	-9.0	96.4	95.95
1.0	-10.0	95.5	95.00

Yes, our stone seems to hang in the air. The numerical solution lags behind the analytical one, and the further we think, the higher the counting error. The error accumulates, as at each step we take a more and more rough approximation. What to do?

First, you can reduce the step. Let's say 10 times, let $$$\ Delta t = 0.01$seconds

Time with	Speed, m / s	Height, m	Exact solution, m
0.0	0.0	100.0	100.00
0.1	-1.0	99.96	99.95
0.2	-2.0	99.81	99.80
0.3	-3.0	99.57	99.55
0.4	-4.0	99.22	99.20
0.5	-5.0	98.78	98.75
0.6	-6.0	98.23	98.20
0.7	-7.0	97.59	97.55
0.8	-8.0	96.84	96.80
0.9	-9.0	96.00	95.95
1.0	-10.0	95.05	95.00

Already better, the error at the end of the account does not exceed 0.05 meters, and this is 10 times less than the previous value. It can be assumed that by reducing the step another 10 times, we will get an even more accurate solution. I cheated by displaying values ​​for only 10 points with a step of 0.1, in fact, to get such a table, 100 iterations are needed instead of 10. At step 0.001, a thousand iterations will be needed, and the result will be

Time with	Speed, m / s	Height, m	Exact solution, m
0.0	0.0	100.0	100.00
0.1	-1.0	99.9505	99.95
0.2	-2.0	99.8010	99.80
0.3	-3.0	99.5515	99.55
0.4	-4.0	99.2020	99.20
0.5	-5.0	98.7525	98.75
0.6	-6.0	98.2030	98.20
0.7	-7.0	97.5535	97.55
0.8	-8.0	96.8040	96.80
0.9	-9.0	95.9545	95.95
1.0	-10.0	95.0050	95.00

If you have tried to perform these calculations manually, then you understand now how monotonous and time consuming they are if you need high accuracy. That is why the flowering of numerical simulation coincided with the advent of computers. They are just needed in order to quickly perform many uniform operations on numbers.

The Euler method is the simplest known method for integrating differential equations. From our simple example, it can be seen that the error of the method is directly proportional to the step of integration, and this is true. Such methods are called methods of the 1st order of accuracy.

The accuracy of calculations, even at step 0.1, can be improved if we apply, say, the 4th order Runge-Kutta method. But this is a separate story.

Conclusion

Think about it ... We looked at a very simple example. We did not even use a computer, but we already understand the principle by which the most powerful supercomputers in the world work, which model the early stages of the life of the Universe. Of course, everything is much more complicated, but the principle is the same.

Imagine what a powerful tool you have in your hands. This last article where we will not use a computer. I promised Octave. Next time it will be him.