鵲橋仙 -- 秦觀(guān)

纖云弄巧，飛星傳恨，銀漢迢迢暗度。金風(fēng)玉露一相逢，便勝卻人間無(wú)數(shù)。

柔情似水，佳期如夢(mèng)，忍顧鵲橋歸路。兩情若是久長(zhǎng)時(shí)，又豈在朝朝暮暮。

《古詩(shī)十九首》云：“河漢清且淺，相去復(fù)幾許？盈盈一水間，脈脈不得語(yǔ)。

牛郎織女之所以不能朝夕相處，琴瑟和鳴，最主要的原因是大河相隔，因此要解決這個(gè)問(wèn)題的一個(gè)重要方法時(shí)修橋。 技術(shù)發(fā)展日新月異，實(shí)際上修建一個(gè)現(xiàn)代化的大橋?qū)ε＠煽椗纳钣泻艽蟮母纳?。在修橋的時(shí)候，模擬仿真軟件發(fā)揮著重要的作用。 做打油詩(shī)一首，

下面，我簡(jiǎn)單模擬一下，橋由于隨機(jī)的初始速度，然后橋隨著時(shí)間的抖動(dòng)情況。100米的大橋被簡(jiǎn)化成7個(gè)節(jié)點(diǎn)，用牛頓第二定律和胡克定律來(lái)建模， 用verlet 方法數(shù)值求解常微分方程組。

使用hangingBridge.m將初始位置初始化為平衡位置，并使用寬度dv=0.5m/s的高斯隨機(jī)分布初始化初始速度(gauss分布的寬度就是根方差)。還定義一個(gè)函數(shù)force =@（r）bridgeForces（r，m，k），并調(diào)用所提供的verlet函數(shù)，以獲得tmax=10s期間節(jié)點(diǎn)的位置，間隔dt=0.1s? 畫(huà)這座橋各個(gè)節(jié)點(diǎn)的位置隨著時(shí)間的變化。?

北太天元的代碼如下：

% 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('初始時(shí)刻的節(jié)點(diǎn)位置', ['時(shí)間=', num2str((nt-1)*dt), ?'時(shí)的節(jié)點(diǎn)位置']);

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