Learning Numerical Algorithms with Playground

In middle school, there was nothing more terrifying than math drawings.

Many problems have no ready-made analytical solutions, and some problems can be solved but require a lot of computational work. In this case, you need algorithms to help.

I hope you didn’t throw up on this algorithm. If you did, look on the bright side: you can have another lunch :]

In this tutorial, you will learn about algorithms and how to use them to solve problems that are difficult to analyze. With the playground, it is easy to visualize the solution and make it easy to see.

Even if you are not a math lover, or have little interest in physics or computer science, you will still find something interesting in this tutorial. All you need is some knowledge of calculus and basic physics.

You will learn how to solve two difficult problems using numerical algorithms. However, for the sake of clarity, both problems can also be solved analytically. Algorithms are ideal when analytical methods don't work, but even so, by comparing the two methods, it will be easier to understand their essence.

What are numerical algorithms?

Simply put, numerical algorithms are a class of methods for solving mathematical problems that do not rely on closed-form analytical solutions.

The question is, what is a closed-form analytical solution?

If there is a formula that allows us to use known numbers and a finite number of steps to obtain an accurate result, it is called a closed-form analytical solution.

To put it more simply: if you can use algebra to find an expression to solve a problem with an unknown quantity, and substitute certain known numbers to get the result, it means that you have obtained a closed-form analytical method.

When to use numerical algorithms?

Many problems have no ready-made analytical solutions, and some problems can be solved but require a lot of computational work. In this case, you need algorithms to help.

For example, suppose you need to write a physics engine that calculates the behavior of many objects in a finite amount of time. You can use numerical algorithms to get results faster.

This has the side effect that faster calculations mean the results will not be quite as accurate, but in most cases this is good enough.

Weather forecasting benefits from numerical calculations. The speed at which weather changes and the number of influencing factors are staggering. This is a very complex system, and only numerical simulation can accomplish the task of predicting the future.

It may be due to the lack of these algorithms that your iPhone always tells you that it is going to rain while the sun is shining outside.

[[142780]]

start

As a warm-up, let's play a game. Next, you will calculate the square root of a given number. Both methods will use the bisection method. The amazing thing is that you may already know this method, but you don't know its name.

Think back to when we were children, we may have played a game like this: choose a number within 100, and let the other person guess. You will only hint that the number he guessed is higher or lower.

The game starts. Xiao Ming tells Xiao Qiang to start guessing. Xiao Qiang says he guesses 1, Xiao Ming says it’s too small. Xiao Qiang guesses 2, Xiao Ming says it’s too small. Xiao Qiang then guesses 3, 4, and finally chooses 5, which is finally correct.

It guessed it right in 5 steps, which is good. But if Xiao Ming had chosen 78, it would have taken some time to guess it right.

If this game is played using the bisection method, the guessing speed will be much faster.

dichotomy

We know that the numbers are between 1 and 100, so instead of guessing one by one or just guessing blindly, we divide the numbers into two intervals: a = [1,50], b = [51, 100].

Next, we determine which interval the number is in. This is very simple. We just need to ask whether the number is 50. If the number is less than 50, then interval b is not considered. Next, we continue to divide interval a in half and repeat this step.

For example:

The number is 60. The intervals are: a = [1,50], b = [51, 100].

In the first step, Xiaoqiang said 50 (which is the upper limit of interval a), and Xiaoming said it was too small. At this time, Xiaoqiang knew that the number must be in interval b.

The second step is to decompose interval b into: c=[51,75], d=[76,100]. He still chooses the upper limit of interval c, 75, and Xiaoqiang tells him that it is too high. This means that the number must be in interval c.

Continue to decompose...

Using this method, you can get the correct answer in 7 attempts, while trying one by one would require 60 attempts.

1. 50 -> Smaller

2. 75 -> Bigger

3. 62 -> Bigger

4. 56 -> Smaller

5. 59 -> Smaller

6. 61 -> Bigger

7. 60 -> That’s right!

The process is similar to calculating the square root of x. The square root of a number x is between 0 and x. It can also be written as: (0, x]. If the number x is greater than or equal to 1, it can be written as: [1, x].

Decomposing the interval, we obtain a=(0, x/2], b=(x/2, x].

