In 5.4, youll actually get your hands dirty in writing an ODE solver
Youll write forward Euler and backward Euler
Forward Euler is also very commonly called explicit Euler because its an explicit solver
You will then use the explicit Euler solver you wrote to generate trajectories of the simple harmonic oscillator for different parameter sets
Heres my Matlab implementation of forward Euler
It takes in an initial condition x0, a step size h, a final time, and a showPlot parameter which I just need later on as a flag whether or not to plot the trajectories
You can pretty much just ignore this one if you need to
The first step is just to make xCurrent, which is a helper variable, equal to the initial condition, and make the x0 the first element in the trajectory
Notice that, since Matlab starts array indexing at one, x0 will actually be the first element in the trajectory
That is, x0 will be trajectory[1], so x1 will be trajectory[2]
These off-by-one errors can often be very confusing
We then just have this main loop
All it does is take xCurrent, which is that helper variable I used above, starts from that position, and takes one step in the direction of the differential equation
So the next element in the trajectory is simply the current position plus a tiny step in the direction of the ODE
We then update xCurrent and loop for as many times as we need to get to the final time
We then may or may not plot the trajectory, depending on that showPlot parameter
The simple harmonic oscillator is given by these two differential equations, and has three parameters: k, m, and beta
We then calculate x based on these differential equations evaluated at xCurrent and return x
Notice that k = 2, m = 0.5, and beta = 0 are precisely the parameters that we used in both Quiz 5.3 and will use in Homework 5.4 for the first problem
Its always a good idea when you write new software to desk check the software
That means to do the first iteration by hand and make sure that what you got by hand is what you get with the software
In Quiz 5.3, you used forward Euler to solve the simple harmonic oscillator for one time step using a time step of 0.1, starting at [-1, -2]
And if we do that, we see that we get -1.2 and -1.6, which is precisely what we got doing the desk check
So I would say that our software is good to go
The first problem on Homework 5.4 after weve developed this software of course is to generate a trajectory for 0.5 seconds starting at k = 2, m = 0.5, and beta = 0 from the initial condition [-1, -2]
Lets see what we get
So, starting from initial condition [-1, -2], using a step size of 0.1, and going to time, or final time, 0.5, and in this case we want to not plot the final trajectory, we get that the last point in the trajectory, or the trajectory at time 0.5 using this solver, is position -1.5283 and velocity 0.6246
That is the answer to question 1(a)
For part b, we want to keep all the parameters the same, and the same initial condition
We want to generate a 200-point trajectory instead of a 5-point trajectory, and well generate this trajectory using a time step of 0.1 as well as a time step of 0.11
If we do this, we get the following plot
In this plot, the blue curve is generated using a time step of 0.1, and the red curve is using a time step of 0.11
As you can see, and as youd expect from the lecture, the red curve spirals out faster
That is, if you take a bigger time step, youre going to step farther off the trajectory
And in this case, youll spiral out faster
See the lecture for a very nice visualization of this
From this plot, we can see that using the delta t of 0.11 spiraled out faster
Question 2 asks us to write backward Euler
This is also called implicit Euler, as it is an implicit solver
My implementation is a Matlab implementation that takes the exact same inputs as forward Euler did: namely, an initial condition, a step size, a final time, and this plot parameter
My backward Euler code is identical to my forward Euler code except for this one key step
Remember, with an implicit solver, you approximate the next step and then you use that as your current guess
So in the case of backward Euler, this line is identical, except for, instead of evaluating the simple harmonic oscillator at xCurrent, you use xCurrent to do a forward Euler step, save that, and then you update the trajectory by using the slope at the forward Euler step
This can be a little bit confusing, but youre approximating a step into the future, pulling the direction at that new step back in time, and then using the direction of the new step to go forward in time
Liz provided a really nice visualization of whats happening here in the lecture, and if youre confused about this I encourage you to go back and see that
While implicit solvers seem like they may be a little bit confusing, and you dont seem like youre gaining much, theyll become a lot more important when we do industrial-strength integrators in Unit 6
Question 2, part a asks us to find what x at time 0.5 would be using a time step of 0.1, starting at [-1, -2] using our backwards Euler solver
If we do this, we see that x at time 0.5 will be -1.2451 and the velocity will be 0.6136
This is precisely the answer to question 2(a)
Question 2(b) asks us to repeat problem 1(b), but instead using the backward Euler solver we just created
Again, in this plot, the blue trajectory uses a step size of 0.1, and the red trajectory uses a step size of 0.11
As we would expect, the red trajectory using a larger step size spirals in faster
So this is the answer to question 2(b): the delta t of 0.11 spirals in faster than using a time step of 0.1
Problem 2(c) asks us to use both a forward and backward Euler solver, using a time step of 0.1, to generate two 50-point trajectories, and plot them both together to see what the difference is
Heres that plot: the blue uses backward Euler, and the red uses forward Euler
Clearly the difference here is that one is spiraling in and one is spiraling out
This should not be surprising from what we just did
Whats interesting here is that this is an undamped simple harmonic oscillator, yet one of these is growing in amplitude, and one of these is shrinking in amplitude
So the numerical solver is introducing dynamics that dont actually exist
But all you need to know for this problem is that the difference between these two is that one is spiraling in and one is spiraling out
This is the answer to question 2(c)