Simulation of Wave Equation, introduction and discovery

Simple simulator, run, click on the screen and see waves propagate!

Hi! This project was intended to learn Bevy and try by myself to simulate the wave equation. I decided to change a little the document, and instead of having it just for myself, I'll write what I have learned from Bevy, computers and math on this small project for other people.

The wave equation is a math function. I won't go deep explaining each term; instead I'll try to give you the intuition behind each piece, so you can get what I want to say instead of just filling with technicisms.

I'll start with the coding side, then move to the math behind the wave equation, and finally connect the math with computing and the challenges I have found.

Why not just read books and apply them? Because the idea is to learn, to try to be creative, to look for new solutions. Innovation comes when we have limited resources; in this case I'm limiting the knowledge. The solutions you reach now don't need to be the same as the book's, but it will be hard to get new ones if you already know the book's answers.

I took my time playing with this, because I did it in my free time, while at work we need to adapt to the time limit.

There will be a lot of math, I strongly suggest reading it, because all the post-computer parts are about how to integrate the math into a computer, which is a very big topic by itself.

Coding and computer

First: if you've seen the code, it's not perfect! I didn't intend to have code that uses best practices; it would take more time, and the priority was learning Bevy and simulations and their challenges, so I stayed focused on that.

The code is trivial, and I used Bevy. You may think it's overkill for this, but the reason is simple: learning more UI frameworks would require a lot of time. Bevy works here, works in other things, is flexible, organized, and will cover 99.9% of my use cases! So I prefer to keep and learn one well.

To fully understand this section I suggest reading the math section further down, but you can also follow along at a high level. Let's start:

The Space type handles different constants that we need to represent the simulation. The minimum values are the width and height. We need f32 because every point on the Bevy screen uses it, and usize because we need to access an array. Keeping them here helps me keep them consistent, though I know the f32 types should not be here in the physics representation but in a higher layer.

#[derive(Debug, Clone, Copy)]
pub struct Space {
        pub width: f32,
        pub width_usize: usize,
        pub height: f32,
        pub height_usize: usize,
        pub speed: [f32; 2],
}

The Frame represents the actual state of something, maybe speed, position or another thing. Since Space is also cheap, I just keep a copy of it here; I don't change it after it's set, so there's no problem. Also, because each frame should keep its own space properties, they shouldn't belong to a higher layer.

I'm storing the information as a Vec just because Bevy stores the image data this way.

#[derive(Debug, Clone)]
pub struct Frame {
        pub data: Vec<f32>,
        space: Space,
}

The Data struct is the one that keeps all the information needed for the simulation to work. As I explain later in the math section, we need at least the last frame and the previous one to calculate the next one. We also store the space information and some values of the last frame, like the min/max values, to play with the normalization.

I'm recording the currently needed values and the next ones just to avoid allocating everything each time, so I can always use the same memory space to store and perform the calculations for the next frame.

The struct is a Resource; it's the one that is accessed from Bevy to work on.

#[derive(Resource)]
pub struct Data {
        pub velocity: Frame,
        pub position: Frame,
        pub next_velocity: Frame,
        pub next_position: Frame,
        pub max_now: f32,
        pub min_now: f32,
        pub space: Space,
        pub img_cache: Vec<u8>,
}

Challenges and annoying things

There have been some sections that were a headache during this coding. While ECS makes everything very easy, some sections are... let's see them!

Some of these are not really challenges; some are just things I found annoying to handle.

Duplicated numbers

When we want to use grids and values, we have 3 types of them to represent the same concept at different levels, and the shared concept is the value of height and width.

• usize: to access RAM
• f32: to interact with the screen, Bevy uses f32 to represent it
• simulation type: the internal type used to perform the calculation, it can be f64 or f32, depends on your needs

I'm not a fan of having so many types to represent the same thing! And you must also keep the same values consistent. It's not hard, but I don't like the idea of converting values all the time, so in most cases I prefer to have them already converted somewhere to access them directly.

Still, having to centralize one value at different levels of logic is annoying; ideally all of them would be independent.

These values will be used in the simulator, in Bevy to interact, and in Bevy to show the results.

Bandwidth

ECS is cool! Bandwidth is one of the big reasons to pick ECS, especially as a communication protocol.

