pendulum_period

function pendulum_period
% pendulum waves-寻找合适的摆长
clc;clear;close all
g=10;         % 重力加速度
Tend=200;     % 大于最长摆长的周期/4的一个时间
T=80;        % 蛇摆周期
theta=pi/6;   % 最大摆角
L=0.5:0.01:3; % 摆长范围
% ode设置
opts=odeset('RelTol',1e-13,'AbsTol',1e-14, ...
    'MaxStep',1e-3,'InitialStep',1e-4, ...
    'Stats','off','Events',@(t,u)evns(t,u));
for ii=1:length(L)
    [t,~]=ode45(@(t,u)pw(t,u,g,L(ii)), ...
        [0 Tend],[theta 0],opts);
    w(ii)=T/4/t(end); % w:各摆在T内摆动次数
end
w_max=floor(max(w)); % 最大摆动次数，整数
w_min=ceil(min(w));  % 最小摆动次数，整数
% 计算整数摆动次数对应摆长
L=interp1(w,L,w_min:w_max);
% 保存变量用以模拟调取
save matlab.mat L T theta g
%% 微分方程与事件函数函数
function du=pw(~,u,g,L)
du=[u(2);-g/L*sin(u(1))];
function [position,isterminal,direction] = evns(~,u)
position = u(1);isterminal = 1;direction = -1;