matlab模拟

matlab模拟

For today’s recreational coding exercise, we will look at the gravitational N-body problem. We will create a simulation of a dynamical system of particles interacting with each other gravitationally. Such a system may describe the orbits of planets in the Solar System or stars in our Galaxy.

对于今天的娱乐编码练习，我们将研究重力N体问题。 我们将创建一个粒子相互重力相互作用的动力学系统的模拟。 这样的系统可以描述太阳系中的行星或银河系中的恒星的轨道。

You may find the accompanying Matlab code on github. (And if you prefer to use Python, please see my Python version of the article)

您可以在github上找到随附的Matlab代码。 (如果您更喜欢使用Python，请参阅本文的Python版本)

But first, below is a gif of what running our simulation looks like:

但是首先，下面是模拟运行情况的图像：

力计算 (Force Calculation)

We will assume a system of N point particles, indexed by i=1,…,N. Each particle has a:

我们将假定一个由i = 1，…，N索引的N点粒子系统。 每个粒子都有一个：

mass mᵢ,

质量mᵢ

position rᵢ = [ xᵢ, yᵢ, zᵢ ] ,

位置rᵢ = [ xᵢ ， yᵢ ， zᵢ ]，

velocity vᵢ = [ vxᵢ, vyᵢ, vzᵢ ]

速度vᵢ= [VXᵢ，Vyᵢ，Vzᵢ]

Each particle feels the gravitational attraction of all the other particles according to Newton’s law of universal gravitation (the famous ‘inverse square-law’). That is, each particle feels an acceleration:

根据牛顿的万有引力定律(著名的“反平方律”)，每个粒子都感受到所有其他粒子的引力吸引。 也就是说，每个粒子都会感觉到加速度：

where G=6.67×10^-11 m³/kg/s² is the Gravitational constant.

其中G = 6.67×10 ^ -11m³/ kg /s²是引力常数。

We can write a Matlab function to perform the calculation on an input N×3 matrix of positions:

我们可以编写一个Matlab函数来对输入的N ×3位​​置矩阵进行计算：

function [ a ] = getAcc( pos, mass, G, softening ) %GETACC Calculate the acceleration on each particle due to Newton's Law % pos is an N x 3 matrix of positions % mass is an N x 1 vector of masses % G is Newton's Gravitational constant % softening is the softening length % a is N x 3 matrix of accelerations N = size(pos,1); a = zeros(N,3); for i = 1:N for j = 1:N if i ~= j dx = pos(j,1) - pos(i,1); dy = pos(j,2) - pos(i,2); dz = pos(j,3) - pos(i,3); inv_r3 = (dx^2 + dy^2 + dz^2 + softening^2)^(-1.5); a(i,1) = a(i,1) + G * (dx * inv_r3) * mass(j); a(i,2) = a(i,2) + G * (dy * inv_r3) * mass(j); a(i,3) = a(i,3) + G * (dz * inv_r3) * mass(j); end end end end

The softening parameter in the code is a small number added to avoid numerical issues when 2 particles are close to each other, in which case the acceleration from the ‘inverse square-law’ goes to infinity. In real life, masses are not exactly point particles and have finite extent.

代码中的软化参数是添加的一个小数字，以避免2个粒子彼此靠近时的数值问题，在这种情况下，来自“反平方律”的加速度将变为无穷大。 在现实生活中，质量并不完全是点粒子，并且范围有限。

The performance of the above code can actually be improved in Matlab by vectorization. That is, formulating the problem in terms of vector and matrix operations. I highly recommend this approach in general with interpreted languages, as it often can lead to 100× speedup. It also makes the code more readable. The downside is that storing intermediate calculations inside matrices uses significant memory. Here is a vectorized version of the function that computes the acceleration on all the particles:

在Matlab中通过矢量化实际上可以改善上述代码的性能。 也就是说，用向量和矩阵运算来表达问题。 我强烈建议在解释型语言中使用这种方法，因为它通常会导致100倍的加速。 这也使代码更具可读性。 缺点是将中间计算结果存储在矩阵中会占用大量内存。 这是该函数的矢量化版本，可计算所有粒子上的加速度：