If the number x is 9, the interval is [1, 9], the decomposed interval is a=[1, 5], b=(5, 9], and the middle value m is (1+9)/2 = 5.

Next, we check if m * m - x is greater than the desired precision. Here, we check if m * m is greater than or less than or equal to x. If so, we continue checking using the interval (0, x/2], otherwise, we use the interval (x/2, x].

Let's take a look at the actual steps, initializing m=5, and expecting the accuracy value to be 0.1:

1. Calculate m * m - x: 5 * 5 - 9 = 11

2. Check if the result is less than or equal to the expected value: 25 - 11 <= 0.1?

3. Obviously not satisfied, continue to the next interval: Is m * m greater than 9? Yes. Next, use the interval [1, 5], and the test value is (1+5)/2=3.

4. Calculate m * m - x: ３ * 3 - 9 = 0

5. Check whether the expected value is met: 9 - 9 <= 0.1?

6. Done.

Note: You may be wondering what the brackets mean? Intervals have a fixed format (lower bound and upper bound). Different brackets represent different interval boundaries. Round brackets indicate that the boundary is outside the interval (i.e., open interval), while square brackets indicate that the boundary is within the interval (closed interval). For example, (0, a] includes a but does not include 0. In the above example: the interval (0, a/2] includes a/2 but does not include 0; the interval (a/2, a] represents all numbers greater than a/2 and less than a.

Using the Dichotomy Method on Playground

Now, it's time to apply the theory we learned. Let's implement the binary search algorithm ourselves. Create a new playground file and add the following code:

 func bisection(x: Double) -> Double {
 //1  
 var lower = 1.0  
 //2  
 var upper = x
 //3  
    var m = (lower + upper) / 2  
 var epsilon = 1e- 10  
 //4  
 while (fabs(m * m - x) > epsilon) {
 //5  
        m = (lower + upper) / 2  
 if m * m > x {
 upper = m
 } else {
 lower = m
 }
 } 
   
 return m
 }

1. Set the lower limit of the interval to 1

2. Set the upper limit of the interval to x

3. Find the middle value and define the expected accuracy

4. Check whether the operation can meet the accuracy requirements

5. If not satisfied, find a new middle value, define a new interval, and continue searching.

Add the following code to test the function:

 let bis = bisection( 2.5 )

To the right of the line m = (lower + upper) / 2, you can see that the code is executed 35 times, which means it takes 35 steps to find the answer.

#p#

Right away we can see the benefit of a lovely feature of the playground: the full history of values ​​is available.

The binary search method can successfully calculate an approximation of the actual result. By looking at the numerical history graph, you can see how the numerical algorithm approaches the correct solution.

Press option+cmd+enter to open the assistant editor and click the round button to the right of the m = (lower + upper) / 2 line to add history to the assistant editor.

You will see the method jump to the correct result little by little.

Classical mathematics can also be done

The next algorithm goes back to ancient times, it originated in ancient Babylon in 1750 BC and was described in Hero of Alexandria’s book Metrica 100 AD. That’s why it’s called “Hero’s method”. You can learn more about Hero at this link.

This method uses a function, and here you want to calculate the square root of a. If you can find the "zero point" of the function curve, where the function value is 0, then solving for the x value will give you the square root of a.

To do this, we first pick an arbitrary x value as the starting value, calculate the tangent value (the tangent line of the function) at that value, and then look at the 0 point on the tangent line (that is, the intersection of the line and the x-axis). We then use this 0 point as the starting value again, and repeat this process until we meet the accuracy requirement.

Because each time the tangent value is calculated, it will get closer to the true value of 0.

This process will gradually approach the true result. The following figure shows the square root of a when solving a=9, and the starting value is 1.

The starting point x0 = 1 generates a red tangent line, then the x1 point is generated, and then the purple line is generated; then the x2 point is generated, connected with the blue line, and finally the correct answer is found.

When classical mathematics meets Playground

We have a lot of good stuff that the ancient Babylonians didn't have, such as playground. Add the following code to the playground and check it out:

 func heron(x: Double) -> Double {
 //1  
 var xOld = 0.0  
 var xNew = (x + 1.0 ) / 2.0  
 var epsilon = 1e- 10   
   
 //2  
 while (fabs(xNew - xOld) > epsilon) {
 //3  
 xOld = xNew
    xNew = (xOld + x / xOld) / 2  
 } 
   
 return xNew
 }

