FPUT问题：重复FPUT回归（MATLAB代码）

fengjr_lzu

已于 2022-07-25 17:13:13 修改

阅读量1k

点赞数

文章标签： matlab

于 2022-07-25 17:00:05 首次发布

本文链接：https://blog.csdn.net/yqx_1128/article/details/125976294

版权

FPUT是四个人名的缩写：Fermi-Pasta-Ulam-Tsingou，该问题是指某些非线性系统不会分散它们的能量（能量均分状态），而是回到初始激发状态，这一违背常理的现象被称为FPUT问题。FPUT问题直接导致了非线性科学的诞生，并引领人们进入了使用计算机进行科学计算的时代。

1955年费米(Fermi)，帕斯塔(Pasta)和乌拉姆(Ulam)为了验证经典统计物理中能量均分原理而设计了一个数值试验，他们将64个质点用非线性弹簧连成一个振动弦, 开始能量只集中在第一模态, 按照能量均分原理, 由于弱非线性相互作用, 长时间以后应该导致能量有涨落的平均分布到所有模态从而达到能量均分状态。开始的计算结果确实表明能量在向其他模态转移, 但是出于意料的事, 长时间演化以后几乎全部能量又回到了初始分布状态, 这和能量均分原理完全不一致。这一“ 违背常理” 的现象后来以三位发现者姓名的首字母命名为FPU重现(FPU recurrent)，这一问题被称为Fermi-Pasta-Ulam-Tsingou问题。此外, 由于该数值试验是由Mary Tsingou编程实现并将结果最终以图表的形式展示出来, 考虑到在第一台计算机上设计第一个数值实验在当时的难度, Dauxois呼吁应该将其称为Fermi-Pasta-Ulam-Tsingou问题，以此表示对Tsingou在该研究中所做贡献的尊重。

将若干个质点用非线性弹簧连成一个振动弦

图片来源：http://www.scholarpedia.org/

系统哈密顿量

$H=\sum_{i=1}^N \frac{x_i^2}{2}+\sum_{i=0}^N[\frac{1}{2} (x_{i+1}-x_i )^2+\frac{\alpha}{3} (x_{i+1}-x_i )^3+\frac{\beta}{4} (x_{i+1}-x_i )^4]$

$\omega_i=2sin\frac{i\pi}{2(N+1)}$

$\psi_{ij}=\sqrt{\frac{2}{N+1}}sin\frac{ij\pi}{N+1}$

可移动粒子数N=31，即共有33个粒子

初始位移：第一个本征模式

$x_i=sin\frac{i\pi}{32}$

初始速度：所有原子均为0

模拟步长：0.01

模拟步数：1000000

计算时要注意的几个要点

模拟结果

源代码（MATLAB）

clear;
clc;
N=31;
dt=1;ddt=0.01;%步长0.01
totalt=1000000;%步数1000000
alpha=0.25;
beta=0;%只计算到2阶
omega=zeros(N+2,1);
psi=zeros(N+2,N+2);
xx=zeros(N+2,totalt);
vv=zeros(N+2,totalt);
aa=zeros(N+2,totalt);
%第一个时刻：初始化
for ii=2:N+1
    omega(ii)=2*sin((ii-1)*pi/(2*(N+1)));
    xx(ii,1)=sin((ii-1)*pi/32);
    vv(ii,1)=0;
end
for ii=2:N+1
    for jj=2:N+1
        psi(ii,jj)=sqrt(2/(N+1))*sin((ii-1)*(jj-1)*pi/(N+1));
    end
end
%用牛顿运动定律，粗略算出第二个时刻的位置、速度
for ii=2:N+1
    aa(ii,1)=(xx(ii+1,1)+xx(ii-1,1)-2*xx(ii,1))+alpha*((xx(ii+1,1)-xx(ii,1))^2-(xx(ii,1)-xx(ii-1,1))^2)+beta*((xx(ii+1,1)-xx(ii,1))^3-(xx(ii,1)-xx(ii-1,1))^3);
    xx(ii,2*dt)=xx(ii,1)+vv(ii,1)*ddt+0.5*aa(ii,1)*ddt^2;
    vv(ii,2*dt)=vv(ii,1)+aa(ii,1)*ddt;
end
%verlet算法，用前两个时刻的位置算出第三个时刻的位置，以此类推，算出所有时刻的位置、速度、加速度
for t=2:dt:totalt
    for ii=2:N+1
        aa(ii,t)=(xx(ii+1,t)+xx(ii-1,t)-2*xx(ii,t))+alpha*((xx(ii+1,t)-xx(ii,t))^2-(xx(ii,t)-xx(ii-1,t))^2)+beta*((xx(ii+1,t)-xx(ii,t))^3-(xx(ii,t)-xx(ii-1,t))^3);
        xx(ii,t+dt)=2*xx(ii,t)-xx(ii,t-dt)+aa(ii,t)*ddt^2;
        vv(ii,t)=(xx(ii,t+dt)-xx(ii,t-dt))/(2*ddt);
    end
end

ii=1:N+2;
%h=plot(ii,xx(:,34));
qq=zeros(1,N+2);pp=zeros(1,N+2);energy=zeros(N+2,totalt);
sum_energy=zeros(totalt,1);

for t=1:dt:totalt
    qq=psi*xx(:,t);%将普通坐标变换为简正坐标
    pp=psi*vv(:,t);
    energy(:,t)=0.5*(pp.*pp+omega.*omega.*qq.*qq);
end
%取前五个简正模式，绘图
plot(1:dt:totalt,energy(1:5,:))