Не получаются правильные расчеты
Задача следующая:
10 частиц распределены вдоль линии по Z, т.е. по направлению движения
Расстояние между частицами одинаково.
Q = 10 нКл – заряд каждой частицы.
Масса каждой частицы равна массе электрона.
L = 0,004 см – длина последовательности частиц.
N= 10 - кол-во частиц,
W = 10 МеВ – начальная полная энергия частиц.
Частицы начинают двигаться вдоль положительного направления оси Z
Найти распределение энергии между частицами после прохождения 1 см пути, учитывая только
Кулоновское взаимодействие.*
Для решения задачи я написал следующий код в Matlab (версия R2020b)
clc
clear
N=10;
dt=1*10^-12;
q=10^-9;
eps0=8.85*10^-12;
c=3*10^8;
pi=3.1415926;
m=9.1*10^-31;
k=1/(4*pi*eps0);
W=[1.6*10^-13 1.6*10^-13 1.6*10^-13 1.6*10^-13 1.6*10^-13 1.6*10^-13 1.6*10^-13 1.6*10^-13 1.6*10^-13 1.6*10^-13];
R=[0 0 0 0 0 0 0 0 0 0;0 0 0 0 0 0 0 0 0 0; 4*10^-5 8*10^-5 12*10^-5 16*10^-5 20*10^-5 24*10^-5 28*10^-5 32*10^-5 36*10^-5 4*10^-4];
dp=zeros(3,N);
dr=zeros(3,N);
V=zeros(3,N);
F_Q_2=zeros(3,N,N);
F_Q=zeros(3,N);
F_Q2=zeros(3,N);
Z=zeros(1,N);
p=zeros(3,N);
X=[];
Y=[];
Z=[];
R=vpa(R);
while R(3,10)<4*10^-4+10^-2
F_Q = ForceQ(F_Q, R,N,k,q);
for i=1:N
p(3,i)=sqrt((W(1,i)^2/c^2)-(m^2*c^2));
V(3,i)=((c^2/W(1,i))*p(3,i));
dp=F_Q(1:3,i)*dt;
p(1:3,i)=p(1:3,i)+dp;
pn=sqrt(p(1,i).^2+p(2,i).^2+p(3,i).^2);
W(1,i)=c*sqrt(pn^2+m^2*c^2);
V(1:3,i)=(c^2/W(1,i))*p(1:3,i);
dr(1:3,i)=V(1:3,i).*dt;
R(1:3,i)=R(1:3,i)+dr(1:3,i);
R = real(R);
Z(1,i) = R(3,i);
end
end
colormap hot
plot(Z(1:N),W(1:N), '.');
ylabel('\itW, Дж','fontsize',14)
xlabel('\itZ, м','fontsize',14)
function F_Q = ForceQ(F_Q2, R,N,k,q)
R = circshift(R,1,2);
for count=1:N
F_Q_2=zeros(3,N,N);
R = circshift(R,-1,2);
for var=2:N
r=sqrt((R(1,var)-R(1,1))^2+(R(2,var)-R(2,1))^2+(R(3,var)-R(3,1))^2);
force=(k*q^2)/r^2;
alpha=asin(real(R(3,var)-R(3,1))/r);
rxy=r*cos(alpha);
gamma=asin(real(R(2,var)-R(2,1))/rxy);
gamma2=asin(real(R(1,var)-R(1,1))/rxy);
F_Q_2(1,1,var)=F_Q2(1,count)-(force*1*sin(gamma2));
F_Q_2(2,1,var)=F_Q2(2,count)-(force*1*sin(gamma));
F_Q_2(3,1,var)=F_Q2(3,count)-(force*sin(alpha));
end
F_Q1 = sum(F_Q_2,3);
F_Q(:,count) = F_Q1(:,1);
end
end
Опуская физику, чисто математически мне нужно построить зависимость W от Z, причем частицы должны менять координаты таким образом, что первые 5 двигаются в отрицательном направлении, оставшиеся 5 - в положительном. В корректности формул я уверен, однако расчеты получаются не совсем верными. Проверял в debug, первые два прохода в цикле while считаются корректно, на 3 проход уже при вызове функции F_Q на выходе значения численно верные, но знаки у 5 и 6 частиц инверсно меняются (должно быть -число у 5 и +число у 6, а на выходе наоборот). Причем функцию F_Q проверял в отдельном файле с циклом for, считается корректно. Соответственно, из-за этого портятся все остальные расчеты координат и на выходе к концу цикла while получается совершенно не то, что должно быть. Пробовал также использовать другой вариант расчета F_Q без использования функции, через циклы, однако там проблема встает при расчете проекций вектора. Код вот (прикрепляю уже сам цикл, тк начальные условия до цикла те же):
while R(3,10)<4*10^-2
for i=1:N
p(3,i)=sqrt((W(1,i)^2/c^2)-(m^2*c^2));
V(3,i)=((c^2/W(1,i))*p(3,i));
for var=i+1:N
r=sqrt((R(1,var)-R(1,i))^2+(R(2,var)-R(2,i))^2+(R(3,var)-R(3,i))^2);
force=(k*q^2)/r^2;
alpha=asin(abs(R(3,var)-R(3,i))/r);
rxy=r*cos(alpha);
gamma=asin(abs(R(2,var)-R(2,i))/rxy);
F_Qx(1,i)=F_Qx(1,i)+(force*cos(alpha)*cos(gamma));
F_Qy(1,i)=F_Qy(1,i)+(force*cos(alpha)*sin(gamma));
F_Qz(1,i)=F_Qz(1,i)+(force*cos(alpha));
F_Qx(1,var)=F_Qx(1,var)-(force*cos(alpha)*cos(gamma));
F_Qy(1,var)=F_Qy(1,var)-(force*cos(alpha)*sin(gamma));
F_Qz(1,var)=F_Qz(1,var)-(force*cos(alpha));
F_Q(1:3,i)=[F_Qx(1,i); F_Qy(1,i); F_Qz(1,i)];
F_Q(1:3,var)=[F_Qx(1,var); F_Qy(1,var); F_Qz(1,var)];
end
end
for i=1:N
F=F_Q(1:3,i);
dp=F*dt;
p(1:3,i)=p(1:3,i)+dp;
pn=sqrt(p(1,i).^2+p(2,i).^2+p(3,i).^2);
W(1,i)=c*sqrt(pn^2+m^2*c^2);
V(1:3,i)=(c^2/W(1,i))*p(1:3,i);
dr(1:3,i)=V(1:3,i).*dt;
R(1:3,i)=R(1:3,i)+dr(1:3,i);
R = real(R);
Z = R(1:N);
end
end
colormap hot
plot(Z(1:N),W(1:N), '.');
ylabel('\itW, Дж','fontsize',14)
xlabel('\itZ, м','fontsize',14)
Я не понимаю в чем причина неправильных расчетов, мои предположения в том, что либо я неправильно оформил порядок в цикле, либо расчет портит circshift (эту функцию я использовал в другой похожей программе, расчет F_Q велся аналогично, и в целом там все работало корректно).