-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path12.RMS_Error_Midpoint-RK45_ODE.c
200 lines (158 loc) · 6.14 KB
/
12.RMS_Error_Midpoint-RK45_ODE.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
/****
This program find the RMS Error between the values obtained by solving the ODE using Midpoint Method and 4th Order Runge Kutta Method
****/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define gamma -1.76e11
#define Hz 7.96e4 //100e-3*(1/Mu0)
#define omega (gamma*(100e-3)*-1) //*-1 to get positive T0(time for one precession) omega = Mu0*Hz*gamma
#define T0 ((2*M_PI)/omega)
#define alpha 0.05
double step_size;
double lambda=0.0;
double* theta; //stores values for midpoint method
double* phi;
double* rktheta; //stores values for RK45 method
double* rkphi;
double tk1=0,tk2=0,pk1=0,pk2=0; //stores values for midpoint method
double rktk1=0,rktk2=0,rktk3=0,rktk4=0,rkpk1=0,rkpk2=0,rkpk3=0,rkpk4=0; //stores value for RK45 Method
double t_half_step=0.0; //half step in theta
double p_half_step=0.0; //halfstep in phi
double mod = 0.0;
double tvalue=0.0;
double pvalue=0.0;
double step=0.0;
int flg = 1;
double val = M_PI / 180.0;
int n=0; //to scale theta between 0 and 360
double temp=0.0; //double to turn negative angles into positive
//double error_theta_phi=0.0; //double to store error in angle if we want to find it.
double error_m=0.0;
double dthetadt(double theta1,double phi1)
{
tvalue=lambda*Hz*sin(theta1*val);
return tvalue;
}
double dphidt(double theta2,double phi2)
{
pvalue=dthetadt(theta2,phi2)*-1*(1/(alpha*sin(theta2*val)));
return pvalue;
}
int main(int argc,char* argv[])
{
step_size=atof(argv[1]);
lambda=alpha*(gamma/(1+pow(alpha,2)));
theta = malloc(1*(sizeof(double)));
phi = malloc(1*(sizeof(double)));
rktheta = malloc(1*(sizeof(double)));
rkphi = malloc(1*(sizeof(double)));
int theta_initial= 179 ; //in degrees
int phi_initial = 1 ; //in degrees
/* Using the relation:
mx[0]=(sin(theta[0]*val))*(cos(phi[0]*val));
my[0]=(sin(theta[0]*val))*(sin(phi[0]*val));
mz[0]=cos(theta[0]*val);
*/
theta[0]=theta_initial;
phi[0]=phi_initial;
rktheta[0]=theta_initial;
rkphi[0]=phi_initial;
//initial error is zero
//CALCULATING VALUES OF THETA AND PHI
long long unsigned int size=2;
int i=0;
while(flg>0)
{
theta=realloc(theta,(size)*sizeof(double));
phi=realloc(phi,(size)*sizeof(double));
rktheta=realloc(rktheta,(size)*sizeof(double));
rkphi=realloc(rkphi,(size)*sizeof(double));
//Using Midpoint Method Method
tk1=(dthetadt(theta[i],phi[i]))*step_size;
pk1=(dphidt(theta[i],phi[i]))*step_size;
t_half_step=(0.5)*tk1;
p_half_step=(0.5)*pk1;
tk2=(dthetadt(theta[i]+t_half_step, phi[i]+p_half_step))*step_size;
pk2=(dphidt(theta[i]+t_half_step, phi[i]+p_half_step))*step_size;
theta[i+1]=theta[i]+(tk2);
phi[i+1]=phi[i]+(pk2);
//Using RK45 or Fourth Order Runge-Kutta Method
rktk1=(dthetadt(rktheta[i],rkphi[i]))*step_size;
rkpk1=(dphidt(rktheta[i],rkphi[i]))*step_size;
t_half_step=(0.5)*tk1;
p_half_step=(0.5)*pk1;
rktk2=(dthetadt(rktheta[i]+t_half_step, rkphi[i]+p_half_step))*step_size;
rkpk2=(dphidt(rktheta[i]+t_half_step, rkphi[i]+p_half_step))*step_size;
t_half_step=(0.5)*tk2;
p_half_step=(0.5)*pk2;
rktk3=(dthetadt(rktheta[i]+t_half_step, rkphi[i]+p_half_step))*step_size;
rkpk3=(dphidt(rktheta[i]+t_half_step, rkphi[i]+p_half_step))*step_size;
rktk4=(dthetadt(rktheta[i]+rktk3, rkphi[i]+rkpk3))*step_size;
rkpk4=(dphidt(rktheta[i]+rktk3, rkphi[i]+rkpk3))*step_size;
rktheta[i+1]=rktheta[i]+((rktk1+(2*rktk2)+(2*rktk3)+rktk4)/6);
rkphi[i+1]=rkphi[i]+((rkpk1+(2*rkpk2)+(2*rkpk3)+rkpk4)/6);
//scaling values down of phi from 0 to 360
if(phi[i+1]>0&&phi[i+1]>360)
{
n=floor(phi[i+1]/360);
phi[i+1]=phi[i+1]-(n*360);
}
if(phi[i+1]<0&&fabs(phi[i+1]>360))
{
temp=fabs(phi[i+1]);
n=floor(temp/360);
temp=temp-(n*360);
phi[i+1]=360-temp;
}
//scaling values down of rkphi from 0 to 360
if(rkphi[i+1]>0&&rkphi[i+1]>360)
{
n=floor(rkphi[i+1]/360);
rkphi[i+1]=rkphi[i+1]-(n*360);
}
if(rkphi[i+1]<0&&fabs(rkphi[i+1]>360))
{
temp=fabs(rkphi[i+1]);
n=floor(temp/360);
temp=temp-(n*360);
rkphi[i+1]=360-temp;
}
/* Using the following relations to get <mx>,<my>,<mz>
mx[i+1]=(sin(theta[i+1]*val))*(cos(phi[i+1]*val));
my[i+1]=(sin(theta[i+1]*val))*(sin(phi[i+1]*val));
mz[i+1]=cos(theta[i+1]*val);
*/
//FINDING AND PRINTING ERROR
/* //error in theta,phi (if needed to find)
double error_theta = (rktheta[i+1]-theta[i+1]);
double error_phi = (rkphi[i+1]-phi[i+1]);
error_theta_phi = error_theta_phi + pow(error_theta,2);//ignoring error in phi since it is negligible
*/
//error in <mx>,<my>,<mz>
double error_mx = ((sin(rktheta[i+1]*val))*(cos(rkphi[i+1]*val)))-((sin(theta[i+1]*val))*(cos(phi[i+1]*val)));
double error_my = ((sin(rktheta[i+1]*val))*(sin(rkphi[i+1]*val)))-((sin(theta[i+1]*val))*(sin(phi[i+1]*val)));
double error_mz = ((cos(rktheta[i+1]*val))-(cos(theta[i+1]*val)));
error_m = error_m + (pow(error_mx,2)+pow(error_my,2)+pow(error_mz,2));
//time taken = (i+1)*step_size
//flag to stop calculating when dipole alligns with H i.e. when angle between M and H is less than 1 degree
if(theta[i+1]<=1)
{flg=0;}
if(theta[i+1]>1)
{
i++;
size++;
}
}
//PRINTING ERROR FOR GIVEN STEP SIZE
//printf("%le %le\n",step_size,error_theta_phi); //If needed to find error in angle
//printf("%le %le\n",step_size,error_m); //If needed to find the squared error in <m>
printf("%le %le\n",(step_size/pow(10,-15)),sqrt(error_m/size)); //prints RMS Error in <m>
//step size in femto seconds
//freeing memory
free(theta);
free(phi);
free(rktheta);
free(rkphi);
return 0;
}