1
00:00:03,000 --> 00:00:10,000
For this homework, the first question is to write an implementation of the so-called trapezoidal method, which is a second order implicit Runge-Kutta method
2
00:00:10,000 --> 00:00:16,000
This method is effectively averaging a forward Euler step and a backward Euler step
3
00:00:16,000 --> 00:00:23,000
While generally, both backward and forward Euler provide a lot of error in the solution, trapezoidal does a little bit better
4
00:00:23,000 --> 00:00:31,000
As some of you may have noticed, it is called trapezoidal method due to its similarities with the trapezoidal integration method you learned in calculus
5
00:00:31,000 --> 00:00:34,000
Heres my Matlab implementation of the trapezoidal method
6
00:00:34,000 --> 00:00:38,000
As inputs, it takes an initial condition, a step size, and a final time
7
00:00:38,000 --> 00:00:43,000
Functionally, this code is very similar to both backward and forward Euler, its just what we do with each step
8
00:00:43,000 --> 00:00:47,000
In particular, this line right here
9
00:00:47,000 --> 00:00:56,000
This takes the average of a forward Euler step and a backward Euler step
10
00:00:56,000 --> 00:01:00,000
Notice that you have to do a forward Euler step to get x to the next time step in order to see the slope there
11
00:01:00,000 --> 00:01:02,000
This is very common in implicit methods
12
00:01:02,000 --> 00:01:11,000
Notice that while this line of code is just simply the code for forward Euler, I could have just as easily called the forward Euler implementation that I wrote for a previous homework
13
00:01:11,000 --> 00:01:16,000
Below the trapezoidal code, for completeness, I provide the simple harmonic oscillator code as well
14
00:01:16,000 --> 00:01:18,000
This is also provided in earlier homework solutions
15
00:01:18,000 --> 00:01:23,000
Note that k =2, m = 0.5, and beta = 0 are the parameters that well use throughout this homework
16
00:01:23,000 --> 00:01:33,000
Part a asks us to find what x and v will be at time 0.5, starting from [-1, -2], using a timestep of h = 0.1, and the given parameters
17
00:01:33,000 --> 00:01:43,000
If we do this experiment using the code I just showed, you get x = -1.3811 and v = 0.6211 at time 0.5
18
00:01:43,000 --> 00:01:52,000
For part b, we needed to generate a 500-point trajectory of the same ODE system with a timestep of 0.01, and see which of the following describes the trajectory best
19
00:01:52,000 --> 00:01:58,000
This plot is precisely that trajectory, and as you can see, an ellipse is a good description of this
20
00:01:58,000 --> 00:02:11,000
For part c, were supposed to keep in mind that there is no damping in the system that is, that friction is equal to 0 and decide whether the trajectories generated with forward and backward Euler in Homework 5.4 are more accurate or less accurate than the trajectory generated with the trapezoidal method
21
00:02:11,000 --> 00:02:15,000
The best way to approach this is to simply plot all three on the same plot
22
00:02:15,000 --> 00:02:17,000
If we do that, we get this plot
23
00:02:17,000 --> 00:02:20,000
All three trajectories started at this point
24
00:02:20,000 --> 00:02:25,000
The blue trajectory is backward Euler, and it quickly spirals in to the origin
25
00:02:25,000 --> 00:02:33,000
The red trajectory is trapezoidal, and then does the elliptical pattern it should, and the green trajectory is forward Euler, and spirals out very quickly
26
00:02:33,000 --> 00:02:39,000
So from this picture, it does seem that the trapezoidal method is providing a more accurate solution, or a better trajectory for these dynamics
27
00:02:39,000 --> 00:02:46,000
You can actually formally prove that the trapezoidal method is order h cubed locally, whereas backward and forward Euler are both order h squared
28
00:02:46,000 --> 00:02:49,000
So the red really is a better representation of the trajectory
29
00:02:49,000 --> 00:03:00,000
For parts d and e, we need to generate a 5,000-point trajectory of the same ODE system and the same initial conditions using a timestep of 0.01, as we did in part b
30
00:03:00,000 --> 00:03:04,000
And we want to see qualitatively if theres a difference between this plot and the plot you generated using 500 points
31
00:03:04,000 --> 00:03:10,000
Describe that difference, if there is one; and if there is a difference, we need to try to diagnose that difference
32
00:03:10,000 --> 00:03:16,000
In this plot, I have plotted both the 500-point trajectory, in red, and the 5,000-point trajectory in blue
33
00:03:16,000 --> 00:03:22,000
As you can see, the trajectory has slightly fattened; that is, the trajectory has gotten a little bit wider and takes up a little bit more of phase space
34
00:03:22,000 --> 00:03:29,000
Whats happening here is that the numerical error caused in the trapezoidal method, while much, much smaller than either forward or backward Euler, still exists
35
00:03:29,000 --> 00:03:40,000
This numerical error is violating the conservation of energy properties of a simple harmonic oscillator with no friction, and since the system is a conservative system, a symplectic integrator should be used, or some type of integrator that conserves energy
36
00:03:40,000 --> 00:03:50,000
If youre brave enough to do the suggested problem, and code up RK4, you can see that this problem, where energy is either being gained or lost in the system, persists even for RK4 and RK5
37
00:03:50,000 --> 00:03:56,000
What this should show you is that, while RK4 is the go-to workhorse in dynamical systems, its not appropriate for every dynamical system
38
00:03:56,000 --> 00:04:02,000
Even a simple system like the simple harmonic oscillator with no damping can give RK4 grief
39
00:04:02,000 --> 00:04:11,000
Just because you have order h to the fourth global error propagation doesnt mean youll get accurate solutions, especially in the cases where special properties need to be preserved like conservation of energy
40
00:04:11,000 --> 00:04:19,000
This is something to be aware of, and a reason to know how to code up an ODE solver, and not just use the out-of-the-box method that comes with Mathematica or Matlab or Numpy, for example
41
00:04:19,000 --> 00:04:25,000
So the answer to question d is that the trajectory fattened, taking up more space
42
00:04:25,000 --> 00:04:33,000
And the answer to question e is that the numerical error is violating the conservation of energy property of a simple harmonic oscillator with beta = 0
43
00:04:33,000 --> 00:04:39,000
Since the system is a conservative system, a symplectic integrator should have been used