What does this code do?

1. Create a variable to store the current result. xOld is the result of the last calculation, and xNew is the actual result.

• Finding a way to initialize xNew is a good starting point

• Epsilon is the desired accuracy

• In this example, we calculate the square root to 10 decimal places.

2. Use a while loop to check whether the desired accuracy is achieved.

3. If the accuracy cannot be achieved, set xNew to xOld and continue to the next iteration.

Use the following code to verify the function:

      
 let her = heron( 2.5 )

The Hero method only needs 5 iterations to find the correct result.

Click the rounded button to the right of the line xNew = (xOld + x/xOld)/2 to add the value history and you can see that a good approximation is found in the first iteration.

Simulating the motion of a harmonic oscillator

In this section, we look at how to use numerical integration algorithms to simulate the motion of a simple harmonic system - a basic type of dynamical system.

This system can describe many phenomena, such as the swing of a pendulum or the vibration of a spring. In particular, it can describe the scene where an object moves a certain amount of displacement over a period of time.

In this example, we assume that there is a mass attached to a spring. To simplify the problem, we ignore the drag and gravity. Therefore, the only force acting is the elastic force of the spring, which pulls the object back to its original position.

Under these assumptions, you only need to use two laws of mechanics to deal with the problem:

• Newton's second law of motion, which describes the relationship between the forces acting on an object and its acceleration.

• Hooke's law, which describes the proportional relationship between the spring force and the displacement of an object. Here k is the spring coefficient and x is the displacement of the object.

Simple Oscillator Equation

Since the elastic force is the only force acting on the object, we can rearrange the two equations above into:

Simplified writing:

It is also recorded as , which is the square of the resonant frequency.

A more precise equation is:

Where A is the amplitude, which is the offset of the object, called the phase difference. Both values ​​are set to fixed values ​​when initialized.

If you set time t = 0, then, and, you can easily calculate the amplitude and phase difference:

Let's look at an example. We have an object with a mass of 2kg connected to a spring with a spring constant of k=196N/m. At the initial time t=0, the spring moves 0.1 meters. We can calculate the amplitude, phase difference, and resonant frequency using the formula.

Try this out in a Playground:

Before we can use this formula to calculate the value at any point in time, we need to add some code.

Go back to the playground file and add the following code at the end:

 //1  
 typealias Solver = (Double, Double, Double, Double, Double) -> Void 
   
 //2  
 struct HarmonicOscillator {
 var kSpring = 0.0  
 var mass = 0.0  
 var phase = 0.0  
 var amplitude = 0.0  
 var deltaT = 0.0   
   
  init(kSpring: Double, mass: Double, phase: Double, amplitude: Double, deltaT: Double) {
 self.kSpring = kSpring
 self.mass = mass
 self.phase = phase
    self.amplitude = amplitude
 self.deltaT = deltaT
 } 
   
 //3  
  func solveUsingSolver(solver: Solver) {
    solver(kSpring, mass, phase, amplitude, deltaT)
 }
 }

What is happening in these code blocks?

1. Define a typealias of a function type. The function has five parameters of type Double and returns void.

2. Create a struct to define the resonator.

3. This function simply creates a Clouseur to solve the vibration problem. (There is no actual calculation code)

Let’s look at the precise solution.

The code is as follows:

 func solveExact(amplitude: Double, phase: Double, kSpring: Double, mass: Double, t: Double) {
 var x = 0.0  
 //1  
  let omega = sqrt(kSpring / mass) 
   
 var i = 0.0   
   
 while i < 100.0 {
 //2  
    x = amplitude * sin(omega * i + phase)
 i += t
 }
 }

This method contains all the parameters needed to solve the equations of motion.

1. Calculate the resonant frequency

2. Calculate the current position of the object in the while loop, and use i as the self-incrementing variable for the next calculation.

Add the following test code:

 let osci = HarmonicOscillator(kSpring: 0.5 , mass: 10 , phase: 10 , amplitude: 50 , deltaT: 0.1 )
 osci.solveUsingSolver(solveExact)

