[C language] Calculate the motion of charged particles in a uniform magnetic field by the Runge-Kutta method

Motion of charged particles in a uniform magnetic field

Equation of motion

m\frac{d \mathbf{v} }{dt} = q \mathbf{v} × \mathbf{B}

$ \ mathbf {B} $ is uniform and in the $ z $ direction, q> 0. Initial condition: $ t = 0 and x = y = 0, v_x = v_0, v_y = 0 $

$ x $ component

\frac{dv_x}{dt} = \frac{qB}{m}v_y

$ y component $

\frac{dv_y}{dt} = -\frac{qB}{m}v_x

To make it easier to calculate, replace it as follows.

u_x = \frac{v_x}{v_0}
\\
u_y = \frac{v_y}{v_0}
\\
τ = \frac{qBt}{m}
\\
X = \frac{qBx}{mv_0}
\\
Y = \frac{qBy}{mv_0}

Using these, the calculation is as follows.

\frac{du_x}{dτ} = u_y \・ ・ ・ ①
\\
\frac{du_y}{dτ} = -u_x \・ ・ ・ ②
\\
\frac{dX}{dτ} = u_x
\\
\frac{dY}{dτ} = u_y

Initial condition: $ τ = 0 and X = Y = 0, u_x = 1, u_y = 0 $

sample

The motion of a positively charged charged particle emitted from the origin in a uniform magnetic field in the $ x $ direction at a velocity of $ v_0 $ is calculated by the Lungekutter method and illustrated.

sample.c

#include <stdio.h>
#include <math.h>

int main(void) {

	FILE *data;

	double x1, x2, y1, y2, vx1, vx2, vy1, vy2, t, dt, x, y, kx1, kx2, kx3, kx4, kx;
	double ky1, ky2, ky3, ky4, ky, hx1, hx2, hx3, hx4, hx, hy1, hy2, hy3, hy4, hy, x0, y0;

	data = fopen("data-sample.csv", "w");

	x1 = 0;
	y1 = 0;
	vx1 = 1;
	vy1 = 0;
	dt = 0.01;

	for (int i=0; i<=1000; i++) {
        t = i*dt;

        kx1 = dt*vx1;
        ky1 = dt*vy1;
        hx1 = dt*vy1;
        hy1 = -dt*vx1;

        kx2 = dt*(vx1 + hx1/2.);
        ky2 = dt*(vy1 + hy1/2.);
        hx2 = dt*(vy1 + hy1/2.);
        hy2 = dt*( -(vx1 + hx1/2.) );

        kx3 = dt*(vx1 + hx2/2.);
        ky3 = dt*(vy1 + hy2/2.);
        hx3 = dt*(vy1 + hy2/2.);
        hy3 = dt*( -(vx1 + hx2/2.) );

        kx4 = dt*(vx1 + hx3);
        ky4 = dt*(vy1 + hy3);
        hx4 = dt*(vy1 + hy3);
        hy4 = dt*( -(vx1 + hx3) );

        kx = (kx1 + 2*kx2 + 2*kx3 + kx4)/6.;
        ky = (ky1 + 2*ky2 + 2*ky3 + ky4)/6.;

        x2 = x1 + kx;
        y2 = y1 + ky;

        hx = (hx1 + 2*hx2 + 2*hx3 + hx4)/6.;
        hy = (hy1 + 2*hy2 + 2*hy3 + hy4)/6.;
        vx2 = vx1 + hx;
        vy2 = vy1 + hy;

        x0 = sin(t);
        y0 = cos(t) - 1;

        if (fmod(i, 10)<0.01) {
            printf("t = %5.1f,Calculated value-> x=%f, y=%f,Theoretical value-> x0=%f, y0=%f\n", t, x1, y1, x0, y0);
            fprintf(data,"%f, %f, %f, %f\n", x1, y1, x0, y0);
        }

        x1 = x2;
        y1 = y2;
        vx1 = vx2;
        vy1 = vy2;

    }
    fclose(data);

    return 0;
}

Execution result

problem

The motion of a negatively charged charged particle emitted from the origin in a uniform magnetic field in the $ y $ direction at a velocity of $ v_0 $ is calculated by the Lungekutter method and illustrated.

---

Since it is a negative charge, the signs of equations (1) and (2) are interchanged.

problem.c

#include <stdio.h>
#include <math.h>

int main(void) {

	FILE *data;

	double x1, x2, y1, y2, vx1, vx2, vy1, vy2, t, dt, x, y, kx1, kx2, kx3, kx4, kx;
	double ky1, ky2, ky3, ky4, ky, hx1, hx2, hx3, hx4, hx, hy1, hy2, hy3, hy4, hy, x0, y0;

	data = fopen("data-problem.csv", "w");

	x1 = 0;
	y1 = 0;
	vx1 = 0;
	vy1 = 1;
	dt = 0.01;

	for (int i=0; i<=1000; i++) {
        t = i*dt;

        kx1 = dt*vx1;
        ky1 = dt*vy1;
        hx1 = -dt*vy1;
        hy1 = dt*vx1;

        kx2 = dt*(vx1 + hx1/2.);
        ky2 = dt*(vy1 + hy1/2.);
        hx2 = dt*( -(vy1 + hy1/2.) );
        hy2 = dt*(vx1 + hx1/2.);

        kx3 = dt*(vx1 + hx2/2.);
        ky3 = dt*(vy1 + hy2/2.);
        hx3 = dt*(-(vy1 + hy2/2.) );
        hy3 = dt*(vx1 + hx2/2.);

        kx4 = dt*(vx1 + hx3);
        ky4 = dt*(vy1 + hy3);
        hx4 = dt*( -(vy1 + hy3) );
        hy4 = dt*(vx1 + hx3);

        kx = (kx1 + 2*kx2 + 2*kx3 + kx4)/6.;
        ky = (ky1 + 2*ky2 + 2*ky3 + ky4)/6.;
        x2 = x1 + kx;
        y2 = y1 + ky;

        hx = (hx1 + 2*hx2 + 2*hx3 + hx4)/6.;
        hy = (hy1 + 2*hy2 + 2*hy3 + hy4)/6.;
        vx2 = vx1 + hx;
        vy2 = vy1 + hy;

        x0 = sin(t) - 1;
        y0 = cos(t);

        if (fmod(i, 10)<0.01) {
            printf("t = %5.1f,Calculated value-> x=%f, y=%f,Theoretical value-> x0=%f, y0=%f\n", t, x1, y1, x0, y0);
            fprintf(data,"%f, %f, %f, %f\n", x1, y1, x0, y0);
        }

        x1 = x2;
        y1 = y2;
        vx1 = vx2;
        vy1 = vy2;

    }
    fclose(data);

    return 0;
}

Execution result