I did all the code on the CPU. But let's think first; let's take as an assumption a 1920x1920 grid, and we'll use f64 to perform the calculations.

$$ 1920_1920_8\text{ bytes} = 29{,}491{,}200\text{ bytes} \approx 29.5\text{ MB} $$

Each frame will contain 29.5 MB. Sounds small? Well, if we want to achieve 60 FPS, you'll get 1.77 GB/s!! This is not a trivial bandwidth to work with, especially because this is not even our simulation bandwidth; it would be just the bandwidth needed to retrieve the amount of data to process to show 60 FPS!

Each simulation has two sections:

• Simulation: this code is executed most of the time to perform the simulation; this is not shown, it's just the internal code to get the best precision based on our preferences
• Show results: at some point we need to pick the current value calculated on the simulation and show it to the user!

So, if we want the same FPS for simulation and showing results, the bandwidth just doubles, and even more if you want more precision from the simulation. The bandwidth will just keep increasing a lot.

When we perform tracing, actually one of the most expensive operations is converting the actual simulation into Bevy. This is not trivial; let's see the workflow:

1. Perform the simulation calculation n times (at least 1.77 GB/s)
2. Extract a calculation at a specific FPS rate; while we do this, the simulation stops (29.5 MB/s + stop simulation)
3. Convert the simulation type into a Bevy type, prepare the results and show them (at least 1.77 GB/s)

I know, I know, the CPU was never meant to be used like this! But let's solve each thing. There are two main challenges: speed up everything and avoid stopping the simulation when we want to see it.

Bandwidth for speed

OK, this is trivial! We just need to use the GPU, and for that we can use WGPU, a comfy framework that can be executed almost everywhere.

Here is a very interesting thing about ECS, and why it is used as a communication protocol.

Unless we want to write everything in WGPU, you must talk between CPU <-> GPU; each one has its own memory, and if we follow the code I wrote, you should be moving a lot of information between them, and again the bandwidth would just collapse. So here is the interesting part:

1. Bevy makes a UI thing, anything: it's a component loaded into VRAM
2. If you want to interact and perform operations, you make a new component that will be attached to your entity
3. When you want to perform an operation or do something, you use the small component as a communication device to tell the GPU to run anything

A very practical case is the Transform component:

1. Load an image into the GPU
2. Attach a Transform component to the image
3. When you want to scale or rotate, you only touch the Transform component and it will update the values on the GPU
4. The GPU will run its own fast code according to the instruction

Bandwidth for extracting value from simulation

If we want to read information, we need nothing using our current memory section. So if we only have one memory to perform the simulation, there is no way to avoid this!

So, there are several ways. I would choose something like a section of the VRAM with a copy, where we can access and read it when we want.

The system would be simple: a component where:

1. The simulation has its own VRAM
2. A VRAM section with a frame to be extracted
3. A flag to know whether the frame was read
4. Each time we request to read from the VRAM, we get the frame and mark whether it was read
5. The VRAM check, every X ms, looks at whether the frame was read: if it was not read, increase the check time; if it was read, update the value

In the end there will be a balance point, where the GPU will only copy to the new location at the rhythm at which the CPU can handle it.

Logic isolation

While writing the code has been easy enough with Bevy, I'm worried about the code isolation. A simulator and Bevy don't belong to the same logic territory; performing calculations and interactions should at some point have its own API, an interface between them.

But all this is handled by the design pattern, ECS in this case. While some functions can be general for any pattern, we still need to make a way to implement it...

I haven't dug into whether there is an abstract crate for simulation and ECS that we could mix for Bevy, but there are also mixed components:

• The simulation
• Connection to the GPU/WGPU
• Custom components for optimization
• Handling showing the results

And some more, so even if it seems ideal to achieve, this is harder than anticipated.

Color

Showing the results is very hard!

As we talked about, simulations can run on very high and low values, and in real time we have no way to properly know the upper and lower bounds to know how the values will change.

Some simulations will work on small changes, others will work on high scales!

Color on RGB just works on 3 × 255 values; if we want an intuitive way, we would mostly rely on one unique color, which leads us to a 255 color range, which is... far from enough to represent f32 and f64 when we don't know in which range of values the simulation will be working.