The method in this solution is a bit strange: it takes parameters but returns no data and does not display anything.

#p#

What are the benefits of doing this?

The purpose of using this function is to calculate the specific dynamic values ​​during the vibration process in the while loop. In the playground, you can observe the history of these dynamic values.

By adding a value history at x = amplitude * sin(omega * i + phase) we can see the trajectory of the motion.

[[142803]]

Now that the first exact algorithm has been verified, let's take a look at the numerical calculation scheme.

Euler method

Euler's method is a method for numerical integration. The method was first proposed in 1768 by Leonhard Euler in his book Institutiones Calculi Integralis. For more information on this method, please refer to: http://en.wikipedia.org/wiki/Euler_method

The core idea of ​​Euler's method is to approximate a curve by using a broken line.

The specific method is to calculate the slope of a given point, and then draw a polyline with the same slope. At the end of this polyline, continue to calculate the slope and draw another line. As you can see, the accuracy of this algorithm depends on the length of the polyline.

Do you want to know what deltaT does?

The numerical algorithm uses a step size, which is important for the accuracy of the algorithm. A large step size results in low accuracy, but faster execution. Conversely, accuracy will increase, but speed will decrease.

The deltaT variable represents the step size. We initialize this value to 0.1, which means we calculate the position of the object every 0.1 seconds. In the Euler algorithm, this means that the length of the projection of the polyline on the x-axis is 0.1.

Before writing any code, you need to look at this formula again:

The second-order differential equation can be transformed into two first-order differential equations. They can be written as and

We can use the difference quotient to complete the conversion:

You will also get:

With these equations in hand, we can implement Euler's method directly.

Add the following code after the solveExact method:

 func solveEuler(amplitude: Double, phase: Double, kSpring: Double, mass: Double, t: Double) {
 //1  
  var x = amplitude * sin(phase)
  let omega = sqrt(kSpring / mass)
 var i = 0.0  
 //2  
  var v = amplitude * omega * cos(phase);
 var vold = v
 var xoldEuler = x 
   
 while i < 100 {
 //3  
    v -= omega * omega * x * t
 //4  
 x += vold * t
 xoldEuler = x 
   
 vold = v
 i+=t
 }
 }

Code content:

1. Set the current position and omega value.

2. Calculate the current speed.

3. In the while loop, calculate the new speed through the first-order function.

4. Calculate the new position using the velocity, and at the end, increase the value of i by the step size t.

Now, add the following test code:

 osci.solveUsingSolver(solveEuler)

Add the value history at xoldEuler = x and check it out. You will see that this method shows a sine function with increasing amplitude. Because the Euler method is not an exact algorithm, the large step size of 0.1 here also leads to inaccurate calculation results.

[[142811]]

Here is another function graph with a step size of 0.01, which is clearly better. So, when using Euler's method, you have to keep in mind that the smaller the step size, the better the result. But using a compromise step size is easier to implement.

[[142812]]

Velocity Verlet

The last algorithm discussed is called Velocity Verlet. It works on the same idea as Euler's method, but the way it calculates the new position is slightly different.

Euler ignores the actual acceleration when calculating the new position, and the formula is:, while the velocity Verlet algorithm takes acceleration into account when calculating:. This gives better results when the step lengths are equal.

Add new code after the solveEuler method:

 func solveVerlet(amplitude: Double, phase: Double, kSpring: Double, mass: Double, t: Double) {
 //1  
  var x = amplitude * sin(phase)
 var xoldVerlet = x
  let omega = sqrt(kSpring / mass) 
   
  var v = amplitude * omega * cos(phase)
 var vold = v
 var i = 0.0   
   
 while i < 100 {
 //2  
    x = xoldVerlet + v * t + 0.5 * omega * omega * t * t
    v -= omega * omega * x * t
 xoldVerlet = x
 i+=t
 }
 }

Code content:

1. All the code before the loop is the same as in Euler's method.

2. Calculate the new position based on the velocity and increase the value of i.

Test the code in a Playground:

 osci.solveUsingSolver(solveVerlet)

[[142815]]

What to do next?

You can download the completed project from this link.

I hope you have fun on this journey through the world of numbers. It would be helpful to know a little about algorithms, or even just some interesting history of classical mathematics.