序言
过了好几天,终于有时间把这个写出来啦!!!
接着上一篇文章继续写。
Precondition
the situation is special,the initial condition is shown in the following figure
Part2
Motion Of Three Particles
E = 0.3; %决定粒子往x轴运动的情况
B = 0.65; %决定粒子往y轴旋转的情况
kappa = 9.6;
const = 0.125;
dt = 0.01;
x1 = zeros(1,100);
y1 = zeros(1,100);
z1 = zeros(1,100);
x2 = zeros(1,100);
y2 = zeros(1,100);
z2 = zeros(1,100);
x3 = zeros(1,100);
y3 = zeros(1,100);
z3 = zeros(1,100);
vx1 = zeros(1,100);
vy1 = zeros(1,100);
vz1 = zeros(1,100);
vx2 = zeros(1,100);
vy2 = zeros(1,100);
vz2 = zeros(1,100);
vx3 = zeros(1,100);
vy3 = zeros(1,100);
vz3 = zeros(1,100);
t = zeros(1,100);
x2(1) = 1;
y3(1) = 1;
i = 1;
rrr12 = zeros(1,100);
rrr13 = zeros(1,100);
rrr23 = zeros(1,100);
rrr12(1) = 1;
rrr13(1) = 1;
rrr23(1) = 2^0.5;
while i<round(10/dt)
vx1(i+1) = vx1(i) + dt*kappa*(E+const*((x1(i)-x2(i))/rrr12(i) + (x1(i)-x3(i))/rrr13(i))-B*vz1(i));
vy1(i+1) = vy1(i) + dt*kappa*const*((y1(i)-y2(i))/rrr12(i) + (y1(i)-y3(i))/rrr13(i));
vz1(i+1) = vz1(i) + dt*kappa*(const*((z1(i)-z2(i))/rrr12(i) + (z1(i)-z3(i))/rrr13(i))+B*vx1(i));
x1(i+1) = x1(i) + dt*vx1(i+1);
y1(i+1) = y1(i) + dt*vy1(i+1);
z1(i+1) = z1(i) + dt*vz1(i+1);
vx2(i+1) = vx2(i) + dt*kappa*(E+const*((x2(i)-x1(i))/rrr12(i) + (x2(i)-x3(i))/rrr13(i))-B*vz2(i));
vy2(i+1) = vy2(i) + dt*kappa*const*((y2(i)-y1(i))/rrr12(i) + (y2(i)-y3(i))/rrr13(i));
vz2(i+1) = vz2(i) + dt*kappa*(const*((z2(i)-z1(i))/rrr12(i) + (z2(i)-z3(i))/rrr13(i))+B*vx2(i));
x2(i+1) = x2(i) + dt*vx2(i+1);
y2(i+1) = y2(i) + dt*vy2(i+1);
z2(i+1) = z2(i) + dt*vz2(i+1);
vx3(i+1) = vx3(i) + dt*kappa*(E+const*((x3(i)-x1(i))/rrr12(i) + (x3(i)-x2(i))/rrr13(i))-B*vz3(i));
vy3(i+1) = vy3(i) + dt*kappa*const*((y3(i)-y1(i))/rrr12(i) + (y3(i)-y2(i))/rrr13(i));
vz3(i+1) = vz3(i) + dt*kappa*(const*((z3(i)-z1(i))/rrr12(i) + (z3(i)-z2(i))/rrr13(i))+B*vx3(i));
x3(i+1) = x3(i) + dt*vx3(i+1);
y3(i+1) = y3(i) + dt*vy3(i+1);
z3(i+1) = z3(i) + dt*vz3(i+1);
rrr12(i+1) = ((x1(i+1)-x2(i+1))^2+(y1(i+1)-y2(i+1))^2+(z1(i+1)-z2(i+1))^2)^1.5;
rrr13(i+1) = ((x1(i+1)-x3(i+1))^2+(y1(i+1)-y3(i+1))^2+(z1(i+1)-z3(i+1))^2)^1.5;
rrr23(i+1) = ((x3(i+1)-x2(i+1))^2+(y3(i+1)-y2(i+1))^2+(z3(i+1)-z2(i+1))^2)^1.5;
t(i+1) = t(i) + dt;
i = i + 1;
end
粒子的运动轨迹如下图
说明一下此图仅描述了粒子一的运动情况
总结
有了初步建模的方法,以后会再接再励模拟更复杂更有趣的系统,为兴趣而活,模拟全世界!
还有数据分析也是自己的一个发展方向。日后也要多家练习。