Our only option for a real time simulation is to choose a range, or to analyze all the data to know how it changes to be able to choose what and how to show in the results. But this has a big downside: the computational cost is too high; still each frame is close to 30 MB, so doing this from time to time is the best.

I tried normalizing the functions per frame; while we can see clearly how everything changes in the video, with time the intuition of what each color represents makes it really hard to know what each value is. Also, normalization makes it so that a change in the color does not guarantee that the change represents the pixel going up or down, because the value is in relation to each frame independently.

I think with just 255 values, or a very limited range, there will be no good or global solution for this, and the implementation will change with each use case.

For the current code, where I want a perfect physics simulator, I have the following options:

• Sidebar showing the current height values and their colors
• We can use any static range and use zoom; the zoom, instead of increasing resolution, will open a new window that increases the color resolution, mapping the selected sector into a new color scale
• We can make clusters of pixels and create color breaks with two values: one will be the average of the cluster, and the other will represent the internal variation. This would help to know where it would be useful to zoom

Introduction to the wave equation

Note: I know every function will depend on time. Just to avoid filling that variable and make it easier to read, I won't write it unless necessary.

Waves, essentially sound, or just dropping a rock on a pool of water, will give you waves.

We will start with the fundamental definition, which is: you have particles, and you join each particle with springs.

In a string, you have any particle $j$, then you connect a spring with "j - 1" and another one at $j + 1$, so the particle $j$ will react to the forces from its sides.

In a plane you will just connect a particle to each side; you will use 4 springs, 2 for the $x$ axis and 2 for the $y$ axis.

There is a core simplification in the wave equation: the particles do not move from their position! They can only move up and down, but will never move from there.

In all cases you have the height of each particle called $u$, and if it is a string or a grid, you connect them with springs on its sides.

What we just read can be expressed like this:

$$ Force_j = Force_{\text{j-1 to j}} + Force_{\text{j+1 to j}} $$

Spring forces are described by Hooke's law; the force of each spring is described by two values:

$$ Force_{spring} = -k*(x - x_0) $$

A normal spring has a normal length, where the spring exerts no force; this is called $x_0$. When you compress it, it pushes away; if you stretch it, it pulls inwards.

While $k$ is how hard the spring is: higher values make it harder to compress, and when released it causes a lot of force.

The springs in the wave equation have length 0, which means there is always tension between the particles pulling them closer. If it weren't for the assumption that the particles do not move to the sides, our string would compact.

The final equation is:

$$ Force_j = k*(u_{\text{j + 1}} - u_j) + k*(u_{\text{j - 1}} - u_j) $$

This is the fully discrete version. The continuous one is when you perform some tricks: instead of thinking of a grid, $u$ will represent the height of the string at position $x$ of itself, so it is like $u(x)$. Now it is continuous, so when we look to the right and left of the string we are moving $\Delta x$ forwards and backwards.

I won't go deep on this, but basically when you change $u$ in this way, some terms will change and appear:

$u(x)$: height at position $x$

$K$: Global constant of $k$ (depends on it)

$M$: Total mass of the string

$\Delta x$: Particle we will pick forward or backwards to approximate the current force

The next step is crucial. When we mix everything up, there will be a lot of terms like:

$$ \frac{u(x + \Delta x) - u(x)}{\Delta x} $$

When we want to get the continuous version, we make $\Delta x$ go to zero. This is the definition of a limit, which in more intuitive terms is the definition of speed.

Think for a moment: you are walking, you are at position $10m$, you traveled until position $50m$, and you know it takes $5s$. The speed is:

$$ \frac{50m - 10m}{5s} = \frac{40m}{5s} = 8\frac{m}{s} $$

The trick is that if you look at a very close moment of your walk, a very small traveled distance, you get the instant speed you had.

Now think of a similar case: you want to calculate the slope of a mountain. You are at position 5m with a height of 20m, and you know that at position 6m the mountain has a height of 30m. If we want the slope, it is:

$$ \frac{30m - 20m}{6m - 5m} = \frac{10m}{1m} = 10 $$

As you see, calculating the slope is exactly the same as calculating the speed; this is because both are the same concept: we are calculating how a variable changes in relation to another one.

• Speed: change of position in relation to time, $\frac{m}{s}$
• Slope: change of height in relation to position, a ratio

