在力学的学习过程中经典分析力学是最基本的入门知识，过冷水之前和大家一起学习了两个小车通过弹簧链接起来的做来回摆动运动的运动轨迹学习。推文中直接给了一个微分方程组，然后解出微分方程组就得到了小车的演化轨迹。本期过冷水从零开始构建一个微分方程组，而不是单纯解微分方程。

假设有一个一维的谐振子在三维空间中做运动。m=5 kg；k_x=2 N/m;k_y=7 N/m;k_z=4 N/m。

根据三维空间受力分析可得：

谐振子的势能函数为：

过冷水在这里需要给大家补充一点点初中弹簧做功的知识。W=Fr是做功的函数，已知F=kr;则W=kr²，而这里势能却是W=0.5kr²，这是因为错误使用W=Fr，该公式是恒力做功，而弹簧力是变力做功，需要使用积分来算功。

谐振子的动能为：

谐振子的总能量为：

谐振子的Lagrange函数为：

已知Lagrange方程方程形式：

Lagrange方程的证明较为复杂，在此过冷水只告诉打算使用了该方程而不证明方程的由来，感兴趣的可以找一本力学方面的书进行探究。Lagrange是一个微分方程，针对于我们的情况，可以很容易的就解出原函数。Lagrange方程的通解为：

设初始条件t=0时:

解得：

这样我们就求得了一个谐振子在三维空间的运动方程。可以得到如下的运动轨迹图

原始代码为：

clear
t=linspace(0,50,50000);
x=2.5*cos(sqrt(2/5)*t)+2*sqrt(5/2)*sin(sqrt(2/5)*t);
y=5*cos(sqrt(7/5)*t)+2*sqrt(5/7)*sin(sqrt(7/5)*t);
z=10*cos(sqrt(4/5)*t)+2*sqrt(5/4)*sin(sqrt(4/5)*t);
Coordiate=[x',y',z'];
[m,~]=find(Coordiate(:,1)==x(1)&Coordiate(:,2)==y(1)&Coordiate(:,3)==z(1))figure1 = figure('Name','谐振子运动','NumberTitle','off','Color','w','Menubar','figure','WindowState','maximized');
subplot1 = subplot(2,2,1,'Parent',figure1);
view(subplot1,[-37.5 30]);
subplot2 = subplot(2,2,2,'Parent',figure1);
subplot3 = subplot(2,2,3,'Parent',figure1);
subplot4 = subplot(2,2,4,'Parent',figure1);
hold(subplot1,'on');
hold(subplot2,'on');
hold(subplot3,'on');
hold(subplot4,'on');
gg=1;
for k=1:length(t)%fill3([0 14 14 0],[0 0 14 14],[0 0 0 0],[0 1 1],'Parent',axes1);%hold(axes1,'on')plot3%fill3([0 14 14 0],[14 14 14 14],[0 0 14 14],[0 0 1],'Parent',axes1);%fill3([14 14 14 14],[0 0 14 14],[0 14 14 0],[1 0 0],'Parent',axes1);plot3(x(k),y(k),z(k),'Parent',subplot1,'MarkerFaceColor',[0 1 1],'MarkerEdgeColor',[0 1 1],'MarkerSize',3,'Marker','o','LineStyle','none');xlim(subplot1,[min(x)-1 max(x)+1]);ylim(subplot1,[min(y)-1 max(y)+1]);zlim(subplot1,[min(z)-1 max(z)+1]);box(subplot1,'on');ylabel(subplot1,'$Y$','FontSize',20,'Interpreter','latex');xlabel(subplot1,'$X$','FontSize',20,'Interpreter','latex');zlabel(subplot1,'$Z$','FontSize',20,'Interpreter','latex');title(subplot1,'三维空间轨迹图','FontSize',20,'FontName','楷体');plot(t(k),z(k),'Parent',subplot2,'MarkerFaceColor',[0 1 1],'MarkerEdgeColor',[0 1 1],'MarkerSize',10,'Marker','o','LineStyle','none');ylim(subplot2,[min(z)-1 max(z)+1]);xlim(subplot2,[min(t)-2 max(t)+2])box(subplot2,'on');xlabel(subplot2,'$T$','FontSize',20,'Interpreter','latex');zlabel(subplot2,'$Z$','FontSize',20,'Interpreter','latex');title(subplot2,'TZ图','FontSize',20,'FontName','楷体');plot(t(k),x(k),'Parent',subplot3,'MarkerFaceColor',[0 1 1],'MarkerEdgeColor',[0 1 1],'MarkerSize',10,'Marker','o','LineStyle','none');ylim(subplot3,[min(x)-1 max(x)+1]);xlim(subplot3,[min(t)-2 max(t)+2])box(subplot3,'on');xlabel(subplot3,'$T$','FontSize',20,'Interpreter','latex');zlabel(subplot3,'$X$','FontSize',20,'Interpreter','latex');title(subplot3,'TX图','FontSize',20,'FontName','楷体');plot(t(k),y(k),'Parent',subplot4,'MarkerFaceColor',[0 1 1],'MarkerEdgeColor',[0 1 1],'MarkerSize',10,'Marker','o','LineStyle','none');ylim(subplot4,[min(y)-1 max(y)+1]);xlim(subplot4,[min(t)-2 max(t)+2])box(subplot4,'on');xlabel(subplot4,'$T$','FontSize',20,'Interpreter','latex');zlabel(subplot4,'$Y$','FontSize',20,'Interpreter','latex');title(subplot4,'TY图','FontSize',20,'FontName','楷体');fmat(:,k)=getframe;           %拷贝祯到矩阵fmat中im = frame2im(getframe);[I,map] = rgb2ind(im,256);if gg == 1imwrite(I,map,'lxx.gif','GIF', 'Loopcount',inf,'DelayTime',0.05);gg = gg + 1;elseimwrite(I,map,'lxx.gif','WriteMode','append','DelayTime',0.05);end
end
movie(fmat,1)

根据轨迹图像可以看出谐振子在三维空间中运动轨迹很复杂。这只是动力学上比较简单一个运动问题。过冷水要带大家做的是有很多个相互作用的空间运动问题。不积跬步无以至千里，本期先和大家分享这个简单的案例，下期带大家了解两个谐振子的一维运动图像。想要深入了解学习请持续关注过冷水的推文。

往期回顾>>>>>>

数值优化方法—迭代法&终止条件

数值优化方法及MATLAB实现（一）

十万个matlab编程问题征集，欢迎来问

数值优化—复杂函数重积分计算方法实例演示

MATLAB编程爱好者Q群正式升级为2000人群，欢迎加入！