my problem is, frankly, i'm unsure how works.
i need modify double f() function solve arbitrary differential equation d2θ/dt2 = −ω2sinθ, unsure how proceed.
the rk4 function runge4() understand; don't understand how f() function returns correct values harmonic oscillator.
would please @ least explain logic behind f() function?
original code below.
/* ************************************************************************ * rk4.c: 4th order runge-kutta solution harmonic oscillator * * * * from: "a survey of computational physics" rh landau, mj paez, , cc bordeianu copyright princeton university press, princeton, 2008. electronic materials copyright: r landau, oregon state univ, 2008; mj paez, univ antioquia, 2008; & cc bordeianu, univ bucharest, 2008 support national science foundation * ************************************************************************ */ #include <stdio.h> #define n 2 /* number of equations */ #define dist 0.1 /* stepsize */ #define min 0.0 /* minimum x */ #define max 10.0 /* maximum x */ void runge4(double x, double y[], double step); double f(double x, double y[], int i); int main() { double x, y[n]; int j; file *output; /* save data in rk4.dat */ output = fopen("rk4.dat","w"); y[0] = 1.0; /* initial position */ y[1] = 0.0; /* initial velocity */ fprintf(output, "%f\t%f\n", x, y[0]); for(x = min; x <= max ; x += dist) { runge4(x, y, dist); fprintf(output, "%f\t%f\n", x, y[0]); /* position vs. time */ } printf("data stored in rk4.dat\n"); fclose(output); } /*-----------------------end of main program--------------------------*/ /* runge-kutta subroutine */ void runge4(double x, double y[], double step) { double h=step/2.0, /* midpoint */ t1[n], t2[n], t3[n], /* temporary storage */ k1[n], k2[n], k3[n],k4[n]; /* runge-kutta */ int i; (i=0; i<n; i++) t1[i] = y[i]+0.5*(k1[i]=step*f(x, y, i)); (i=0; i<n; i++) t2[i] = y[i]+0.5*(k2[i]=step*f(x+h, t1, i)); (i=0; i<n; i++) t3[i] = y[i]+ (k3[i]=step*f(x+h, t2, i)); (i=0; i<n; i++) k4[i] = step*f(x + step, t3, i); (i=0; i<n; i++) y[i] += (k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0; } /*--------------------------------------------------------------------*/ /* definition of equations - harmonic oscillator */ double f(double x, double y[], int i) { if (i == 0) return(y[1]); /* rhs of first equation */ if (i == 1) return(-y[0]); /* rhs of second equation */ }
start hooke's law:
f = -kx combine newton's second law differential equation linear harmonic oscillator:
ma = f = -kx mx'' = -kx x'' = -k/m x arbitrarily chose our units k/m == 1, , equation becomes just:
x'' = -x now, introduce dummy variable y = x', , write second-order differential equation two-dimensional first-order system:
x' = y y' = -x the function f in code encodes system; i'm going change variable names clarity:
double f(double t, double v[], int i) { if (i == 0) return(v[1]); if (i == 1) return(-v[0]); } v vector [x,y] 2 dimensional system above. given i, t, , v, function f returns derivative respect t of ith component of v. re-writing 2d system using these names, get:
dv[0]/dt = v[1] dv[1]/dt = -v[0] which function f does.
Comments
Post a Comment