When each term is computed at the closest distance, the closer one position is to the other, or one time to another, this is called the instant change.

These terms of instant change appear a lot in the discrete version of the wave equation, and this has its nomenclature:

$$\frac{\partial u}{\partial t}(t) = \lim_{h_t \to 0} \frac{u(t+h_t) - u(t)}{h_t}$$

Don't be afraid of these terms; this is just speed, but in a math way.

$h$ replaces $\Delta t$. Usually if we are going to do approximations we use $\Delta t$ because it is not zero; when we want to see infinitely close between two moments, we use $h$.

Actually there is also the speed of the speed! $\frac{\partial^2 u}{\partial x^2}(x)$.

Sometimes, we usually say:

$$\frac{\partial u}{\partial t}(t) \approx \frac{u(t+\Delta t) - u(t)}{\Delta t}$$

Which is fine, sometimes.

When we apply this limit to the discrete wave function, we get this:

$$\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}$$

On the left side we have the speed of the speed, the acceleration of how the height changes; while on the right, we have the speed of the speed, but not in relation to time. This time it is the height of the particle in relation to its neighbors, the next one on the $x$ axis.

The new term $c$ is the speed of propagation, which is:

$$c = L \sqrt{\frac{K}{M}}$$

3D Wave equation

Above we have the 2D wave equation. The 3D one is for cases when we have a grid; there are particles to the right, left, above and down. I won't go over everything again; instead I'll post the same logic above but for this case, so we can have a clear nomenclature.

All the force equations follow one rule; the interpretation and input change a little:

$$ Force = \text{Left} + \text{Right} + \text{Above} + \text{Down} $$

Discrete function for a grid, $u_\text{(i, j)}$:

Force Above: $k*(u_{\text{(i, j + 1)}} - u_\text{(i, j)})$

Force Down: $k*(u_{\text{(i, j - 1)}} - u_\text{(i, j)})$

Force Right: $k*(u_{\text{(i + 1, j)}} - u_\text{(i, j)})$

Force Left: $k*(u_{\text{(i - 1, j)}} - u_\text{(i, j)})$

The discrete function for position, which is used to then take the distance between particles to zero:

Force Above: $k*(u(x, y + \Delta h) - u(x, y))$

Force Down: $k*(u(x, y - \Delta h) - u(x, y))$

Force Right: $k*(u(x + \Delta h, y) - u(x, y))$

Force Left: $k*(u(x - \Delta h, y) - u(x, y))$

When we perform the limit, the function that describes the continuous case is:

$$\frac{\partial^2 u}{\partial t^2} = c^2 (\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2})$$

Computer results

I'll write a later section about how I structured the code; I prefer to talk more about how to try to write the previous equations in a computer. You will see why this is key. Still, I tried to optimize almost everything I could find in the code, so we can take this seriously. The main feature I did not integrate, which will be explained later, is WGPU, because it would require me to learn a new language.

The first point here is how to simulate the wave equation, but we need to face the reality about it:

• Discrete version: limited number of particles
• Continuous version: unlimited number of particles

In continuous functions, like a string, where we can know the position as $f(x)$, we have the issue that for any range, even the smallest one of the string, there will always be an infinite number of numbers in it. Like between 0 and 1, we have 0.5; between 0.5 and 1 we have 0.75; you can always continue this infinitely. You might think that with enough particles we could simulate or approximate the continuous formula, but this is one of the major issues: no matter how many particles we put in the simulation, it will never reach infinity. Infinity is underestimated a lot here; any number is infinitely far from infinity!

$$(x - \infty)^2 = \infty$$

This takes a very important place in computers, because computers have limited resources, limited resolution; it is the nature of the digital world. Math, on the other hand, works on a perfect axiomatic world with infinite resolution! We will see how this infinity affects each attempt to approximate the formulas.

Due to this, instead of trying to simulate infinite particles, my first choice was to use the discrete formula. We can use a grid, simulate the springs between them, and try this simple scenario. Still, it did not give the best results.

Discrete function approximation

In order to use the discrete function, we still need to touch a formula. Maybe the number of particles is finite, but the acceleration still works on an infinite $\frac{\partial^2 u}{\partial t^2}$.

We can be brute and just try:

