Advertisement
Y-BH

N_pendulum_symbol

May 8th, 2022
1,034
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. %% N摆的拉格朗日方程
  2. %% ——用MATLAB求解
  3.  
  4. clc;clear
  5. N=3;
  6. syms a [1 N]
  7. syms b [1 N]
  8. syms c [1 N]
  9. syms theta(t) [1 N]
  10. syms t
  11. % 几何约束
  12. for k=1:N
  13.     if k==1
  14.         x(k)=sin(eval(['theta',num2str(k)]));
  15.         y(k)=-cos(eval(['theta',num2str(k)]));
  16.  
  17.     else
  18.     x(k)=x(k-1)+sin(eval(['theta',num2str(k)]));
  19.     y(k)=y(k-1)-cos(eval(['theta',num2str(k)]));
  20.     end
  21.     E(k)=Ek(x(k),y(k),t);
  22. end
  23. % 势能、动能、拉格朗日量
  24. V=sum(y);T=sum(E);L=T-V;
  25. % 拉格朗日方程
  26. for k=1:N
  27.     D2(k)=diff(eval(['theta',num2str(k)]),t,2);
  28.     D1(k)=diff(eval(['theta',num2str(k)]),t);
  29.     eq(k)=diff(diff(L,D1(k)),t)-diff(L,eval(['theta',num2str(k)]));
  30. end
  31. % 获得动能函数,用以验证能量是否守恒
  32. Tf=subs(T,[D1,theta],[c,b]);% 符号替换
  33. Tf=simplify(Tf);
  34. Tf=matlabFunction(Tf);
  35. for k=1:2*N
  36.     if k==1
  37.         x=['u(',num2str(k),')'];
  38.     else
  39.         x=[x,',u(',num2str(k),')'];
  40.     end
  41. end
  42. Tkf=eval(['@(u) Tf(',x,')']);
  43. % 利用拉格朗日方程解出二阶导数表达式
  44. eq=subs(eq,[D1,theta,D2],[c,b,a]);% 符号替换
  45. eq=simplify(eq);% 简化方程
  46. tic;D=solve(eq,a);toc% 求解
  47. % 生成一阶导数的系数矩阵
  48. Dn=[];
  49. for k=1:N
  50.     Dk=eval(['D.a',num2str(k)]);
  51.     Dk=coeffs(Dk,c);
  52.     Dn=[Dn;Dk];
  53. end
  54. tic;Dn=simplify(Dn);toc
  55. % 转为函数脚本,以供数值计算使用
  56. ddn=matlabFunction(Dn);
  57. % 转化为向量为自变量的函数
  58. for k=1:N
  59.     if k==1
  60.         x=['u(',num2str(k),')'];
  61.     else
  62.         x=[x,',u(',num2str(k),')'];
  63.     end
  64. end
  65. f=eval(['@(u) ddn(',x,')']);
  66. save func_N3k.mat f N Tkf
  67. %%
  68. function E=Ek(x,y,t)
  69. % 动能函数
  70. E=1/2*((diff(x,t))^2+(diff(y,t))^2);
  71. end
Advertisement
RAW Paste Data Copied
Advertisement