Not sure if this is the right place: Runge Kutta in C for electron trajectory

  • Thread starter QueenB
  • Start date
  • #1
1
0

Homework Statement


Code to do a 4th order Runge-Kutta to get trajectories for electron paths.
When trying to substitute conditions into my program, it returns rubbish.

Homework Equations



Basically I want to change my lower limit to -0.003, my upper limit to 0.007.
Initial conditions:
y(x=-0.003)=a
dy/dx|(x=-0.003)=0

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

The Attempt at a Solution



Code:
#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 */

main(){
	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 */
	{
   t=(j*dist)-3;
   runge4(t, y, dist);

   fprintf(output, "%f\t%f\n", t, y[0]);
	}

	fclose(output);
}

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 */
}
 

Answers and Replies

Related Threads on Not sure if this is the right place: Runge Kutta in C for electron trajectory

  • Last Post
Replies
4
Views
6K
Replies
3
Views
5K
  • Last Post
Replies
4
Views
7K
  • Last Post
Replies
2
Views
1K
Replies
3
Views
307
Replies
3
Views
651
Top