$$ \frac{\partial u}{\partial t}(t) \approx \Delta u(t) = \frac{u(t + h) - u(t)}{h} $$

Where $h$ in this case is the time resolution we want to use. $h$ is important because it is usually used as a value that we know we want to limit to zero; so if we want to approximate slopes it would be the spatial definition, and if we have a grid you have two $h$'s, $h_x$ the distance between particles in the $x$ axis and $h_y$ the distance between particles in the $h_y$ axis.

$$ \frac{\partial^2 u}{\partial t^2} = k*(u_{\text{j + 1}} - u_j) + k*(u_{\text{j - 1}} - u_j) $$

We know all the heights for all the particles in the present in our grid, so we can know the right part of the equation; while the left one can be approximated with the equation above, so we can solve $u(t + h)$ value for all the current particles.

The $\frac{\partial^2 u}{\partial t^2}$ approximation will require values from now and from the previous state to compute the future, so we need at least two copies. Why two? Because only for velocity you have $u(t + h)$, and if you insert that function into itself it is like $u(t + 2h)$, which depends on both.

Result of discrete function

Now we have all the pieces to write some code! I'll write another section about the code; here it will be mainly about what happens.

In a pool, if we drop a rock we get circular waves, but with the current method and code, we will get... squares! There are several reasons why we failed on this.

Information compression

When we have the discrete wave equation you can notice something: we only put springs on the particles on the sides, so at some point the force depends just on:

$$ Force = \text{Particle on the left} + \text{Particle on the right} + \text{Particle down} + \text{Particle up} $$

So if we just connect each particle to those ones, it is to be expected to get squares... then why does the continuous equation give us circles?

It's because in math and physics performing it like this is a known trick: you usually only need what will happen closely, but not perfectly. Then this is the crucial operation:

$$ \text{Continuous function} = \lim_{h \to 0} \text{Discrete Function}(h) $$

The distance between all particles goes to zero! While in a discrete world, we have particles, and each one is not connected to every other; we have directions: particles to the left, particles up, particles on the diagonals, all of them exist, and each relation between them is independent.

In a continuous world instead, you know that if you throw a ball at 1 m/s on the y axis and 2 m/s on the x axis, like a vector, only by looking at the axes you have all the information to know how that point will change.

In a continuous world, each point has its own position, velocity and acceleration, and only needs to be defined at each axis, because we are looking infinitely small, so small that each particle is surrounded by an infinite number of particles.

While in the discrete world there is a distance between each particle and its diagonals which we did not connect with springs; when we make all distances zero, all diagonals collapse onto the point, all information of forces is compressed lossily into just the axes of the grid, so each force in each direction can be written using only those two forces. In discrete they just don't exist.

This effect of looking at the infinity causes that a rock that would have a square wave in a discrete world has a circular one in the continuous function.

There are still some tricks: we can increase the resolution. Increasing it a lot, and consuming more time and computation, I was able to get something like a circular shape, but only in the exterior of the wave. The interior, no matter how much I increased the resolution or changed the values, the inner part of the wave was always a square.

Depending on the parameters we can get some weird and curious shapes; this is not a bad result, but it changes how we can work on this. The main change would be that the question would be:

Which combination of parameters in the discrete functions is closest to the continuous one's shape?

Stability

I had... a hard time with this concept. Let's think about it: if we drop a rock in a pool, will the wave be the same? Will we be able to see the wave all the time?

This point makes the difference between "seems to work" and "maybe works"; this will also be the difference between "need" and "enough".

Each wave is made of particles. Each particle has a height and a velocity. If the wave does not have any reason to go faster, just as an imaginary example, the total energy is:

$$ \text{Total Energy} = \sum_{i}{\frac{h_i^2}{2}} $$

We know it's a circle, so as an approximation let's get the perimeter of it:

$$ P = \pi * r $$

$$ \text{Total Energy} = P*\sum_{i}{\frac{h_i^2}{2}} $$

Now, if we look at the wave when it expands... the energy increases! Obviously, more time passes, the circle is bigger, the energy is bigger.

Obviously we know the energy should not increase, that's not the exercise. The point is more, when we look at an image of the simulation: if we always see the wave, we know it's not working!

Sometimes we see our simulations; while we see the wave we can think "it's working", and if it disappears we think "it's not working!".

