PK iA7Y Y L Unit_7_-_Advanced_Applications_of_Numerical_Methods/01-Welcome_Unit_7.en.srt1
00:00:00,000 --> 00:00:02,000
Welcome back.
2
00:00:02,000 --> 00:00:04,000
By now you have seen all kinds of models based on differential equations
3
00:00:04,000 --> 00:00:08,000
and all kinds of tricks to solve them on the computer.
4
00:00:08,000 --> 00:00:10,000
If that was all there is to it,
5
00:00:10,000 --> 00:00:12,000
thousands of phD students would be jobless.
6
00:00:12,000 --> 00:00:15,000
In this final unit, we'll touch upon more advanced topics
7
00:00:15,000 --> 00:00:18,000
such as modeling buildings and air and water.
8
00:00:18,000 --> 00:00:21,000
We'll explore why it's hard to predict the weather,
9
00:00:21,000 --> 00:00:24,000
and we'll look into aspects of software engineering,
10
00:00:24,000 --> 99:59:59,000
that is, how to do computer simulations in a proficient way.
PK iAJ1 S Unit_7_-_Advanced_Applications_of_Numerical_Methods/02-Finite_Element_Method.en.srt1
00:00:00,000 --> 00:00:05,000
Let's start with a bird's eye view on the finite-element method--FEM.
2
00:00:05,000 --> 00:00:10,000
These days it's the workhorse of almost all mechanical engineers.
3
00:00:10,000 --> 00:00:16,000
FEM can answer, for instance, what happens when you put a huge truck on a bridge.
4
00:00:16,000 --> 00:00:21,000
To compute this deformation is application of FEM to the static case,
5
00:00:21,000 --> 00:00:24,000
but FEM can also be applied in a dynamic setting,
6
00:00:24,000 --> 00:00:28,000
for instance, to simulate the effect of a crash on a car body.
7
00:00:28,000 --> 00:00:31,000
Here it's important to look at the process of buckling.
8
00:00:31,000 --> 00:00:37,000
Whereas in the static case we're not really interested in how the truck was placed on the bridge.
9
00:00:37,000 --> 00:00:39,000
It just has to be there.
10
00:00:39,000 --> 00:00:43,000
I want to outline three fundamental ideas of the finite element method.
11
00:00:43,000 --> 00:00:46,000
The first is discretization.
12
00:00:46,000 --> 00:00:53,000
The continuous structures are approximated with the help of--guess what--finite elements--
13
00:00:53,000 --> 00:00:59,000
elements of finite size, not infinitesimal size.
14
00:00:59,000 --> 00:01:05,000
When we do so the first question is which geometry these finite elements should have.
15
00:01:05,000 --> 00:01:07,000
Should they be tetrahedra?
16
00:01:07,000 --> 00:01:09,000
Should they be cubes?
17
00:01:09,000 --> 00:01:11,000
Or should they even be curvilinear?
18
00:01:11,000 --> 00:01:15,000
The second fundamental idea is that of interpolation.
19
00:01:15,000 --> 00:01:20,000
Given the finite elements, how do I compute a value at an arbitrary location?
20
00:01:20,000 --> 00:01:26,000
For the static case, a really fundamental idea is that of minimization
21
00:01:26,000 --> 00:01:28,000
of the potential energy.
22
00:01:28,000 --> 00:01:32,000
Think about a ball that rolls on a terrain of mountains and valleys.
23
00:01:32,000 --> 00:01:36,000
Eventually, it's going to come to rest in a valley.
24
00:01:36,000 --> 00:01:39,000
In the static case, potential energy is minimized.
25
00:01:39,000 --> 00:01:42,000
Let's have a closer look at that valley.
26
00:01:42,000 --> 00:01:48,000
In the Gedankenexperiment, it splices this object a little further to the left
27
00:01:48,000 --> 00:01:50,000
or a little further to the right.
28
00:01:50,000 --> 00:01:55,000
Then the energy stays almost the same because we're at the bottom of the valley,
29
00:01:55,000 --> 00:02:01,000
which means that to displace this object in this way require no work.
30
00:02:01,000 --> 00:02:03,000
This is the concept of virtual work.
31
00:02:03,000 --> 00:02:07,000
For all infinitesimal displacements that are allowed--
32
00:02:07,000 --> 00:02:11,000
we can't go down and we can't go up, obviously--
33
00:02:11,000 --> 00:02:14,000
the virtual work equals 0.
34
00:02:14,000 --> 00:02:19,000
In mathematics, this way of posing the problem with the help of virtual work
35
00:02:19,000 --> 00:02:21,000
is called a weak form.
36
00:02:21,000 --> 00:02:25,000
The strong form would be to ask for all forces to compensate.
37
00:02:25,000 --> 00:02:33,000
The weak form asks for the virtual work to be equal to 0 for all allowed virtual displacements.
38
00:02:33,000 --> 00:02:39,000
This weak form results in a finite number of equations that we can solve on the computer.
39
00:02:39,000 --> 99:59:59,000
This finite number, however, may range in the hundred thousands or even millions.
PK iA8ob J Unit_7_-_Advanced_Applications_of_Numerical_Methods/03-Hanging_Rope.en.srt1
00:00:00,000 --> 00:00:03,000
Now let's try out some of these ideas with the flexible rope.
2
00:00:03,000 --> 00:00:06,000
We'll fix both ends of that rope
3
00:00:06,000 --> 00:00:11,000
and want to see what the equilibrium shape is going to be
4
00:00:11,000 --> 00:00:13,000
under the influence of gravity.
5
00:00:13,000 --> 00:00:18,000
The obvious choice of finite elements is springs.
6
00:00:18,000 --> 00:00:23,000
Springs of a given rest length, so that the potential energy of each spring
7
00:00:23,000 --> 00:00:29,000
amounts to 1/2 times the spring constant times the square of the extension of the spring,
8
00:00:29,000 --> 00:00:34,000
which is the distance between the two endpoints minus the rest length of the string.
9
00:00:34,000 --> 00:00:40,000
To model the mass of that rope, we attach mass points to these strings.
10
00:00:40,000 --> 00:00:45,000
Of course, these mass points carry potential energy due to gravity.
11
00:00:45,000 --> 00:00:50,000
Given the constants that we provide, compute the potential energy of that rope.
12
00:00:50,000 --> 00:00:55,000
The one end of that rope will be fixed at x = 0, y = 0.
13
00:00:55,000 --> 00:01:05,000
The other end of that rope will be fixed at x = 1.3 m and y = 0.4 m.
14
00:01:05,000 --> 00:01:08,000
Our code starts with a random initialization
15
00:01:08,000 --> 00:01:13,000
and then applies a pretty simplistic strategy to minimize the energy.
16
00:01:13,000 --> 00:01:16,000
For a certain number of times it's going to pick one of the masses
17
00:01:16,000 --> 00:01:20,000
and change the position of that mass point by a random vector.
18
00:01:20,000 --> 00:01:23,000
If the energy decreases, it keeps that new position.
19
00:01:23,000 --> 00:01:26,000
If it doesn't, it returns to the position before.
20
00:01:26,000 --> 99:59:59,000
Very simple, but highly inefficient.
PK iAeR R S Unit_7_-_Advanced_Applications_of_Numerical_Methods/04-Hanging_Rope_Solution.en.srt1
00:00:00,000 --> 00:00:06,000
In the implementation we are summing up the energy for all points,
2
00:00:06,000 --> 00:00:14,000
which is mass times gravitation acceleration times y position for every point
3
00:00:14,000 --> 00:00:16,000
and the energy for every spring.
4
00:00:16,000 --> 00:00:20,000
The resulting diagram may be regarded as abstract art,
5
00:00:20,000 --> 00:00:22,000
but also shows the evolution of that rope.
6
00:00:22,000 --> 00:00:25,000
We are starting with a random configuration,
7
00:00:25,000 --> 00:00:29,000
but eventually the shape of the rope becomes plausible.
8
00:00:29,000 --> 00:00:35,000
Keep in mind that this is not--I repeat--not temporal evolution,
9
00:00:35,000 --> 00:00:37,000
even though it looks like that.
10
00:00:37,000 --> 00:00:42,000
This shows a path of optimization toward an equilibrium.
11
00:00:42,000 --> 00:00:45,000
We are not talking about forces or acceleration here.
12
00:00:45,000 --> 99:59:59,000
It's a search process, like finding the shortest path to connect 10 cities.
PK iAnS S L Unit_7_-_Advanced_Applications_of_Numerical_Methods/05-Fluid_Dynamics.en.srt1
00:00:00,000 --> 00:00:05,000
A topic highly related to the finite element method is computational fluid dynamics--CFD--
2
00:00:05,000 --> 00:00:08,000
the study of gases and liquids with the help of the computer.
3
00:00:08,000 --> 00:00:14,000
For instance, to find the optimum shape of an airfoil or the optimum shape of a car body
4
00:00:14,000 --> 00:00:16,000
or to design hydraulic machinery.
5
00:00:16,000 --> 00:00:19,000
I've titled this section "Why CFD is hard,"
6
00:00:19,000 --> 00:00:21,000
so now let's look into that.
7
00:00:21,000 --> 00:00:26,000
If we were looking at a single particle of mass m that is subject to a force F
8
00:00:26,000 --> 00:00:30,000
Then the rate of change of the velocity would be proportional to that force.
9
00:00:30,000 --> 00:00:34,000
That's one of Newton's laws--force equals mass times acceleration,
10
00:00:34,000 --> 00:00:37,000
which is the rate of change of velocity,
11
00:00:37,000 --> 00:00:40,000
the derivative of velocity with respect to time.
12
00:00:40,000 --> 00:00:43,000
For the fluid, something similar has to happen,
13
00:00:43,000 --> 00:00:46,000
but now we're not dealing with a single particle.
14
00:00:46,000 --> 00:00:49,000
We're dealing with a virtually infinite amount of particles.
15
00:00:49,000 --> 00:00:54,000
What we are working with is not the velocity of the particle.
16
00:00:54,000 --> 00:00:56,000
It's the velocity field.
17
00:00:56,000 --> 00:00:59,000
For every position in space, we specify the velocity,
18
00:00:59,000 --> 00:01:03,000
so the velocity that we specify is the velocity of that particle
19
00:01:03,000 --> 00:01:05,000
at that instant of time.
20
00:01:05,000 --> 00:01:09,000
Before and after, most probably, this location is going to be occupied
21
00:01:09,000 --> 00:01:12,000
by other particles at other times.
22
00:01:12,000 --> 00:01:16,000
When we write down Newton's equation for this particle--
23
00:01:16,000 --> 00:01:20,000
force equals mass times the derivative of velocity with respect to time--
24
00:01:20,000 --> 00:01:22,000
We have to be a little careful.
25
00:01:22,000 --> 00:01:25,000
Let's look at what happens after a very short time step.
26
00:01:25,000 --> 00:01:28,000
It's mass times 1 over the time step,
27
00:01:28,000 --> 00:01:33,000
and now we have to form the difference of the velocity after that time step
28
00:01:33,000 --> 00:01:36,000
minus the velocity before that time step.
29
00:01:36,000 --> 00:01:38,000
The before part is easy.
30
00:01:38,000 --> 00:01:43,000
That's simply our velocity field at the current time and the current location.
31
00:01:43,000 --> 00:01:45,000
The tricky thing is the after part.
32
00:01:45,000 --> 00:01:50,000
It's the velocity field at the later time--t plus time step.
33
00:01:50,000 --> 00:01:54,000
Now we have to take care of the fact that our particle has moved a little.
34
00:01:54,000 --> 00:01:58,000
We don't need the velocity field at that later time.
35
00:01:58,000 --> 00:02:02,000
At position x it has to be a slightly different position,
36
00:02:02,000 --> 00:02:04,000
namely, how far did we advance?
37
00:02:04,000 --> 00:02:08,000
We advanced by time step times velocity.
38
00:02:08,000 --> 00:02:10,000
Now, this is going to make things ugly.
39
00:02:10,000 --> 00:02:15,000
The velocity field of something that includes the velocity field.
40
00:02:15,000 --> 00:02:17,000
A function applied to itself.
41
00:02:17,000 --> 00:02:19,000
This is what makes things ugly and eventually leads to computational fluid dynamics
42
00:02:19,000 --> 00:02:22,000
and eventually leads to computational fluid dynamics being hard.
43
00:02:22,000 --> 00:02:25,000
If we do the math right, this becomes the following.
44
00:02:25,000 --> 00:02:29,000
First we have to look into the change of the velocity field with time,
45
00:02:29,000 --> 00:02:33,000
so we get its partial derivative with respect to time.
46
00:02:33,000 --> 00:02:36,000
But then we also have to look into its change with position,
47
00:02:36,000 --> 00:02:40,000
which is the partial derivative with respect to x, for instance.
48
00:02:40,000 --> 00:02:44,000
The larger the velocity is, the more effect the spatial derivative has.
49
00:02:44,000 --> 00:02:50,000
What we get in the end is the x component of the velocity times the partial derivative
50
00:02:50,000 --> 00:02:52,000
of the velocity with respect to x.
51
00:02:52,000 --> 00:02:55,000
Of course, the same happens with y and z.
52
00:02:55,000 --> 00:02:58,000
This is going to be the acceleration,
53
00:02:58,000 --> 00:03:00,000
and this is inherently nonlinear.
54
00:03:00,000 --> 00:03:05,000
We have a product of a function that we're looking for--the velocity field--with itself.
55
00:03:05,000 --> 00:03:10,000
This is going to make solving the differential equation that results from this really hard.
56
00:03:10,000 --> 00:03:13,000
Finally, however, even though the resulting equation--
57
00:03:13,000 --> 00:03:17,000
the so-called Navier Stokes equation--is going to look pretty complex,
58
00:03:17,000 --> 99:59:59,000
it's nothing else but Newton's law applied to the velocity field.
PK iAQE E Q Unit_7_-_Advanced_Applications_of_Numerical_Methods/06-Force_from_Pressure.en.srt1
00:00:00,000 --> 00:00:05,000
At every point of our fluid we have a velocity, but also we have pressure.
2
00:00:05,000 --> 00:00:08,000
Velocity has direction. Pressure doesn't.
3
00:00:08,000 --> 00:00:12,000
Now let's try to find out how to incorporate pressure in our equations.
4
00:00:12,000 --> 00:00:15,000
To make things simple, we look at a one-dimensional situation--
5
00:00:15,000 --> 00:00:18,000
a straight pipe filled with fluid.
6
00:00:18,000 --> 00:00:21,000
Capital a is the cross-section area of that pipe.
7
00:00:21,000 --> 00:00:26,000
Let's replace some amount of that fluid with the piston of width Δx,
8
00:00:26,000 --> 00:00:30,000
and let's assume that the pressure to the left of that piston is p₁
9
00:00:30,000 --> 00:00:33,000
and the pressure to the right of that piston is p₂.
10
00:00:33,000 --> 00:00:38,000
Then this force is proportional to the difference of those two pressures.
11
00:00:38,000 --> 00:00:43,000
If p₁, the pressure on the left, is larger than p₂, the pressure on the right,
12
00:00:43,000 --> 00:00:45,000
the force is to the right.
13
00:00:45,000 --> 00:00:48,000
So we have to choose sides like that.
14
00:00:48,000 --> 00:00:51,000
The proportionality constant is simply the area.
15
00:00:51,000 --> 00:00:54,000
If you double the area, you double the force.
16
00:00:54,000 --> 00:00:58,000
If we say that this piston has a volume of v,
17
00:00:58,000 --> 00:01:04,000
then the area, of course, is the volume of the piston divided by its width.
18
00:01:04,000 --> 00:01:09,000
If we want to use this as part of our differential equations, what should we be writing?
19
00:01:09,000 --> 00:01:11,000
What's the force per volume?
20
00:01:11,000 --> 00:01:14,000
The partial derivative of pressure with respect to x?
21
00:01:14,000 --> 00:01:18,000
Minus the partial derivative of the pressure with respect to x?
22
00:01:18,000 --> 99:59:59,000
Half the partial derivative or the pressure divided by the square root of the area?
PK iA'' Q Unit_7_-_Advanced_Applications_of_Numerical_Methods/07-Force_from_Pressure.en.srt1
00:00:00,000 --> 00:00:04,000
It's the second one--minus the partial derivative of pressure with respect to x.
2
00:00:04,000 --> 00:00:07,000
If you look closely into what we've got here,
3
00:00:07,000 --> 00:00:10,000
(p₁ - p₂) / Δx--
4
00:00:10,000 --> 00:00:16,000
if this was (p₂ - p₁) / Δx, we would get the derivative,
5
00:00:16,000 --> 99:59:59,000
but it's the negative of it.
PK iA;,j4b
b
O Unit_7_-_Advanced_Applications_of_Numerical_Methods/08-The_Lorenz_System.en.srt1
00:00:00,000 --> 00:00:02,000
Because fluid dynamics is so hard,
2
00:00:02,000 --> 00:00:06,000
many people try to come up with toy models to at least get a glimpse
3
00:00:06,000 --> 00:00:08,000
of what's happening.
4
00:00:08,000 --> 00:00:11,000
One of these toy models--a particularly interesting one--is the Lorenz system.
5
00:00:11,000 --> 00:00:16,000
By the way, this name is not to be confused with that of the physicist who spelled with "tz."
6
00:00:16,000 --> 00:00:19,000
This is what is described by this model.
7
00:00:19,000 --> 00:00:24,000
There is a layer of incompressible liquid between a surface of constant high temperature
8
00:00:24,000 --> 00:00:26,000
and a surface of constant low temperature.
9
00:00:26,000 --> 00:00:28,000
We've got some melting ice-cubes here.
10
00:00:28,000 --> 00:00:34,000
The velocity field in that liquid is described by three highly abstract parameters,
11
00:00:34,000 --> 00:00:36,000
called x, y, and z.
12
00:00:36,000 --> 00:00:39,000
They describe the relative amount of different sorts of motion.
13
00:00:39,000 --> 00:00:42,000
Even though this is a very sketchy and abstract description,
14
00:00:42,000 --> 00:00:45,000
some people have succeeded in implementing these equations
15
00:00:45,000 --> 00:00:47,000
with an actual water wheel.
16
00:00:47,000 --> 00:00:51,000
The idea about studying these equations is that they may be telling us something
17
00:00:51,000 --> 00:00:55,000
about every system in which convection is present.
18
00:00:55,000 --> 00:00:59,000
When we plot the solutions of the Lorenz system in 3D this is what we get--
19
00:00:59,000 --> 00:01:01,000
a very intricate butterfly-like pattern.
20
00:01:01,000 --> 00:01:05,000
Here we're seeing two solutions--one in green and one in blue--
21
00:01:05,000 --> 00:01:09,000
starting from almost the same point--something to explore in the next section.
22
00:01:09,000 --> 00:01:13,000
The trajectories make some turns on one ring, then jump to the other ring,
23
00:01:13,000 --> 00:01:15,000
jump back to the first ring and so on.
24
00:01:15,000 --> 00:01:20,000
The number of turns they take per ring seems to be almost unpredictable.
25
00:01:20,000 --> 00:01:22,000
There is one more thing.
26
00:01:22,000 --> 00:01:25,000
Even though both trajectories have almost the same starting point,
27
00:01:25,000 --> 00:01:29,000
you can see that they quickly diverge in the course of time.
28
00:01:29,000 --> 99:59:59,000
We're going to study that in the next segment.
PK iAi3N N V Unit_7_-_Advanced_Applications_of_Numerical_Methods/09-Lorenz_and_Forward_Euler.en.srt1
00:00:00,000 --> 00:00:04,000
Now we are going to look into the strange behavior concerning the two solutions
2
00:00:04,000 --> 00:00:06,000
that start at almost the same point.
3
00:00:06,000 --> 00:00:09,000
This is about deterministic chaos.
4
00:00:09,000 --> 00:00:12,000
Implement the forward Euler method for the Lorenz system
5
00:00:12,000 --> 00:00:16,000
and do this computation for two initial conditions
6
00:00:16,000 --> 00:00:19,000
that are only slightly--every so slightly--different from each other.
7
00:00:19,000 --> 00:00:24,000
Our program will eventually plot how the distance between these two trajectories evolves.
8
00:00:24,000 --> 00:00:28,000
If you first want to see the butterfly, uncomment the lines below
9
00:00:28,000 --> 99:59:59,000
and comment the lines above.
PK iAO O V Unit_7_-_Advanced_Applications_of_Numerical_Methods/10-Lorenz_and_Forward_Euler.en.srt1
00:00:00,000 --> 00:00:04,000
This is a straightforward implementation of the forward Euler method.
2
00:00:04,000 --> 00:00:09,000
This x at the initial step is actually 2D vector. S is the y.
3
00:00:09,000 --> 00:00:15,000
The result will again be a 2D vector with those two values for the two trajectories in it.
4
00:00:15,000 --> 00:00:18,000
The resulting diagram may need some explanation.
5
00:00:18,000 --> 00:00:24,000
We see the distance between the x values of the two trajectories over the course of time.
6
00:00:24,000 --> 00:00:28,000
Initially, they are really close--10^-14.
7
00:00:28,000 --> 00:00:33,000
And then they grow and eventually become macroscopic--10^100.
8
00:00:33,000 --> 00:00:37,000
The important thing to be seeing is this linear slop.
9
00:00:37,000 --> 00:00:39,000
But actually it's not linear.
10
00:00:39,000 --> 00:00:42,000
We are working with a logarithmic vertical axis,
11
00:00:42,000 --> 00:00:45,000
which means that the distance grows in exponential fashion.
12
00:00:45,000 --> 00:00:53,000
It takes about 30 units of time to increase the error from 10^-14 to 10^-4.
13
00:00:53,000 --> 00:00:59,000
These are 10 orders of magnitude, a factor of 1 to 10 billion.
14
00:00:59,000 --> 00:01:01,000
This leads to the following conclusion:
15
00:01:01,000 --> 00:01:04,000
whenever there is an uncertainty about the initial conditions--
16
00:01:04,000 --> 00:01:10,000
however tiny--the effect of this uncertainty will quickly grow to a visible size.
17
00:01:10,000 --> 00:01:13,000
So, even those the system is deterministic,
18
00:01:13,000 --> 00:01:16,000
once we have fixed the initial conditions, everything is determined--
19
00:01:16,000 --> 00:01:20,000
even though that is the case, it's behavior looks random.
20
00:01:20,000 --> 00:01:25,000
This sensitive dependence on initial conditions is the hallmark of deterministic chaos.
21
00:01:25,000 --> 00:01:27,000
That is to be expected.
22
00:01:27,000 --> 00:01:32,000
Every system, which includes convection, also inherits this aspect.
23
00:01:32,000 --> 00:01:38,000
If somewhere a butterfly flaps its wings, some days or months afterward
24
00:01:38,000 --> 99:59:59,000
this is going to change the weather at some remote location on the earth.
PK iArL`g L Unit_7_-_Advanced_Applications_of_Numerical_Methods/11-Roundoff_Error.en.srt1
00:00:00,000 --> 00:00:04,000
Up to now we have only dealt with the effect of a finite time step,
2
00:00:04,000 --> 00:00:10,000
but actually also the representation of the values on the x axis is cross-grained.
3
00:00:10,000 --> 00:00:12,000
Now we have a look at what this leads to,
4
00:00:12,000 --> 00:00:15,000
and for simplicity it's again the forward Euler method.
5
00:00:15,000 --> 00:00:21,000
As a test case, we're using that satellite that's almost geostationary from Unit 2.
6
00:00:21,000 --> 00:00:25,000
It takes precisely 24 hours to complete one orbit around the earth.
7
00:00:25,000 --> 00:00:29,000
What's critical about the forward Euler method in terms of round-off error
8
00:00:29,000 --> 00:00:32,000
is this addition and this one.
9
00:00:32,000 --> 00:00:35,000
X(t) is a number of reasonable size.
10
00:00:35,000 --> 00:00:37,000
In this case quite an amount of meters.
11
00:00:37,000 --> 00:00:41,000
I'm just drawing some ghost figures here to give an impression,
12
00:00:41,000 --> 00:00:46,000
but the step size times the velocity is pretty small in comparison,
13
00:00:46,000 --> 00:00:50,000
because we're multiplying by that small step size.
14
00:00:50,000 --> 00:00:55,000
If the velocity is of reasonable size and we multiply by that small step size,
15
00:00:55,000 --> 00:00:58,000
we get a number that's rather small in comparison to x.
16
00:00:58,000 --> 00:01:03,000
Then we form the sum of these two numbers and get something that looks like this,
17
00:01:03,000 --> 00:01:07,000
but actually this number again is computed and is being stored
18
00:01:07,000 --> 00:01:13,000
at the position of x--the loose position in that process.
19
00:01:13,000 --> 00:01:17,000
Actually, you see the value of x didn't have enough position to start with.
20
00:01:17,000 --> 00:01:21,000
The smaller we choose the step size, the smaller this number becomes,
21
00:01:21,000 --> 00:01:24,000
the more grave the effect of round-off error will be.
22
00:01:24,000 --> 00:01:31,000
The effect of the round off error is that we actually are losing the details of velocity and force.
23
00:01:31,000 --> 00:01:35,000
We could have gotten away with not computing these digits anyway.
24
00:01:35,000 --> 00:01:40,000
To see a noticeable effect of round-off errors in a manageable time, I'm cheating a little here.
25
00:01:40,000 --> 00:01:47,000
I'm using standard Euler, but them I'm multiplying by a number of 1 plus tiny random number.
26
00:01:47,000 --> 00:01:53,000
This epsilon is defined to be 2 10⁻⁶ in advance.
27
00:01:53,000 --> 00:01:58,000
This looks as though I have a pretty huge round-off error every time.
28
00:01:58,000 --> 00:02:00,000
This is the result that we get.
29
00:02:00,000 --> 00:02:04,000
First I have to explain why we see several points for each step size.
30
00:02:04,000 --> 00:02:07,000
This is because we're dealing with random numbers
31
00:02:07,000 --> 00:02:10,000
and because I let the computer try out every step size 10 times.
32
00:02:10,000 --> 00:02:12,000
The overall behavior is what we know.
33
00:02:12,000 --> 00:02:17,000
With the forward Euler method the error depends linearly on the step size,
34
00:02:17,000 --> 00:02:20,000
at least for step sizes that are small enough.
35
00:02:20,000 --> 00:02:25,000
But for very small step sizes you see--uh-oh--the error explodes,
36
00:02:25,000 --> 00:02:27,000
and we have many small steps.
37
00:02:27,000 --> 00:02:31,000
Round-off error accumulates and destroys our result.
38
00:02:31,000 --> 00:02:35,000
What you should be learning from this reducing the step size
39
00:02:35,000 --> 00:02:37,000
only takes you so far.
40
00:02:37,000 --> 00:02:40,000
There is a threshold after which the error goes up again
41
00:02:40,000 --> 00:02:42,000
and may become severe.
42
00:02:42,000 --> 00:02:45,000
Let's look at the same data from a different viewpoint.
43
00:02:45,000 --> 00:02:49,000
Let's note that the error depends on the step size,
44
00:02:49,000 --> 00:02:53,000
but let's plot how the error depends on the number of steps.
45
00:02:53,000 --> 00:02:55,000
This is what we get.
46
00:02:55,000 --> 00:02:58,000
For a small number of steps, the step size is large,
47
00:02:58,000 --> 00:03:02,000
and the forward Euler method introduces a huge error,
48
00:03:02,000 --> 00:03:04,000
which then decreases in a linear fashion,
49
00:03:04,000 --> 00:03:09,000
but as the number of steps grows, we get more and more error due to round off.
50
00:03:09,000 --> 00:03:14,000
The simple model for this behavior would be to say that the error introduced by rounding
51
00:03:14,000 --> 00:03:16,000
is proportional to the number of steps
52
00:03:16,000 --> 99:59:59,000
and hence inversely proportional to the step size plus h.
PK iA@ S Unit_7_-_Advanced_Applications_of_Numerical_Methods/12-Simple_Roundoff_Error.en.srt1
00:00:00,000 --> 00:00:03,000
Now that we have an idea about the effect of round-off error,
2
00:00:03,000 --> 00:00:05,000
we can actually compute some estimates.
3
00:00:05,000 --> 00:00:10,000
For the forward Euler method, the total error is composed of two parts--
4
00:00:10,000 --> 00:00:15,000
the part due to the Euler method, the global truncation error, if you remember that term.
5
00:00:15,000 --> 00:00:19,000
A constant times the step size plus the error due to round off,
6
00:00:19,000 --> 00:00:21,000
which seems to be proportionate to the number of steps
7
00:00:21,000 --> 00:00:24,000
and, hence, inversely proportional to the step size
8
00:00:24,000 --> 00:00:27,000
Epsilon denotes the relative size of the error.
9
00:00:27,000 --> 00:00:32,000
It's something like 10⁻⁷ or 10⁻¹⁶.
10
00:00:32,000 --> 00:00:36,000
By which fraction does our number tend to change due to round off?
11
00:00:36,000 --> 00:00:38,000
Again, we should see proportionality here.
12
00:00:38,000 --> 00:00:41,000
If we have twice the amount of error in each step.
13
00:00:41,000 --> 00:00:44,000
We should be seeing twice the amount of total error,
14
00:00:44,000 --> 00:00:50,000
even though, again, the error induced in the earlier steps is growing over time.
15
00:00:50,000 --> 00:00:54,000
D is some proportionality constant that we don't know yet
16
00:00:54,000 --> 00:00:57,000
and that depends on the problem as does C.
17
00:00:57,000 --> 00:01:03,000
So, there is one component to the error that grows from zero to infinity on principle.
18
00:01:03,000 --> 00:01:08,000
There is another component to the error that starts at infinity and decays to zero.
19
00:01:08,000 --> 00:01:10,000
We have to have a minimum in between.
20
00:01:10,000 --> 00:01:12,000
This is what we saw in the experiment.
21
00:01:12,000 --> 00:01:16,000
It's called the step size at that minimum hmin.
22
00:01:16,000 --> 00:01:21,000
It does not make sense to use more steps, that is, to use a smaller step size,
23
00:01:21,000 --> 00:01:24,000
because the error will be growing again due to round-off.
24
00:01:24,000 --> 00:01:26,000
Now, here is a quiz.
25
00:01:26,000 --> 00:01:27,000
For the two types of floating point numbers that we're typically dealing with,
26
00:01:27,000 --> 00:01:31,000
we can specify the value of epsilon.
27
00:01:31,000 --> 00:01:35,000
It's about 10⁻⁷ for single precision.
28
00:01:35,000 --> 00:01:39,000
In C and related languages, that's called float.
29
00:01:39,000 --> 00:01:43,000
It's about 10⁻¹⁶ for double position.
30
00:01:43,000 --> 00:01:48,000
Python works at double position, and the arrays of numpy can be configured
31
00:01:48,000 --> 00:01:50,000
as to which position they should be using.
32
00:01:50,000 --> 00:01:52,000
Now here comes the quiz.
33
00:01:52,000 --> 00:01:54,000
Say we do an experiment at single position
34
00:01:54,000 --> 00:01:59,000
and find that the value of hmin amounts to 10⁻⁴.
35
00:01:59,000 --> 00:02:03,000
If we try to solve the same problem at double position,
36
00:02:03,000 --> 00:02:06,000
what do you expect to be the value of hmin?
37
00:02:06,000 --> 00:02:13,000
Pick one of these: 3 10⁻⁹, 4 10⁻¹⁰, 5 10⁻¹¹?
38
00:02:13,000 --> 00:02:21,000
Here is a hint: at the minimum, it turns out that the contribution of both parts is the same.
39
00:02:21,000 --> 99:59:59,000
Ch will equal Dε/h.
PK iAWXt t \ Unit_7_-_Advanced_Applications_of_Numerical_Methods/13-Simple_Roundoff_Error_Solution.en.srt1
00:00:00,000 --> 00:00:04,000
It's 3 x 10^-9. Let me show why.
2
00:00:04,000 --> 00:00:06,380
At the minimum, both contributions are equal.
3
00:00:06,380 --> 00:00:09,140
I'll explain why that is so in a second.
4
00:00:09,140 --> 00:00:13,490
So we find that the constant, C, x the value h minimum
5
00:00:13,490 --> 00:00:17,260
equals the constant D x ε over h minimum.
6
00:00:17,260 --> 00:00:22,440
From which follows that h minimum is the square root of Dε over C.
7
00:00:22,440 --> 00:00:26,400
ε is multiplied by a factor of 10^-9,
8
00:00:26,400 --> 00:00:30,110
which means that the value of h min, for double position,
9
00:00:30,110 --> 00:00:33,110
is multiplied by a factor of square root of 10^-9.
10
00:00:33,110 --> 00:00:37,620
Starting from the value of h min for single position,
11
00:00:37,620 --> 00:00:41,640
we have to multiply by the square root of 10^-9.
12
00:00:41,640 --> 00:00:49,890
In total, this is 10^-8.5. And that's about 3 x 10^-9.
13
00:00:49,890 --> 00:00:55,510
An important remark: This model was made for the forward Euler method.
14
00:00:55,510 --> 00:01:03,540
If we have a method of higher order, if we have Ch square, or even Ch 4 in here.
15
00:01:03,540 --> 00:01:08,980
So from this result we learn that as we go from single to double position,
16
00:01:08,980 --> 00:01:14,440
we may be using a step size that's about 4 to 5 orders of magnitude smaller,
17
00:01:14,440 --> 00:01:19,560
and hence the error will be smaller by 4 to 5 orders of magnitude.
18
00:01:19,560 --> 00:01:22,990
The price that we would be paying is that we have
19
00:01:22,990 --> 00:01:27,010
in the order of 10^4 to 10^5 more computation
20
00:01:27,010 --> 00:01:29,510
which looks a little prohibitive, but then again,
21
00:01:29,510 --> 00:01:34,490
with a method of higher order, you need half your steps than with the forward Euler method.
22
00:01:34,490 --> 00:01:38,090
So the affect of roundoff errors tend to be way smaller
23
00:01:38,090 --> 00:01:40,090
with more advanced methods.
24
00:01:40,090 --> 00:01:44,250
In case you're interested, let me explain why these 2 contributions
25
00:01:44,250 --> 00:01:46,330
are equal at the minimum.
26
00:01:46,330 --> 00:01:48,940
Otherwise, skip over the rest of this video.
27
00:01:48,940 --> 00:01:54,930
Substitute the step size h, by e^-u, with an appropriate value of u.
28
00:01:54,930 --> 00:02:03,140
Then the total error becomes C x e^-u + Dεe^u.
29
00:02:03,140 --> 00:02:05,930
Because we are dividing by that power.
30
00:02:05,930 --> 00:02:09,639
Now you can see that the 1st part becomes a falling exponential function,
31
00:02:09,639 --> 00:02:12,800
and the 2nd one becomes an increasing exponential function.
32
00:02:12,800 --> 00:02:18,200
Both fall, and due to symmetry, the minimum has to be where these 2 curves meet.
33
00:02:18,200 --> 00:02:21,260
A completely different approach to that problem is, of course,
34
00:02:21,260 --> 99:59:59,000
to compute the derivative and set it equal to 0.
PK iAUL;
;
O Unit_7_-_Advanced_Applications_of_Numerical_Methods/14-Solver_Components.en.srt1
00:00:00,000 --> 00:00:04,120
Choosing the right algorithm and choosing an appropriate floating-point format
2
00:00:04,120 --> 00:00:09,080
are just 2 of the issues we encounter when we implement a solver,
3
00:00:09,080 --> 00:00:11,760
a program to solve a differential equation.
4
00:00:11,760 --> 00:00:16,309
On a higher level, one may think about reusability and maintainability.
5
00:00:16,309 --> 00:00:19,840
This is my replica rendering of the Plug and Play style
6
00:00:19,840 --> 00:00:23,020
of implementation that's used in the book "Numerical Recipes."
7
00:00:23,020 --> 00:00:27,100
With this type of implementation, you can mix and match components,
8
00:00:27,100 --> 00:00:29,920
which comes in very handy for experiments.
9
00:00:29,920 --> 00:00:33,240
Let's have a closer look: 1 minor but vital component
10
00:00:33,240 --> 00:00:37,680
is the right hand side; the right hand side of the differential equation, that is.
11
00:00:37,680 --> 00:00:40,320
The rate of change. It's used by the algorithm,
12
00:00:40,320 --> 00:00:46,060
the algorithm being, for instance, the forward Euler method, Hunt's method, and similar.
13
00:00:46,060 --> 00:00:48,370
The algorithm is called by the stepper.
14
00:00:48,370 --> 00:00:52,420
which controls which steps are being made at which size.
15
00:00:52,420 --> 00:00:55,130
The stepper is responsible to adapt the step sizes
16
00:00:55,130 --> 00:00:59,610
and maybe to reduce the steps when it finds the error has grown too large.
17
00:00:59,610 --> 00:01:03,790
So the algorithm is only able to do 1 single step.
18
00:01:03,790 --> 00:01:06,100
The stepper calls the algorithm.
19
00:01:06,100 --> 00:01:08,810
The driver is the interface to the user,
20
00:01:08,810 --> 00:01:11,590
which, of course, is a program and not a human.
21
00:01:11,590 --> 00:01:14,730
The driver can be sent commands to start and stop,
22
00:01:14,730 --> 00:01:18,400
and eventually the driver returns the data to the user.
23
00:01:18,400 --> 00:01:21,630
To this end, it has a subobject that's called output.
24
00:01:21,630 --> 00:01:26,460
The output may collect values at predefined instances of time,
25
00:01:26,460 --> 00:01:32,630
or it may interpolate the data it receives from the stepper to regular time intervals.
26
00:01:32,630 --> 00:01:36,300
The stepper will pick smaller and larger time intervals.
27
00:01:36,300 --> 99:59:59,000
The output may interpolate them to produce data at regular time intervals.
PK iA!8 T Unit_7_-_Advanced_Applications_of_Numerical_Methods/15-Responsible_Components.en.srt1
00:00:00,000 --> 00:00:04,000
Given the division of labor in the solver as discussed before,
2
00:00:04,000 --> 00:00:08,000
which component would be affected by each of these actions?
3
00:00:08,000 --> 00:00:11,000
To use an implicit method? To stop at time 42?
4
00:00:11,000 --> 00:00:15,000
To determine the value at time 13.000?
5
00:00:15,000 --> 00:00:17,000
And to set a different value for the mass?
6
00:00:17,000 --> 00:00:19,000
Enter the first letter of that component.
7
00:00:19,000 --> 00:00:22,000
The driver, the stepper, the algorithm, the right-hand side,
8
00:00:22,000 --> 99:59:59,000
the output, or is it unclear?
PK iA[E% ] Unit_7_-_Advanced_Applications_of_Numerical_Methods/16-Responsible_Components_Solution.en.srt1
00:00:00,000 --> 00:00:03,000
Setting a different value for the mass is obviously
2
00:00:03,000 --> 00:00:06,000
something that has to happen in the right-hand side.
3
00:00:06,000 --> 00:00:09,000
The differential equations, as such, contains our physical model.
4
00:00:09,000 --> 00:00:12,000
Determining the value at one precise instant of time
5
00:00:12,000 --> 00:00:14,000
is a job for the output.
6
00:00:14,000 --> 00:00:18,000
Most probably it has to interpolate between values before and after,
7
00:00:18,000 --> 00:00:21,000
because there is no step at exactly this instant of time.
8
00:00:21,000 --> 00:00:23,000
So, this is O.
9
00:00:23,000 --> 00:00:25,000
Stopping at the given time--that's a job for the driver.
10
00:00:25,000 --> 00:00:27,000
It will stop the stepper.
11
00:00:27,000 --> 00:00:29,000
Using an implicit method is tricky.
12
00:00:29,000 --> 00:00:31,000
That's unclear.
13
00:00:31,000 --> 00:00:34,000
Somebody would be responsible to solve the equation.
14
00:00:34,000 --> 00:00:36,000
Who would that be?
15
00:00:36,000 --> 99:59:59,000
There is no explicit right-hand side of that equation.
PK iAs P Unit_7_-_Advanced_Applications_of_Numerical_Methods/17-Parallel_Computing.en.srt1
00:00:00,000 --> 00:00:04,510
1 decision that has to be made early on in the interpretation of a server
2
00:00:04,510 --> 00:00:06,760
is whether or not to use parallel computing.
3
00:00:06,760 --> 00:00:09,540
There are a range of ways for doing so
4
00:00:09,540 --> 00:00:12,190
and all of these can be combined with each other.
5
00:00:12,190 --> 00:00:15,310
We have vector units on the CPU,
6
00:00:15,310 --> 00:00:18,350
which are units that can operate on vectors;
7
00:00:18,350 --> 00:00:23,620
apply the same operation to say, 4 different values, at the same time.
8
00:00:23,620 --> 00:00:29,860
Multicore CPUs allow us to do several very complex things at the same time.
9
00:00:29,860 --> 00:00:32,479
The same is true for multiprocessor systems.
10
00:00:32,479 --> 00:00:38,270
GPUs, graphics chips, are optimized for working with huge arrays of pixels,
11
00:00:38,270 --> 00:00:40,270
but they can work with everything else as well.
12
00:00:40,270 --> 00:00:43,990
Which means, of course, they are also good at working with matrices.
13
00:00:43,990 --> 00:00:47,680
It's rather difficult to speed up servers for ordinary differential equations.
14
00:00:47,680 --> 00:00:50,140
The server does 1 step at a time,
15
00:00:50,140 --> 00:00:53,270
and every following step depends on the step before.
16
00:00:53,270 --> 00:00:55,770
How could you do some of these in parallel?
17
00:00:55,770 --> 00:00:57,770
That does not seem to work well.
18
00:00:57,770 --> 00:01:01,700
Speeding up a differential equation server for ordinary differential equations
19
00:01:01,700 --> 00:01:04,599
with the help of parallel computing may work,
20
00:01:04,599 --> 00:01:09,840
for instance in a situation where you have lots of isolated entities
21
00:01:09,840 --> 00:01:15,260
which rarely influence each other, such as the planets orbiting the sun.
22
00:01:15,260 --> 00:01:17,900
As you have seen in unit 6, it's of vital importance
23
00:01:17,900 --> 00:01:20,960
to speed up servers for partial differential equations.
24
00:01:20,960 --> 00:01:24,930
And in this case, parallel computing really comes to help.
25
00:01:24,930 --> 00:01:28,890
1 option for this is to use spatial subdivision.
26
00:01:28,890 --> 00:01:33,200
Each processor of a multicore system can handle its own spatial domain.
27
00:01:33,200 --> 00:01:37,730
Treating the borders between the different domains right requires communication.
28
00:01:37,730 --> 00:01:40,340
This can become an expensive overhead.
29
00:01:40,340 --> 00:01:42,380
We'll look into that in the next quiz.
30
00:01:42,380 --> 00:01:45,710
When we apply finite elements or work with implicit servers
31
00:01:45,710 --> 00:01:49,980
for partial differential equations, we often find huge matrices.
32
00:01:49,980 --> 99:59:59,000
They are an ideal workload for GPUs.
PK iA
x T Unit_7_-_Advanced_Applications_of_Numerical_Methods/18-Communication_Overhead.en.srt1
00:00:00,000 --> 00:00:04,000
Now let's look into the communication overhead required by parallel computing.
2
00:00:04,000 --> 00:00:08,000
Assume that we have n-squared course or processors,
3
00:00:08,000 --> 00:00:13,000
each of which handles its own square of our complete domain.
4
00:00:13,000 --> 00:00:18,000
The complete domain should consist of l l cells.
5
00:00:18,000 --> 00:00:21,000
What would be a reasonable model for the time taken by the computation?
6
00:00:21,000 --> 00:00:25,000
A constant C by the number of course or processors and squared
7
00:00:25,000 --> 00:00:28,000
times the side length of the total domain plus another constant
8
00:00:28,000 --> 00:00:30,000
times n times l?
9
00:00:30,000 --> 00:00:36,000
Or constant times l-squared divided by n-squared plus another constant times (n - 1)?
10
00:00:36,000 --> 00:00:41,000
Or the constant times n times l-squared minus a constant times n times l?
11
00:00:41,000 --> 00:00:48,000
A constant times l-squared over n-squared plus a constant times (n - 1) times l-squared?
12
00:00:48,000 --> 99:59:59,000
Pick one.
PK iA1xN N ] Unit_7_-_Advanced_Applications_of_Numerical_Methods/19-Communication_Overhead_Solution.en.srt1
00:00:00,000 --> 00:00:03,770
The second one would be a good choice, and I'll explain why.
2
00:00:03,770 --> 00:00:08,000
The total amount of computation to be done is proportional to l²,
3
00:00:08,000 --> 00:00:10,360
the number of cells in the domain.
4
00:00:10,360 --> 00:00:12,930
So if we have n² course or processors,
5
00:00:12,930 --> 00:00:17,220
the amount of work for each of these is proportional to l² over n².
6
00:00:17,220 --> 00:00:21,520
Which means that the first term describes the time taken by the computation alone.
7
00:00:21,520 --> 00:00:24,440
But then the processors, or the course,
8
00:00:24,440 --> 00:00:27,270
have to exchange data at the boundaries.
9
00:00:27,270 --> 00:00:32,640
Add that, and the amount of that data is proportional to n minus 1.
10
00:00:32,640 --> 00:00:35,600
We have n minus 1 vertical boundary lines
11
00:00:35,600 --> 00:00:38,660
and n minus 1 horizontal boundary lines,
12
00:00:38,660 --> 00:00:43,880
and it is proportional to l, because we have l cells along each boundary line.
13
00:00:43,880 --> 00:00:47,020
So of these expressions, the second one makes the most sense.
14
00:00:47,020 --> 00:00:49,300
Footnote: When we use this expression,
15
00:00:49,300 --> 00:00:52,540
we'll find that the optimum number of processors
16
00:00:52,540 --> 00:00:57,280
is proportional to the cubed root of the number of cells, l².
17
00:00:57,280 --> 00:01:00,250
This is a little surprising on first sight.
18
00:01:00,250 --> 00:01:03,170
You may expect that the optimum number is proportional
19
00:01:03,170 --> 00:01:06,960
to the number of cells, not proportional to its cubed root.
20
00:01:06,960 --> 99:59:59,000
This relationship is easy to verify by looking at the minimum of this expression.
PK iAO P Unit_7_-_Advanced_Applications_of_Numerical_Methods/20-Ready-Made_Solvers.en.srt1
00:00:00,000 --> 00:00:02,000
Instead of developing your own,
2
00:00:02,000 --> 00:00:05,000
in many cases you can, of course, use somebody else's solver.
3
00:00:05,000 --> 00:00:08,000
There are lots of ready-made software for that.
4
00:00:08,000 --> 00:00:11,000
First of all, there is general mathematical software--
5
00:00:11,000 --> 00:00:15,000
application packages, which can do differential equations,
6
00:00:15,000 --> 00:00:17,000
but which also can do lots of other stuff.
7
00:00:17,000 --> 00:00:22,000
First of all, MATLAB and its clones, such as Octave or SyLab,
8
00:00:22,000 --> 00:00:26,000
Wolfram Alpha if you manage to formulate your problem in one line,
9
00:00:26,000 --> 00:00:29,000
Mathematica, and tons of others.
10
00:00:29,000 --> 00:00:31,000
Second, there is dedicated modeling and simulation software,
11
00:00:31,000 --> 00:00:36,000
such as Simulink, which comes with MATLAB, COMSOL, Multiphysics,
12
00:00:36,000 --> 00:00:40,000
and Wolfram Alpha System Modeler, and PyDS Tool,
13
00:00:40,000 --> 00:00:42,000
which is based on, as the name says, Python.
14
00:00:42,000 --> 00:00:47,000
Then there are lots of code libraries that you can use to build your own application.
15
00:00:47,000 --> 00:00:51,000
One important library to mention here is scipy integrate--
16
00:00:51,000 --> 00:00:53,000
obviously for Python again.
17
00:00:53,000 --> 00:00:58,000
To create plausible motion for games, there are a number of physics engines.
18
00:00:58,000 --> 00:01:03,000
Finally, there are cook books providing snippets of code to be reused or adapted.
19
00:01:03,000 --> 99:59:59,000
The most prominent one of these are the Numerical Recipes.
PK iA)Z[ [ P Unit_7_-_Advanced_Applications_of_Numerical_Methods/21-Entering_Equations.en.srt1
00:00:00,000 --> 00:00:05,000
Most solvers available on the market do not accept the differential equation in this form--
2
00:00:05,000 --> 00:00:09,000
the second-order differential equation, which is not an explicit form.
3
00:00:09,000 --> 00:00:16,000
What you instead have to provide is a function that returns the first derivative, the rate of change.
4
00:00:16,000 --> 00:00:22,000
The trick is to turn this differential equation into two differential equations of first order.
5
00:00:22,000 --> 00:00:28,000
Here is a hint for you--the first equation should be that the derivative of x with respect to time
6
00:00:28,000 --> 00:00:34,000
is some new variable called y, and now you have to provide the second equation
7
00:00:34,000 --> 00:00:36,000
for the derivative of y with respect to time.
8
00:00:36,000 --> 00:00:44,000
Should that be 7-3y+4x or 7-3 times the derivative of x plus 4x
9
00:00:44,000 --> 99:59:59,000
or 1/37 minus the second derivative plus 4x or 1/37 minus the first derivative of y plus 4x.
PK iA1 Y Unit_7_-_Advanced_Applications_of_Numerical_Methods/22-Entering_Equations_Solution.en.srt1
00:00:00,000 --> 00:00:02,000
The first option is the correct one.
2
00:00:02,000 --> 00:00:07,000
You bring 3 times the derivative of x minus 4x over to the right hand side
3
00:00:07,000 --> 00:00:15,000
and write the derivative of x as y, so this becomes y and the second derivative of becomes
4
00:00:15,000 --> 00:00:24,000
the first derivative of y, so we end up with the first derivative of y=7-3y+4x.
5
00:00:24,000 --> 00:00:28,000
This is the standard trick to convert higher order differential equations
6
00:00:28,000 --> 00:00:32,000
into first order systems of differential equations.
7
00:00:32,000 --> 00:00:38,000
You maybe asking--why are we not using the second one. The solver needs the rate of change.
8
00:00:38,000 --> 00:00:44,000
It provides x and y, the current values so that function should depend on x and y
9
00:00:44,000 --> 99:59:59,000
and not on the rate of change, which is to be determined.
PK iAi a Unit_7_-_Advanced_Applications_of_Numerical_Methods/23-Verification_Calibration_Validation.en.srt1
00:00:00,000 --> 00:00:03,000
Writing your own software or using somebody else's software
2
00:00:03,000 --> 00:00:06,000
is one thing--what about getting it right?
3
00:00:06,000 --> 00:00:10,000
Generally, one distinguishes between three different aspects of that.
4
00:00:10,000 --> 00:00:15,000
First is verification, checking whether your code really implements your model.
5
00:00:15,000 --> 00:00:22,000
Second, calibration. Ensuring the parameters of your model so as to reproduce given observations.
6
00:00:22,000 --> 00:00:28,000
And the third one is validation. Check whether the real word reflects the predictions of your model.
7
00:00:28,000 --> 99:59:59,000
Of course, within a certain application domain and to a certain accuracy.
PK iAk97 7 S Unit_7_-_Advanced_Applications_of_Numerical_Methods/24-Domain_Identification.en.srt1
00:00:00,000 --> 00:00:05,000
Verification, calibration, validation--to which of these three aspects
2
00:00:05,000 --> 00:00:07,000
do these three problems belong.
3
00:00:07,000 --> 00:00:14,000
The first one--our actual car has a braking time of more than 6 seconds. Think about Unit 5.
4
00:00:14,000 --> 00:00:19,000
The second one--the orbit of the spacecraft around earth is not closed. Think about Unit 2.
5
00:00:19,000 --> 99:59:59,000
And the third on--without harvesting the amount of fish reaches 10⁷ tons. Think about Unit 4.
PK iAc4 4 \ Unit_7_-_Advanced_Applications_of_Numerical_Methods/25-Domain_Identification_Solution.en.srt1
00:00:00,000 --> 00:00:05,000
If the orbit of the spacecraft is not closed even though it has to be in the model,
2
00:00:05,000 --> 00:00:09,000
we've got some problem with implementations, so that's a verification issue.
3
00:00:09,000 --> 00:00:13,000
The maximum carrying capacity was a model parameter for the fish model,
4
00:00:13,000 --> 00:00:17,000
so in this case the calibration was done right.
5
00:00:17,000 --> 00:00:23,000
And the other one is a prediction of the model doesn't work out well, so that's a validation issue.
PK iAE M Unit_7_-_Advanced_Applications_of_Numerical_Methods/26-Order_of_Solver.en.srt1
00:00:00,000 --> 00:00:05,000
Here comes one idea about how to test the implementation of a differential equation solver.
2
00:00:05,000 --> 00:00:10,000
We provide the solver that we've never covered before that looks pretty complex.
3
00:00:10,000 --> 00:00:14,000
The question is what's the order of that solver. Enter that order here.
4
00:00:14,000 --> 00:00:19,000
The problem that we're using is again the satellite the takes 24 hours to orbit the earth
5
00:00:19,000 --> 00:00:25,000
and we are comparing its final position to its start position and plot the arrow.
6
00:00:25,000 --> 99:59:59,000
From that graph, determine the order of that solver.
PK iA?= = M Unit_7_-_Advanced_Applications_of_Numerical_Methods/27-Order_of_Solver.en.srt1
00:00:00,000 --> 00:00:08,000
Look closely--both axis are logarithmic and you can see as the step size increases by a factor of 10,
2
00:00:08,000 --> 00:00:12,000
the error increases by a factor of 10⁵.
3
00:00:12,000 --> 00:00:19,000
Remember, forward Euler was of order 1, Heun was of order 2. This is the solver of order 5.
4
00:00:19,000 --> 00:00:24,000
Read this the other way around as you shrink the step size by a factor of 10,
5
00:00:24,000 --> 00:00:31,000
the error decreases by a factor of 10⁵--that's a huge gain and you can see that this solver,
6
00:00:31,000 --> 00:00:37,000
when it takes steps of 3 minutes around earth gets back to where it started within just millimeters.
7
00:00:37,000 --> 00:00:42,000
This is quite an efficient method. It's called the Runge-Kotta feedback method.
8
00:00:42,000 --> 00:00:46,000
Now assume that with all these strange numbers we mistyped one--
9
00:00:46,000 --> 00:00:54,000
mainly we didn't enter 56,430 but we entered 56,431.
10
00:00:54,000 --> 00:00:57,000
Let's look at the result. This doesn't look that great anymore.
11
00:00:57,000 --> 00:01:01,000
We've got an error that's well above 1 km, so getting these numbers right,
12
00:01:01,000 --> 99:59:59,000
which among other things belong to verification.
PK iAeYg g Z Unit_7_-_Advanced_Applications_of_Numerical_Methods/28-Prediction_vs._Understanding.en.srt1
00:00:00,000 --> 00:00:06,000
There is one fundamental issue about calibration and validation, namely that
2
00:00:06,000 --> 00:00:10,000
prediction does not always go hand in hand with understanding.
3
00:00:10,000 --> 00:00:14,000
Look at these two models for the motion of a planet in the solar system.
4
00:00:14,000 --> 00:00:18,000
Let's look at the left model--here the earth is fixed in space
5
00:00:18,000 --> 00:00:22,000
and close to the earth--thus the center of some large circle.
6
00:00:22,000 --> 00:00:27,000
This point moves around that circle and a point on that second circle,
7
00:00:27,000 --> 00:00:29,000
which then moves in this direction.
8
00:00:29,000 --> 00:00:32,000
This is called an epicycle, and we can continue that process.
9
00:00:32,000 --> 00:00:37,000
That point could be the center of another circle and so on and so on until finally
10
00:00:37,000 --> 00:00:40,000
we say that the planet is moving around that final circle.
11
00:00:40,000 --> 00:00:43,000
That's one of the antique models and let's have a look at the right one.
12
00:00:43,000 --> 00:00:49,000
Here the sun is at the center, the earth moves in a circle around the sun,
13
00:00:49,000 --> 00:00:51,000
and that planet moves in circle.
14
00:00:51,000 --> 00:00:58,000
Here comes the last quiz for you--which of these models allows a more flexible calibration
15
00:00:58,000 --> 00:01:03,000
and which of these two models allows a more accurate prediction of the apparent motion
16
00:01:03,000 --> 99:59:59,000
of the planet as seen from earth.
PK iAj j Z Unit_7_-_Advanced_Applications_of_Numerical_Methods/29-Prediction_vs._Understanding.en.srt1
00:00:00,000 --> 00:00:06,000
Maybe somewhat surprisingly, the left model that of Ptolemy has so many knobs to turn
2
00:00:06,000 --> 00:00:10,000
that it allows more flexible calibration and that it even allows
3
00:00:10,000 --> 00:00:13,000
more accurate prediction of the apparent motion.
4
00:00:13,000 --> 00:00:19,000
Footnote--in a modern perspective, this stack of epicycles actually is a Fourier series.
5
00:00:19,000 --> 00:00:24,000
Historically, it turned out that circles around the sun are not good enough.
6
00:00:24,000 --> 00:00:30,000
Kepler noticed that you have to use ellipses to make this really work and then we get a high accuracy.
7
00:00:30,000 --> 00:00:35,000
The lesson to be learned from that is you may want to replace a sophisticated model
8
00:00:35,000 --> 00:00:40,000
with accurate predictions by a far simpler model, less accurate predictions,
9
00:00:40,000 --> 00:00:42,000
just for the sake of understanding.
10
00:00:42,000 --> 99:59:59,000
This is Occam's Razor in action, trying to find the most simple explanation that does the job.
PK iA&b H Unit_7_-_Advanced_Applications_of_Numerical_Methods/30-Conclusion.en.srt1
00:00:00,000 --> 00:00:03,000
We started this Unit by looking into some fundamental ideas
2
00:00:03,000 --> 00:00:07,000
of finite elements and computational fluid dynamics.
3
00:00:07,000 --> 00:00:10,000
These are two related and highly application-oriented field
4
00:00:10,000 --> 00:00:13,000
that of particular use in mechanical engineering.
5
00:00:13,000 --> 00:00:18,000
The toy model of fluid dynamics left to the notion of deterministic chaos--
6
00:00:18,000 --> 00:00:22,000
the sensitive dependence on initial conditions, but we saw that randomness
7
00:00:22,000 --> 00:00:24,000
can be produced virtually out of nothing.
8
00:00:24,000 --> 00:00:27,000
Then we look into a range of topics from software engineering--
9
00:00:27,000 --> 00:00:33,000
how can we build modules that can be mixed and matched, how can we reuse software,
10
00:00:33,000 --> 00:00:36,000
and how do we calibrate and test simulations.
11
00:00:36,000 --> 00:00:40,000
These questions eventually led us to a little bit of philosophy,
12
00:00:40,000 --> 99:59:59,000
and now it's time to apply all of these in the final exam.
PK iA7Y Y L Unit_7_-_Advanced_Applications_of_Numerical_Methods/01-Welcome_Unit_7.en.srtPK iAJ1 S Unit_7_-_Advanced_Applications_of_Numerical_Methods/02-Finite_Element_Method.en.srtPK iA8ob J Unit_7_-_Advanced_Applications_of_Numerical_Methods/03-Hanging_Rope.en.srtPK iAeR R S Unit_7_-_Advanced_Applications_of_Numerical_Methods/04-Hanging_Rope_Solution.en.srtPK iAnS S L Unit_7_-_Advanced_Applications_of_Numerical_Methods/05-Fluid_Dynamics.en.srtPK iAQE E Q 5 Unit_7_-_Advanced_Applications_of_Numerical_Methods/06-Force_from_Pressure.en.srtPK iA'' Q 3> Unit_7_-_Advanced_Applications_of_Numerical_Methods/07-Force_from_Pressure.en.srtPK iA;,j4b
b
O I@ Unit_7_-_Advanced_Applications_of_Numerical_Methods/08-The_Lorenz_System.en.srtPK iAi3N N V K Unit_7_-_Advanced_Applications_of_Numerical_Methods/09-Lorenz_and_Forward_Euler.en.srtPK iAO O V N Unit_7_-_Advanced_Applications_of_Numerical_Methods/10-Lorenz_and_Forward_Euler.en.srtPK iArL`g L X Unit_7_-_Advanced_Applications_of_Numerical_Methods/11-Roundoff_Error.en.srtPK iA@ S l Unit_7_-_Advanced_Applications_of_Numerical_Methods/12-Simple_Roundoff_Error.en.srtPK iAWXt t \ J{ Unit_7_-_Advanced_Applications_of_Numerical_Methods/13-Simple_Roundoff_Error_Solution.en.srtPK iAUL;
;
O 8 Unit_7_-_Advanced_Applications_of_Numerical_Methods/14-Solver_Components.en.srtPK iA!8 T Unit_7_-_Advanced_Applications_of_Numerical_Methods/15-Responsible_Components.en.srtPK iA[E% ] Unit_7_-_Advanced_Applications_of_Numerical_Methods/16-Responsible_Components_Solution.en.srtPK iAs P 0 Unit_7_-_Advanced_Applications_of_Numerical_Methods/17-Parallel_Computing.en.srtPK iA
x T R Unit_7_-_Advanced_Applications_of_Numerical_Methods/18-Communication_Overhead.en.srtPK iA1xN N ] N Unit_7_-_Advanced_Applications_of_Numerical_Methods/19-Communication_Overhead_Solution.en.srtPK iAO P Unit_7_-_Advanced_Applications_of_Numerical_Methods/20-Ready-Made_Solvers.en.srtPK iA)Z[ [ P } Unit_7_-_Advanced_Applications_of_Numerical_Methods/21-Entering_Equations.en.srtPK iA1 Y F Unit_7_-_Advanced_Applications_of_Numerical_Methods/22-Entering_Equations_Solution.en.srtPK iAi a Unit_7_-_Advanced_Applications_of_Numerical_Methods/23-Verification_Calibration_Validation.en.srtPK iAk97 7 S % Unit_7_-_Advanced_Applications_of_Numerical_Methods/24-Domain_Identification.en.srtPK iAc4 4 \ Unit_7_-_Advanced_Applications_of_Numerical_Methods/25-Domain_Identification_Solution.en.srtPK iAE M { Unit_7_-_Advanced_Applications_of_Numerical_Methods/26-Order_of_Solver.en.srtPK iA?= = M Unit_7_-_Advanced_Applications_of_Numerical_Methods/27-Order_of_Solver.en.srtPK iAeYg g Z @ Unit_7_-_Advanced_Applications_of_Numerical_Methods/28-Prediction_vs._Understanding.en.srtPK iAj j Z Unit_7_-_Advanced_Applications_of_Numerical_Methods/29-Prediction_vs._Understanding.en.srtPK iA&b H Unit_7_-_Advanced_Applications_of_Numerical_Methods/30-Conclusion.en.srtPK $