Edited at

[C言語] 一様磁場中の荷電粒子の運動をルンゲクッタ法で計算

More than 3 years have passed since last update.


一様磁場中の荷電粒子の運動

運動方程式

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

$\mathbf{B}$は一様で$z$方向、q>0とする。

初期条件:$t = 0 で x = y = 0, v_x = v_0, v_y = 0$

$x$成分

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

$y成分$

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

計算しやすいように、以下のように置き換える。

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}

これらを使って、計算すると以下のようになる。

\frac{du_x}{dτ} = u_y \ ・・・①

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

初期条件:$τ = 0 で X = Y = 0, u_x =1, u_y = 0$


サンプル

一様磁場中の原点から$x$方向へ速度$v_0$で放たれた正の電荷を持つ荷電粒子の運動をルンゲクッタ法で計算し、図示する。


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, 計算値 -> x=%f, y=%f, 理論値 -> 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;
}



実行結果

sample-result.png

sample-graph.png


問題

一様磁場中の原点から$y$方向へ速度$v_0$で放たれた負の電荷を持つ荷電粒子の運動をルンゲクッタ法で計算し、図示する。


負の電荷なので、式①と②の符号が入れ替わっている。


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, 計算値 -> x=%f, y=%f, 理論値 -> 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;
}



実行結果

problem-result.png

problem-graph.png