Computers, heuristics, approximations have the big downside that they do not guarantee conservation. They don't guarantee that the energy will increase nor that it will go down, nor whether it'll be one or the other! Our methods could just increase the energy in some parts and reduce it in others.

So what we see in a simulation can be wrong, or can be right: if we see something it might not be working; if we see nothing it might be working. The first step is to recognize that what we see is not proof of whether the result is working or not.

Be careful: it is hard to know whether any value or metric stability means the simulation works right.

Energy Stability

If we have methods where there is no guarantee for energy stability, then we can trick it, at least to guarantee the right total of energy. This really helps with solving methods, computer approximations, and numeric robustness.

Not being able to have stable energy is... a disaster. If the energy goes up, all our values will go NaN due to overflow to positive or negative; if the energy goes down, we will always end up with static objects, a frozen system.

Instead of trying to make each method perfect, let's try to just have the energy be stable on each frame of the simulation, to minimize the impact of the rest of the system.

First, when a simulation starts, we need the total of the energy. Then on each user interaction we need to record how the total energy changed. The simulator allows clicking, so on each click we check: are we changing the speed? the height? And on each interaction we calculate again the total energy of the system. This is very efficient because only matters if something externally increases or decreases the energy; internally it will remain constant.

With the $\text{Total Energy}$ value, we perform all the methods and functions to get the next frame, then we can scale it to conserve the energy the most we can.

We know the total energy of the system is:

$$ \text{Total Energy} = \sum_{i}{\frac{h_i^2 + v_i^2}{2}} $$

With $h$ as the height of the particle and $v$ as the velocity. They don't need to be the same on each frame, so after getting the new frame we can use that formula to get the new one. Let's say height and velocity are scaled by $k$.

$$ \text{Total Energy Post} = \sum_{i}{\frac{(h_i_k)^2 + (v_i_k)^2}{2}} = k^2*\sum_{i}{\frac{h_i^2 + v_i^2}{2}} $$

So:

$$ \text{Total Energy Post} = k^2*\sum_{i}{\frac{h_i^2 + v_i^2}{2}} = k^2*\text{Total Energy} $$

With this we can get a $k$ value that will make the system always be closer to the right energy it should have.

$$ k^2 = \frac{\text{Total Energy Post}}{\text{Total Energy}} $$

Energy, approximations and infinity

I have pointed out methods, ways to trick and approximate a wave equation even at basic levels, but there is one crucial point I have not talked about that makes a lot of approximations a hard topic: the derivative definition.

As we saw above, velocity is:

$$ f'(x) = \lim_{h \to 0}{\frac{f(x + h) - f(x)}{h}} $$

The tricky point, more than understanding the definition, is the convergence of two things.

$f(x + h) - f(x)$ will always be zero, and at some point computable even as an approximation.

$\frac{1}{h}$ will always be infinity, no way to be computable. Remember, Inf from Floats is not infinite; it is just an overflow, but it is a really big number! Like 1.7976931348623157 × 10³⁰⁸!!

Right, F64 can have very big numbers, but how far is that from a true infinity?

$$ \frac{1.7976931348623157*10^308}{\infty} = 0 $$

No matter how big a real number is, as long as the value itself is not an infinity type (there are several types of infinity), there is no way a computer will be able to truly get the values.

Maybe you may think that they still work, at some point will show us the values, and that's right, but this only works in one way.

Throw a ball up, it will make a parabola, it will go up and down. We can sample the curve, calculate the approximation of the values, but if we pick these values and put them on the physics machine they will stop working, because each step of the ball going up and down is meant to be executed an infinite number of times.

This matters a lot in the wave equation. After all, all the methods must handle $\frac{\partial^2 u}{\partial t^2}$; time is continuous, we can't expect using approximations like the ones I used here to work and fix this infinity.

This is not a problem, it's the nature of math being perfect on its axioms and computers working on limited resources. Accepting infinity as a gap between them allows us to find better solutions.

F64 and Math

When it comes to resolution, we can say Floats and in particular F64 are great; they have a lot of digits and a clever design to adapt to several magnitudes. And still they don't work very well when it comes to math.

The main issue, which I have been repeating a lot throughout the document:

