Thread: Runge Kutta for electron trajectories

  1. #1
    Registered User
    Join Date
    May 2007

    Runge Kutta for electron trajectories

    Seriously desperate as deadline is the end of the week. I am having a complete nightmare with my code to do a 4th order Runge-Kutta to get trajectories for electron paths. I have written some code that at least outputs something, however if I try and vary ANYTHING it cocks up.

    Basically I want to change my lower limit to -0.003, my upper limit to 0.007.
    Initial conditions:

    My original equation is
    d^2y/dx^2 + ((e[B0exp(-x^2/0.000001)])/(8mV))
    where B0 is 1, 5 or 15


    #include <stdio.h>
    #include <math.h>
    #define N 2			/* number of first order equations */
    #define dist 0.1		/* stepsize in t*/
    #define MAX 10		/* max for t */
    #define charge 1.60217646e-19	/* elementary charge*/
    #define mass 9.1093897Eľ31	/* rest mass of electron */
    FILE *output;			/* internal filename */
    	double t, y[N];
    	int j;
    	void runge4(double x, double y[], double step);	/* Runge-Kutta function */
    	double f(double x, double y[], int i);		/* function for derivatives */
    	output=fopen("traject.dat", "w");			/* external filename */
    	y[0]=0.001;					/* initial position */
    	y[1]=0;							/* initial velocity */
    	fprintf(output, "-3\t%f\n", y[0]);
    	for (j=1; j*dist<=MAX ;j++)			/* time loop */
       runge4(t, y, dist);
       fprintf(output, "%f\t%f\n", t, y[0]);
    void runge4(double x, double y[], double step){
    	double h=step/2.0,			/* the midpoint */
    	t1[N], t2[N], t3[N],		/* temporary storage arrays */
    	k1[N], k2[N], k3[N],k4[N];	/* for Runge-Kutta */
    	int i;
    	for (i=0;i<N;i++) t1[i]=y[i]+0.5*(k1[i]=step*f(x, y, i));
    	for (i=0;i<N;i++) t2[i]=y[i]+0.5*(k2[i]=step*f(x+h, t1, i));
    	for (i=0;i<N;i++) t3[i]=y[i]+    (k3[i]=step*f(x+h, t2, i));
    	for (i=0;i<N;i++) k4[i]=		step*f(x+step, t3, i);
    	for (i=0;i<N;i++) y[i]+=(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
    double  f(double x, double y[], int i){
    	if (i==0) return(y[1]);						/* derivative of first equation */
    	if (i==1) return(-y[1]-y[0]);
    	//if (i==1) return(-(charge*(exp(pow(-t,2)/pow(-.001,2))/(18*1e5))*y[1]-y[0]);	
    		/* derivative of second equation */

  2. #2
    Mad OnionKnight's Avatar
    Join Date
    Jan 2005
    Umeň, Sweden
    Break your code up into more manageable pieces, not only is it easier to read but you'll be able to insert debugging statements in between, also having more than one assignment in a statement is kinda bad.
    Put a bunch of printf()s in your code to print the values of your variables as the function continues, then compare your output with manual calculations to see where it went wrong.

Popular pages Recent additions subscribe to a feed