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)