• Math: perfect, infinite resolution, perfect operations, everything works over its axioms, it can even work on different infinities!
• Computers: the art of limits, limited resolution, limited space, limited computing power, limited energy, limited everything

The reality about this is that most things in each one of those worlds are not portable to the other. A math equation can hardly be used correctly on computers, and values computed on a computer can make no sense in math because they usually come with approximation. So you can't validate computer results. Want to use an error range? Maybe some epsilon to say the values are right? In math there is a saying, "prove it!"; if you can't, it just seems to work.

Physics and simulations follow very strict rules, like the one we saw above about energy conservation; it is not minor, it cannot be changed. Any approximation, any value out of place, and all guarantees will be lost and the effects can change, as we talked about with square waves and weird things.

There is a whole topic in math, to solve problems numerically, because of this gap. There is a lot of work! And we still must be cautious, because a lot of this work is done by math people, which means it uses math, and again math doesn't work right on computers!

Example: I was looking a little, and I found this, it's called the CFL condition. If you follow it, in theory the wave equation should at least not explode with NaN:

$$ \frac{c*\Delta t}{\Delta x} <= 1 $$

where $c$ is the speed of the wave, $\Delta t$ our time between frames, $\Delta x$ the distance between particles.

Now let's think a little. Floats... Floats are not numbers, there are no numbers. Let's remember its definition: Floats are approximations. That's why $\frac{0}{0}$ is $NaN$, because you can't guarantee that the $0$ is close to zero from the positive or negative side. So how can we really get the values?

Before, you would say $1/0$ is $\infty$ on Floats; let's remember that they have a positive zero and a negative zero. They just mean "numbers very close to zero from the negative or positive side", but none of them really means zero!

The right interpretation is: Floats are ranges. You don't know the next digit value, at least if we pick a good source, so we need to get the max value of $\frac{c*\Delta t}{\Delta x}$ to know if they are lower. Since each value is a range we would have:

$$ \frac{c_\text{max}*\Delta t_\text{max}}{\Delta x_\text{min}} <= 1 $$

We're right... right? Well, at least we were able to use that formula properly, but there is another issue.

This one comes from following a method for the wave equation. The method itself is not important; the right question is whether that method is compatible with Floats, whether applying the result in computers will really give the result that was proposed using math. This is not trivial.

Each algorithm, each method, each theorem follows its own rule. Simulations are very sensitive to any perturbation at any scale if we want a good one. Also remember, computers accumulate errors; this means any perturbation at any scale will propagate over time and space ($t$, $x$, $y$ axes).

There is no easy way to handle this, but there are a lot more methods to handle each scenario. This is more illustrative, to show the challenges.

Simulation type

We've traveled through a lot of things; there have been a lot of conclusions, but there is one that will define our path: which type of simulation do you need?

You may need:

• Approximation even if it does not converge
• Values or ranges
• Real physics simulator
• Consistency on the physics results

There may be more options, and when you look for methods to approximate the wave equation, you also should read the instructions: usually each method can give you a different result, different guarantees, different output. So if you look at what exists, look at what you need and which method suits you better. If there is none, you must take a tradeoff: more imperfection somewhere, or time to develop something new.

I want to give a special mention to the word "consistency". We don't always want a real physics simulator for our world; sometimes you only need the world to be consistent, like with energy preservation, no weird behavior, all right. This opens the door to a very interesting option: discrete physics.

A lot of physics works on a continuous world, like $\frac{\partial{u}}{\partial{t}}$ where time is perfectly continuous. But in physics most things follow a rule: minimize the energy. While in continuous it is to optimize an equation, this can be moved to discrete too. There the word action is used.

We can, and these things exist in physics, take all these continuous formulas and rules into a discrete version that minimizes action: for springs, for particles, for everything. When everything becomes discrete, a lot of the current challenges disappear! Only handling F64 resolution remains, but this method really helps.

Although it has its downsides, it will be consistent in the computer, but it will not represent our world or scale. Now $\frac{\partial{u}}{\partial{t}}$ will have a different value, but at least it will be computable. The result will differ, but while the previous approximation broke the rules, this one will keep things right. It is perfect if you are looking for your own internal simulator rather than a reality simulator; a videogame is a good example where you don't need to follow our reality's rules but need consistency.