Untar the package of C code on the class website named 'diffeq.tar'. This will create a subdirectory named 'diffeq' under the current directory. Change into that 'diffeq' subdirectory.
In this section we will simulate a simple pendulum with a large displacement angle and find the error in the small angle approximation for pendulums. Begin by opening the source code file pendulum.c in a text editor:
gedit pendulum.c &
Follow the code in the main progam as it extracts five parameters off of the command line, in order: the length of the pendulum in meters, the initial values for y0 and y1, and finally the stop and increment values of the independent variable x, which in this case is time in seconds. Note where the rk4 function is called near the end of the main program to invoke the fourth order Runge-Kutta algorithm to solve the differential equation. The function whose pointer is passed to the rk4 routine is defined earlier in the source code file, and is named pendulum(). Complete the source code lines indicated in the pendulum function to implement the pendulum equations. Save your source code file and compile with:
gcc -Wall pendulum.c diffeq.c -o pendulum -lm
First try running your compiled executable program and viewing its help line output:
./pendulum
Now run the program to simulate a 1 meter long pendulum that has been given an initial angular displacement of some small angle, say 0.1 radians, with:
./pendulum 1 0.1 0 10 0.01 >| pendulum.dat
Note this simulation will run for a simulated time of 10 seconds, and with a fairly fine time step of 10 msec. Since the system being solved is a simple one, we can afford a fine time step and still have very fast execution time.
Start another terminal window, change to the diffeq directory, and start a gnuplot process. The rk4 function will write to standard output as many state variables as there are in the system of differential equations that have been passed to it to solve, so in this case the pendulum.dat file will comprise three columns, the first time, the second the angular displacement, and the third the angular velocity. Use gnuplot to plot both the angular displacement and the angular velocity. You will need to make use of the 'using' keyword on the plot command line to select the pairs of columns that will form the x and y coordinates of each curve. You should end up with something like:
Back in the terminal running the bash shell, rerun the program with a few different initial angle values, which will generate new versions of pendulum.dat. Replot these new solutions in the gnuplot window with the 'replot' command.
Now we want to use the tool code from Homework 9 to measure the simulated period of the pendulum. In that homework assignment you had expanded interpolate.c to include functions to perform inverse interpolation. The solution code for that exercise has been included in this lab tar file so that you can link it to the interpolate_solve.c main program from Homework 9. Build the interpolate_solve tool in this directory with:
gcc -Wall -DNEWTON interpolate_solve.c interpolate.c data_file.c -o interpolate_solve -lm
Before we can measure the period, we have a slight incompatibility between the tools that will need to be resolved with a text editor. As detailed above, there are three columns of data output by the rk4 differential equation solving tool, but the interpolate_solve tool expects a two-column data file. So open the pendulum.dat file with the gedit editor and edit it down so that the interpolation tool will find the exact simulated time for the last event when the displacement angle crosses zero from a negative angle to a positive angle. First find two simulated times with the angle (in the second column) is negative on one line and positive on the line just below it. Delete all the lines before and after this event, but leave lines for about 5 time steps before and after the zero crossing. Finally use the editor to also delete the third column with the angular velocity data, and save the data file. Let the interpolate_solve tool interpolate an exact zero angle crossing time with:
./interpolate_solve pendulum.dat 0
This will display the inverse interpolation polynomial divided difference table and the interpolation of the time the polynomial crosses zero.
To find the period, naturally we need two zero crossing times to be able to subtract and find the elapsed time between them. So rerun the simulation with the ./pendulum command above to rewrite the pendulum.dat file, and open the rewritten data file in the gedit editor again. This time edit the lines down to the next to the last time the angle crosses zero from negative to positive, and again leaving only two columns for the interpolation tool. Rerun the interpolate_solve tool to converge on this time, and use a calculator to subtract the two times. Compare this simulated period with what you would calculate using the small angle approximation for a pendulum.
Now try this same procedure but use an angle that is closer to 90 degrees, such as 1 radian, for the initial angular displacement of the pendulum. What period do you measure? How much error is there in the small angle approximation for a pendulum that is swinging at larger angles?
In this section we will solve the second order differential equation that describes a damped harmonic oscillator with a sinusoidal forcing function. The solution of this differential equation is the damped sinusoid that characterizes spring, mass and damper mechanical systems, electrical systems with capacitance, inductance and resistance, and other systems.
First we will look at solutions with no forcing sinusoid, but a nonzero initial condition that will ring the system until it damps out. The solution of the second order differential equation for this case will typically look something like:
The red curve is the solution of the variable that the differential equation is expressed in, for instance the displacement of a mechanical system, and the green curve is the other state variable, the first derivative of the displacement, or the velocity.
Begin by opening the source code file damped_oscillator.c in a text editor:
gedit damped_oscillator.c &
Note the overall structure of this program is similar to the pendulum simulator in section 1. Follow the code in the main progam as it extracts eight parameters off of the command line, in order: the natural frequency of the system, the resonance quality factor 'Q', the forcing frequency and amplitude, the initial values for y0 and y1, and finally the stop and increment values of the independent variable x. Note where the rk4 function is called near the end of the main program to invoke the fourth order Runge-Kutta algorithm to solve the differential equation. The function whose pointer is passed to the rk4 routine is defined earlier in the source code file, and is named harmonic_osc(). Complete the two lines of code indicated so that the function models the harmonic oscillator.
Save your source code file and compile with:
gcc -Wall damped_oscillator.c diffeq.c -o damped_oscillator -lm
First try running your compiled executable program and viewing its help line output:
./damped_oscillator
To run a simulation of the oscillatory system with no forcing function, set the amplitude and frequency of the forcing sinusoid to zero. Set the natural frequency to 2 π and the Q factor to 10 for a fairly "ringy" system. Give it an initial amplitude of 1 and an initial velocity of 0 to initialize the system with some energy and to start the simulation at the peak of a cosine phase solution. Run the simulation up to 10 seconds, and with a fairly fine time step, say 10 msec. Assemble your command line and rerun the damped_oscillator program with its standard output redirected into a data file damped_oscillator.dat.
Now start another terminal window and run gnuplot. Plot the time functions of both the displacement and velocity together in one plot, noting that the data file written by the Runga-Kutta solver will be three columns, the first being the time value and the second and third being the two state variables. The gnuplot keyword 'using' included in the plot command will let you select which pairs of columns form the x and y coordinates for each plot curve. You should end up with a plot resembling that above.
We can see the effect of a forcing sinusoid at various frequencies using the other command line arguments. Assemble another command line that will set a forcing amplitude of 1 and at a frequency equal to the natural frequency of the resonator, 2 π. For this simulation you can set the initial amplitude and velocity both to zero, as the forcing function will provide energy to the system over time. Rerun the simulator and replot in gnuplot to see the resonant oscillations gradually building up over time. Now try different forcing frequencies above and below the resonant frequency of the system. What happens to the oscillation amplitude?
In this section we will solve the second order vector differential equation that describes a the motion of a ballistic object under the influence of quadratic air drag. This differential equation is solved in two dimensions with four state variables.
Begin by opening the source code file ballistic.c in a text editor:
gedit ballistic.c &
Note the overall structure of this program is similar to the pendulum simulator in section 1. Follow the code in the main progam as it extracts six parameters off of the command line, in order: the mass of the object in kg, the initial launch velocity in meters/second, the launch angle relative to the horizontal in degrees, the drag coefficient, and finally the stop and increment values of the independent variable t in seconds. Note where the rk4 function is called near the end of the main program to invoke the fourth order Runge-Kutta algorithm to solve the differential equation. The function whose pointer is passed to the rk4 routine is defined earlier in the source code file, and is named ballistic(). Complete the lines of code indicated in both the main program and the ballistic() function so that the state variables are properly initialized and the ballistic() function models the ballistic motion. Note that the computation of the global variable 'theta' includes the conversion factor DEGTORAD (from constants.h) that allows the launch angle to be entered more conveniently in degrees on the command line.
Save your source code file and compile with:
gcc -Wall ballistic.c diffeq.c -o ballistic -lm
First try running your compiled executable program and viewing its help line output:
./ballistic
Start simulations with a nominal case of a 1 kg mass launched at 100 m/sec at an angle of 45 degrees. For these initial simulations set the drag coefficient alpha to a small number such as 1.0e-6 to make the drag influence insignificant for now. Use a simulation time of 30 seconds and a time step of 10 msec. Assemble your command line and rerun the ballistic program with its standard output redirected into a data file ballistic.dat.
Now start another terminal window and run gnuplot. Plot the relationship betweeen the y and x spatial coordinates of the object, noting that the data file written by the Runga-Kutta solver will be five columns, the first being the time value and the second through the fifth being the four state variables. The gnuplot keyword 'using' included in the plot command will let you select which pairs of columns form the x and y coordinates for each plot curve. You should end up with a plot that gives the simple parabolic path expected with no air drag. An easy way to restrict the plot to the realistic path of the object is to set a fixed xrange to 0 to 1200 meters and the yrange to 0 to 500 meters.
Now try some simulations that introduce some air drag by increasing the alpha parameter to 1e-3 and replotting. The solution of the differential equation will be distorted from the ideal parabola. Assume that the launch angle remains fixed at 45 degrees, and find the increased velocity that the object must be launched with to reach a target distance of 1000 meters. Try different launch velocities until your gnuplot looks like:
A good way to hunt for the optimum launch velocity is to use a binary search similar to how the bisection algorithm searches for the root of a nonlinear function. (This process in fact forms the basis for the "shooting method" for solving boundary value problems numerically with iterative differential equation algorithms such as Runge-Kutta.)
If you are working on this lab exercise on a Linux system that has a script interpreter for Tcl/Tk installed, such as our classroom systems, verify your result by running the GUI tool:
./ballistic_gui
Enter your settings and results from above in the GUI sliders. Note: If you get an error message when trying to execute the script file, you may need to set the file's execute permission bits on by first running:
chmod a+x ballistic_gui