function [ a ] = getAcc_vec( pos, mass, G, softening ) %GETACC Calculate the acceleration on each particle due to Newton's Law % pos is an N x 3 matrix of positions % mass is an N x 1 vector of masses % G is Newton's Gravitational constant % softening is the softening length % a is N x 3 matrix of accelerations % positions r = [x,y,z] for all particles x = pos(:,1); y = pos(:,2); z = pos(:,3); % matrix that stores all pairwise particle separations: r_j - r_i dx = x' - x; dy = y' - y; dz = z' - z; % matrix that stores 1/r^3 for all particle pairwise particle separations inv_r3 = (dx.^2 + dy.^2 + dz.^2 + softening.^2).^(-3/2); ax = G * (dx .* inv_r3) * mass; ay = G * (dy .* inv_r3) * mass; az = G * (dz .* inv_r3) * mass; % pack together the acceleration components a = [ax ay az]; end

时间整合(Time Integration)

The positions and velocities are updated using a leap-frog scheme (‘kick-drift-kick’). For each timestep Δt, each particle receives a half-step ‘kick’:

位置和速度使用跳越方案('kick-drift-kick')更新。 对于每个时间步长Δt，每个粒子都接收一个半步“踢”：

followed by a full-step ‘drift’:

然后是整步“漂移”：

followed by another half-step ‘kick’.

接下来是半步“踢”。

The evolution is performed in the code using a For-loop and our function for the acceleration from earlier:

使用For循环在代码中执行演化，而我们的函数则用于从早期开始进行加速：

%% Simulation Main Loop for i = 1:Nt % (1/2) kick vel = vel + acc * dt/2; % drift pos = pos + vel * dt; % update accelerations acc = getAcc( pos, mass, G, softening ); % (1/2) kick vel = vel + acc * dt/2; % update time t = t + dt; end

Simple, isn’t it?

很简单，不是吗？

The ‘kick-drift-kick’ is a highly accurate second-order technique that does a good job at preserving the total energy of the system.

“踢踢踢”是一种高度精确的二阶技术，可以很好地保持系统的总能量。

初始条件 (Initial Conditions)

The only thing left to carry out a simulation is to specify the initial positions and velocities of the particles at time t=0. Check out the code to play around with different setups. I initialize a random Gaussian distribution of positions and velocities for N=100 particles, but feel free to construct your own.

剩下要做的唯一模拟是指定在时间t = 0时粒子的初始位置和速度。 查看代码以使用不同的设置。 我为N = 100个粒子初始化了位置和速度的随机高斯分布，但可以自行构建。

能源诊断 (Energy Diagnostics)

Physically, the total energy of the system:

从物理上讲，系统的总能量为：

is conserved under time evolution. The first sum in this equation is the kinetic energy (KE) and the second is the potential energy (PE). Our code computes these quantities and keeps track of the total energy to make sure it is being approximately conserved by the numerical method. It is generally good practice to do this as a way to validate your code.

在时间演化下是保守的。 该方程式中的第一个和是动能(KE)，第二个是势能(PE)。 我们的代码将计算这些数量并跟踪总能量，以确保通过数值方法可以大致节省总能量。 通常，这样做是验证代码的一种好习惯。

Running the code allows you to visualize the simulation in real time and will yield the figure:

运行代码可以使您实时可视化仿真，并产生图：

Computational astrophysicists use this type of simulation every day as part of very detailed simulations of the Universe. Newton’s law of gravity is responsible for setting a lot of the structure you see in the Universe:

计算天体物理学家每天都使用这种类型的模拟作为非常详细的宇宙模拟的一部分。 牛顿的万有引力定律负责设定您在宇宙中看到的许多结构：

TNG Simulations TNG模拟中模拟宇宙中恒星的分布

演示地址

Download the Matlab code on github for our N-body simulation to visualize the simulation in real time and play around with the initial conditions. Enjoy!

在github上下载Matlab代码以进行N体仿真，以实时可视化仿真并尝试初始条件。 请享用！

翻译自: https://medium.com/swlh/create-your-own-n-body-simulation-with-matlab-22344954228e

matlab模拟

相关资源：matlab实现的人体识别与跟踪算法