// ******************************************

// ** Soution of coupled differential **

// ** equations (simple pendulum) using **

// ** 2nd order Runge-Kutta method **

// ** COMPUTING II (F31SC6) PJM 11/11/97 **

// ******************************************

#define g 9.8

#include <stdio.h>

#include <math.h>

float d_omega_dt(float);

float Calculate_Energy(float, float);

float time_step=0.01; // Time step in Runge-Kutta method

float l=9.8; // length of pendulum

int main(void)

{

FILE *fp_position, *fp_energy;

float omega[5000], theta[5000],time=0.0,energy,k1,theta_half,omega_half;

int t;

theta[0]=M_PI/18.0; // Starting angle is 10 degrees

omega[0]=0.00; // Starting velocity is 0 rad/s

fp_position=fopen("position.dat","w");

fp_energy=fopen("energy.dat","w");

for (t=0;t<4999;t++)

{

// *****************************

// ** 2nd order Runge-Kutta **

// ** method is given in next **

// ** four lines of code **

// *****************************

// First, evaluate values of theta and omega at centre of

// interval (using Euler's method with delta_t replaced

// by delta_t/2)............

theta_half=theta[t]+((time_step)/2.0)*omega[t];

omega_half=omega[t]+(d_omega_dt(theta[t]))*time_step/2.0;

// Then, calculate new values of theta and omega (i.e.

// values of theta and omega at t+delta_t)........

omega[t+1]=omega[t]+d_omega_dt(theta_half)*time_step;

theta[t+1]=theta[t]+(time_step*omega_half);

// Calculate energy

energy=Calculate_Energy(theta[t],omega[t]);

// Write values of energy and position to a file (6 dp)

fprintf(fp_position, "%f,%.6f \n",time,theta[t]);

fprintf(fp_energy, "%f,%.6f \n",time,energy);

time+=time_step;

}

fclose(fp_position);

fclose(fp_energy);

}

// Function to calculate dv/dt using

// formula dv/dt = -(g/l)*(theta)

// (Small angle approximation)

float d_omega_dt(float angle)

{

float dveldt;

dveldt=-(g/l)*(angle);

return(dveldt);

}

float Calculate_Energy(float angle, float vel)

{

float total_energy,kinetic_energy,pot_energy;

float lin_vel;

lin_vel=l*vel; // v=rw

kinetic_energy=0.5*lin_vel*lin_vel; // mass = 1 kg

// You should be able to show that the potential energy is

// given by U=(mgl/2)*(theta*theta)..........

pot_energy=(g*l/2.0)*angle*angle;

total_energy=kinetic_energy+pot_energy;

return (total_energy);

}