中国情人节-七夕-鹊桥--北太天元为牛郎织女修桥助力
鹊桥仙 -- 秦观
纤云弄巧，飞星传恨，银汉迢迢暗度。金风玉露一相逢，便胜却人间无数。
柔情似水，佳期如梦，忍顾鹊桥归路。两情若是久长时，又岂在朝朝暮暮。
《古诗十九首》云：“河汉清且浅，相去复几许？盈盈一水间，脉脉不得语。
牛郎织女之所以不能朝夕相处，琴瑟和鸣，最主要的原因是大河相隔，因此要解决这个问题的一个重要方法时修桥。 技术发展日新月异，实际上修建一个现代化的大桥对牛郎织女的生活有很大的改善。在修桥的时候，模拟仿真软件发挥着重要的作用。 做打油诗一首，
下面，我简单模拟一下，桥由于随机的初始速度，然后桥随着时间的抖动情况。100米的大桥被简化成7个节点，用牛顿第二定律和胡克定律来建模， 用verlet 方法数值求解常微分方程组。
使用hangingBridge.m将初始位置初始化为平衡位置，并使用宽度dv=0.5m/s的高斯随机分布初始化初始速度(gauss分布的宽度就是根方差)。还定义一个函数force =@（r）bridgeForces（r，m，k），并调用所提供的verlet函数，以获得tmax=10s期间节点的位置，间隔dt=0.1s 画这座桥各个节点的位置随着时间的变化。
北太天元的代码如下：
% hangingBridge.m 适用于北太天元
% Solution to problem 2 MT2 Phy277 Spring 2017.
% Model of a hanging bridge.
% MVFS Aprl 2017
% Set parameters
clear
clf
close all
n = 7; % number of bridge nodes, including end ones
m = 1e2; % mass of each node (kg)
g = 9.8; % gravitational acceleration (m/s^2)
k = 5e2; % spring constant (N/m)
L = 100; % total bridge length (m)
dv = 0.5; % average initial velocities (m/s)
dt = 0.02; % time interval (s)
tmax = 10; % time span (s)
% Find the equilibrium bridge shape by solving the n equations
% x(1) = 0
% -k*(x(i)-x(i-1)) - k*(x(i)-x(i+1)) = 0, i=2:n-1
% x(n) = L
% whose solution is simply x(i)=(i-1)*L/(n-1). Equally
% y(1) = 0
% -k*(y(i)-y(i-1)) - k*(y(i)-y(i+1)) = m*g
% y(n) = 0
% for the n unknowns y(i). We write these eqs. as a*y=b
a = zeros(n,n);
b = zeros(n,1);
for i = 2:n-1
a(i,i) = -2*k;
a(i,i-1) = +k;
a(i,i+1) = +k;
b(i) = m*g;
end
a(1,1) = 1;
a(n,n) = 1;
y = a\b;
y = y.' ; % row vector
x = linspace(0,L,n); % uniformly separated x coordinates
% Check that forces are zero for debugging
% f = bridgeForces( [x;y], m, k )
% Plot cable and nodes
plot(x,y,'b-o','LineWidth',2);
%axis equal % to show real shape
%axis([0,L,-50,10]) % consistent with movie window
xlabel('x (m)')
ylabel('y (m)')
r = [x;y];
force = @(r) bridgeForces(r,m,k);
v0 = dv* randn(size(r));
r0 = r;
mass = repmat(m, 1, n);
[r,v,a] = verlet( force, mass, r0, v0, dt, tmax );
np = size(r0, 2);
N = size(r,2)/ np;
load_plugin("time");
for nt = 1:N
if( mod((nt-1)*dt, 0.1) ~= 0 )
continue;
end
ind = (nt-1)*np+1 : nt*np;
plot(r0(1, :), r0(2,:), 'k-.', 'LineWidth', 1);
hold on
plot(r(1, ind), r(2,ind), 'r-*', 'LineWidth', 3);
nt
(nt-1)*dt
num2str((nt-1)*dt)
legend('初始时刻的节点位置', ['时间=', num2str((nt-1)*dt), '时的节点位置']);
pause(0.5);
hold off
end
%% bridgeForces.m
function f = bridgeForces( r, m, k )
% Forces on a hanging bridge
% Input:
% r(2,n) : coordinates of bridge nodes (m)
% m : nodes’ mass (kg)
% k : cable's spring constant (N/m)
% Output:
% f(2,n) : force on each node, with f(:,1)=f(:,n)=0
g = [0;-9.8]; % gravitational acceleration vector (m/s^2)
n = size(r,2); % number of nodes
% Find forces
f = zeros(size(r));
for i = 2:n-1
f(:,i) = - k*(r(:,i)-r(:,i-1)) - k*(r(:,i)-r(:,i+1)) + m*g(:);
end
end
%% verlet.m
function [r,v,a] = verlet( force, mass, r0, v0, dt, tmax )
% Solves Newton's equation using velocity Verlet method for np particles
% in nd dimensions. MVFS Aprl 2017
%
% Input:
% force : function to calculate the force as f=force(r), in N
% mass(np) : particle masses of the np particles in kg
% r0(nd,np) : intial positions of np particles in nd dimensions, in m
% v0(nd,np) : initial velocities in m/s
% dt : integration time interval in s
% tmax : final integration time in s
% Output:
% r(nd,np*nt) : position vectors at each time (0:dt:tmax), in m
% v(nd,np*nt) : velocities at each time, in m/s
% a(nd,np*nt) : accelerations at each time, in m/s^2
% Check array sizes
[nd,np] = size(r0);
if (numel(mass)~=np)
'verlet: ERROR: numel(mass)~=np'
end
if (size(v0)~=size(r0))
'verlet: ERROR: size(v0)~=size(x0)'
end
% Set auxiliary variables and arrays
t = (0:dt:tmax); % integration times
nt = numel(t); % number of times
% Assign a mass to each coordinate
m = repmat(mass,nd,1);
% Integrate trajectory
r=zeros(nd,np*nt); v=r; a=r; % set size of output arrays
r(:,1:np) = r0; % positions at t=0
v(:,1:np) = v0; % velocities at t=0
a(:,1:np) = force(r0)./m; % accelerations at t=0
for it = 2:nt
r(:, np_it(np,it) ) = r(:,np_it(np, it-1) ) + ...
v(:,np_it(np, it-1) )*dt + ...
a(:,np_it(np, it-1) )*dt^2/2; % positions at t
a(:,np_it(np, it) ) = force(r(:,np_it(np, it) ))./m; % accelerations at t
v(:,np_it(np, it) ) = (r(:,np_it(np, it) )-r(:,np_it(np, it-1) ))/dt + ...
a(:,np_it(np, it) )*dt/2; % velocities at t
end
end % function verlet
function ind = np_it( np, it)
ind = ( (it-1)*np+1 ) : ( it*np );
end
