Programming Assignment 2: Numerical Convergence ESM 3124 - Fall 2020 Due: 11-13-2020 Two results of your first programming assignment are functions that can be called with: dSdt = my_ode(t, S, g, L) [S, t] = eulers_method(@(t, S) my_ode(t, S, g, L),tspan, S0, delta_t) Several students observed that the total energy of the system, calculated using their numerical solution, grew with time. This apparent violation of the conservation of energy has nothing to do with the pendulum being studied, and instead depends on the details of the numerical solution. Notice that eulers_method(...,delta_t) is a function of delta_t, the time step chosen to simulate the motion. This assignment will focus on determining exactly how the solution depends on delta_t. The assignment is broken into 2 main parts. All submissions need to be in a single .zip folder labeled 'Your_Name_Pendulum_Convergence.zip'. 0.1 Convergence Analysis (10 pts) A necessary condition for calculating an accurate numerical solution is that your results are convergent. Conceptually, this means that as the time step ∆t used to simulate motion becomes smaller, we expect that the solution converges to a single state. If we have two solutions θi(t) and θi+1(t) that are calculated using time steps ∆ti and ∆ti+1, we can calculate the magnitude of the difference in the solutions at a fixed point in time ∆θi = |θi+1(t0)− θi(t0)|. If the solutions converge to a single value, then ∆θi → 0 as ∆ti → 0. A solution satisfying this property is said to be convergent. Now that you've finished reading this paragraph re-read this paragraph, it explains the entire assignment. While this may make us want to choose the smallest possible ∆t, accuracy comes with the trade-off of increased run time. Your goal in this section is to design and write a function that determines if your solution is convergent, and times how long it takes eulers_method() to run for numerous values of ∆t. Assess convergence using 2 fairly simple metrics 1) the final angle θf (t = tf ), and the final system energy Ef (t = tf ). Use your code from the first programming assignment as the starting point for this, but remove any parts that are not used in this assignment (e.g., the animation section). 1. Write a local function: [del_S, del_E, run_time] = convergence_analysis(..., delta_t_array) that measures the convergence of your eulers_method function by looking at the change in the final angle and energy as a function of delta_t. Your function must also measure run_time, how long it takes to simulate each trajectory. (See the help menu for tic and toc to learn how to time your code) * You will need to decide what most of the function inputs will be. However, one input argu- ment to convergence_analysis() must be a columnn vector delta_t_array that contains all values of ∆ti that you will simulate motion for. 0.2 Determine Valid ∆t (10 pts) Use the function written above to simulate trajectories for 30 logarithmically spaced values from ∆t = 10−5 seconds to ∆t = 10−1 seconds (pull up the help menu for logspace() for assistance). Plot how the solution converges as a function of ∆t, and what the associated run_times are. Using your plots, conclude what values of ∆t are useful for simulating motion. 1. For each value of ∆t, simulate 10 seconds of pendulum motion, then save the final values of the angle theta_f and total energy E_f. After calculating these values for each delta_t in delta_t_array, calculate the change in each value as the delta_t increases. Simulate a pendulum with L = 1m, M = 1kg, and initial conditions θ0 = pi/1.1 and θ˙0 = 0. 1 Figure 1: Example of convergence results 2. Plot the result of your calculations in a single figure with two subplots. One plot must show the convergence of the final angle ∆θf vs ∆t, and the second the final energy ∆Ef vs ∆t. (See loglog() for plotting logarithmically spaced values, and yyaxis for creating different left / right vertical axes). 3. Print a 200 dpi .png image of your plot called convergence_plot.png. Submission Format Make sure your code is properly commented to explain what each portion does, and how. In a single .zip folder named `Your_Name_Pendulum_Convergence.zip' you will submit: 1. A single .m file that contains your code and local functions described in Part 0.1. The TA and Professor Domann should be able to load this code into Matlab, press Run, and have your code execute successfully. 2. An image called convergence_plot.png that contains an image of your convergence analysis. 2 Recommended Code Update: Compare how you used ode23() and eulers_method() in the first programming assignment. odefun = @(t,S) my_ode(t,S,g,L); [t, S] = ode23(odefun, tspan, S0, options); [t, S] = eulers_method(odefun, tspan, S0, delta_t); These functions are structured to do the same thing, numerically simulate an ODE. Note that the inputs to both functions are nearly identical. The only information ode23() receives about the ODE it is simulating is odefun = @(t, S) my_ode(t, S, g, L), and the initial state S0. In other words, ode23() is flexible enough to make calculations using only those properties. Now consider how most people wrote eulers_method() function [t, S] = eulers_method(odefun, tspan, S0, delta_t)| %% Function Setup ... %% Run Problem Specific Loop for i = 1 : numel(t)-1 % rate of change of state dsdt = odefun(t(i), S(i,:)); % update angle and angular velocity theta(i+1) = theta(i) + dsdt(1)*delta_t; omega(i+1) = omega(i) + dsdt(2)*delta_t; % update state S(i+1,1) = theta(i+1); S(i+1,2) = omega(i+1); end %end loop end %end function While that code is very readable, it only works for a problem with 1 position variable (θ) and one velocity variable (ω). What if you want to simulate a 2 degree of freedom system like a double pendulum? Your code would no longer work. Fortunately there is a very easy change. Note that the code for % update angle ... merely serves as temporary storage, and is not actually needed. As a result the code can simultaneously be condensed, and made more flexible. function [t, S] = eulers_method(odefun, tspan, S0, delta_t)| %% Function Setup ... %% Run More Flexible Loop for i = 1 : numel(t)-1 % rate of change of state dsdt = odefun(t(i), S(i,:)); % update state S(i+1,:) = S(i, :) + dsdt * delta_t; end %end loop end %end function As long as you initialize S so it has the same number of columns that dsdt has elements, then this code will work. Note that the code S(i, :) is read `S at row i, all columns'. If you make this change, then when you want to simulate a new system the only code you need to change is the definition of odefun(), where you write your system of 1st order ODEs. In other words, your numerical method eulers_method() no longer depends on a single problem, and is instead a tool capable of simulating any properly structured ODE. I recommend doing this now before the third and final programming assignment. 3
欢迎咨询51作业君