고급물리학 19-3-4, 29-3-2번 전산모사코드(RK4사용)

19-3-4.cpp

#include<stdio.h>
#include<math.h>
#define h 0.04
double a = 1; // alpha : (eB_0)/m
double e = 0.6666; // collison coefficient
double f1(double x,double y,double xv,double yv){ // dx/dt
	return xv;
}
double f2(double x,double y,double xv,double yv){ // dy/dt
	return yv;
}
double f3(double x,double y,double xv,double yv){ // d(xv)/dt
	return a*yv;
}
double f4(double x,double y,double xv,double yv){ // d(yv)/dt
	return -a*xv;
}
int main()
{
	freopen("output.csv","w",stdout);
	double x,y,xv,yv;
	x=0;
	y=1;
	xv=0.8;
	yv=0;
	for(int t=0;t<2000;t++){
		printf("%lf,%lf,%lf,\n",h*t,x,y);
		double k11,k12,k13,k14,k21,k22,k23,k24,k31,k32,k33,k34,k41,k42,k43,k44;
		
		k11=h*f1(x,y,xv,yv);
		k12=h*f2(x,y,xv,yv);
		k13=h*f3(x,y,xv,yv);
		k14=h*f4(x,y,xv,yv);
		
		k21=h*f1(x+0.5*k11,y+0.5*k12,xv+0.5*k13,yv+0.5*k14);
		k22=h*f2(x+0.5*k11,y+0.5*k12,xv+0.5*k13,yv+0.5*k14);
		k23=h*f3(x+0.5*k11,y+0.5*k12,xv+0.5*k13,yv+0.5*k14);
		k24=h*f4(x+0.5*k11,y+0.5*k12,xv+0.5*k13,yv+0.5*k14);
		
		k31=h*f1(x+0.5*k21,y+0.5*k22,xv+0.5*k23,yv+0.5*k24);
		k32=h*f2(x+0.5*k21,y+0.5*k22,xv+0.5*k23,yv+0.5*k24);
		k33=h*f3(x+0.5*k21,y+0.5*k22,xv+0.5*k23,yv+0.5*k24);
		k34=h*f4(x+0.5*k21,y+0.5*k22,xv+0.5*k23,yv+0.5*k24);
		
		k41=h*f1(x+k31,y+k32,xv+k33,yv+k34);
		k42=h*f2(x+k31,y+k32,xv+k33,yv+k34);
		k43=h*f3(x+k31,y+k32,xv+k33,yv+k34);
		k44=h*f4(x+k31,y+k32,xv+k33,yv+k34);
		
		x +=(1.0/6.0)*(k11+2.0*k21+2.0*k31+k41);
		y +=(1.0/6.0)*(k12+2.0*k22+2.0*k32+k42);
		xv+=(1.0/6.0)*(k13+2.0*k23+2.0*k33+k43);
		yv+=(1.0/6.0)*(k14+2.0*k24+2.0*k34+k44);
		if(y<0){
			y=0;
			//y*=(-1);
			yv*=(-e);
		}
	}
}

29-3-2.cpp

#include<stdio.h>
#include<math.h>
#define h 0.1
double a = 1; // alpha : B=B_0(1+alpha*y)
double b = 1.00; // beta : (eB_0)/m
double f1(double x,double y,double xv,double yv){ // dx/dt
	return xv;
}
double f2(double x,double y,double xv,double yv){ // dy/dt
	return yv;
}
double f3(double x,double y,double xv,double yv){ // d(xv)/dt
	return -b*(1+a*y)*yv;
}
double f4(double x,double y,double xv,double yv){ // d(yv)/dt
	return b*(1+a*y)*xv;
}
int main()
{
	freopen("output.csv","w",stdout); // writes data on csv(excel) file. 
	double x,y,xv,yv; // x,y,v_x,v_y
	x=0; // Initial values
	y=0;
	xv=1;
	yv=0;
	for(int t=0;t<1000;t++){ // Runge-Kutta 4th method
		printf("%lf,%lf,%lf,\n",h*t,x,y);
		double k11,k12,k13,k14,k21,k22,k23,k24,k31,k32,k33,k34,k41,k42,k43,k44;
		
		k11=h*f1(x,y,xv,yv);
		k12=h*f2(x,y,xv,yv);
		k13=h*f3(x,y,xv,yv);
		k14=h*f4(x,y,xv,yv);
		
		k21=h*f1(x+0.5*k11,y+0.5*k12,xv+0.5*k13,yv+0.5*k14);
		k22=h*f2(x+0.5*k11,y+0.5*k12,xv+0.5*k13,yv+0.5*k14);
		k23=h*f3(x+0.5*k11,y+0.5*k12,xv+0.5*k13,yv+0.5*k14);
		k24=h*f4(x+0.5*k11,y+0.5*k12,xv+0.5*k13,yv+0.5*k14);
		
		k31=h*f1(x+0.5*k21,y+0.5*k22,xv+0.5*k23,yv+0.5*k24);
		k32=h*f2(x+0.5*k21,y+0.5*k22,xv+0.5*k23,yv+0.5*k24);
		k33=h*f3(x+0.5*k21,y+0.5*k22,xv+0.5*k23,yv+0.5*k24);
		k34=h*f4(x+0.5*k21,y+0.5*k22,xv+0.5*k23,yv+0.5*k24);
		
		k41=h*f1(x+k31,y+k32,xv+k33,yv+k34);
		k42=h*f2(x+k31,y+k32,xv+k33,yv+k34);
		k43=h*f3(x+k31,y+k32,xv+k33,yv+k34);
		k44=h*f4(x+k31,y+k32,xv+k33,yv+k34);
		
		x +=(1.0/6.0)*(k11+2.0*k21+2.0*k31+k41);
		y +=(1.0/6.0)*(k12+2.0*k22+2.0*k32+k42);
		xv+=(1.0/6.0)*(k13+2.0*k23+2.0*k33+k43);
		yv+=(1.0/6.0)*(k14+2.0*k24+2.0*k34+k44);
	}
}