大二数学物理方法之单摆---matlab研究二阶常微分方程

本人大二。这是我的数学物理方法课程作业，参考教材是顾樵的书。代码在文末 参考视频【官方双语】微分方程概论-第一章】 https://www.bilibili.com/video/BV1tb411G72z/?share_source=copy_web&vd_source=4d4acc16b479cfd36366b0862b8e8dea

function  dxdt= odefun( t,x,a,b )
%化二阶常微分方程为一阶
%   
dxdt=zeros(2,1);
dxdt(1)=x(2);
dxdt(2)=-a*sin(x(1))-b*x(2);

end

function  dxdt= odefun2( t,x,a,b )
%化二阶常微分方程为一阶sinx代替x
%   
dxdt=zeros(2,1);
dxdt(1)=x(2);%相当于dx(1)/dt,等号右边是x1x2
dxdt(2)=-a*x(1)-b*x(2);%相当于dx(2)/dt

end

 开头这两个封装好的函数需要单独新建两个脚本分别存储，在主程序才能正常调用 

a=3;%系统本征值g/l的平方项，g为重力加速度，l为单摆摆长
 b=2; %阻尼系数b>=0，b越大，阻力项f=b*dxdt越大
 tspan=[0 30];%单摆运动时间段
 x0=[pi/12;0];%一次实验的两个初始参数[初角位移x(1)，初角速度x(2)]，x(2)=dx(1)/dt
[t1,x1]=ode45(@(t,x)odefun(t,x,a,b),tspan,x0);%求解关于单摆精确的的角位移与时间的常微分方程，该odefun，能输入两个初始变量，这些变量都是t的函数
[t2,x2]=ode45(@(t,x)odefun2(t,x,a,b),tspan,x0);%求解小角度线性系统拟合，令x代替sinx
%hold on
subplot(2,1,1);
plot(t1,x1(:,1),'-',t2,x2(:,1),'-');%输出的t*2数组中，第一列为对应初始参数x(1)的函数值，第二列为对应初始参数x(2)的函数值
title(['Solution of A Pendulum with Initial AD=',num2str(x0(1)),'rad'])%title，label等修饰要跟在plot后面，使用字符串元胞数组
xlabel('time/s');
ylabel('Angular Displacement/rad'); 
legend('nonlinear','linear');
%画出两次实验角位移（rad)对时间(s)的函数
subplot(2,1,2)
plot(t1,x1(:,2),'-',t2,x2(:,2),'-')
title(['Solution of A Pendulum with Initial AD dot=',num2str(x0(2)),'rad/s'])%title，label等修饰要跟在plot后面
%hold off;
xlabel('time/s');
ylabel('Angular Speed/(rad/s)'); 
legend('nonlinear','linear');

%相图绘制
%a=3;%系统本征值g/l的平方项，g为重力加速度，l为单摆摆长
%b=1; %阻尼系数b>=0，b越大，阻力项f=b*dxdt越大
f=@(t,x)odefun(t,x,a,b);%f为关键函数传导中介
%Y(2)为阻尼项，该定义方式中，每项系数固定为单位一
y1 = linspace(-2,8,20);%生成-2到8的20个节点，组成向量;y1可作为横轴x轴轴指代单摆距离平衡位置的角位移
y2 = linspace(-4,4,20);%y2作为纵轴y轴指代角速度，详情可见推荐3b1b视频
% creates a matrix  for all the x-values and for
% all the y-values on the grid. Note that x and y are matrices of the same
% size and shape, in this case 20 rows and 20 columns
[x,y] = meshgrid(y1,y2);%y1,y2为向量,生成20*20的x、y坐标空间
u = zeros(size(x));%size(x)=[20 20]设置初始向量的x分量为零
v = zeros(size(x));

% we can use a single loop over each element to compute the derivatives at each point (y1, y2)
t=0; % we want the derivatives at each point at t=0, i.e. the starting time（力场与时刻无关）
for i = 1:numel(x)%numel返回元素总数（400个零）/length返回size中较大数
    Ydot = f(t,[x(i); y(i)]);%400对坐标根据f求得的x(1)、x(2)获得初始位矢及角速度
    u(i) = Ydot(1);%分发x、y轴分量
    v(i) = Ydot(2);
end
quiver(x,y,u,v,'r'); figure(gcf)%,'r'调整颜色
grid
xlabel('角位移y_1')
ylabel('角速度y_2')
title(['单摆相图 阻尼系数=',num2str(b)])
axis tight equal;%Image Compacted