Matlab求解傅科摆动
过冷水诚挚邀请你加入Matlab仿真秀官方交流群进行Matlab学习、问题咨询、 Matlab相关资料下载，**：927550334
傅科摆是科学史上一个著名的实验，这个实验是人类第一次不借助地球以外的参照物证明地球自转的伟大尝试。借助理论力学知识容易列出傅科摆的动力学方程，但要从方程中得到其运动轨迹，则需要求解复杂的微分方程， 通常情况下不借助计算机而直接求解微分方程是十分复杂的，过冷水最近一直学习微分方面知识，微分在实际工业中应用中最普遍的就是解动力学方程，本期过冷就和大家分享用不同的方法解傅科摆动的运动轨迹。
设傅科摆摆长为l，摆锤质量为m，悬挂于北纬ϕ度处。由于摆长远大于摆锤尺寸，当摆锤做小角度摆动时，可认为摆锤在水平面内运动。以摆锤平衡位置为原点O，Ox指向正南方向，Oy指向正东方向，建立直角坐标系。
摆锤受重力、科氏惯性力和摆线张力，假设摆线张力F≈ ~mg，可得出动力学微分方程为:
g为重力加速度; ω 为地球自转角速度;φ为实验所在地的纬度。这样的方程组求解非常烦琐，用数学方法直接求出其通解有难度，通解的形式也很复杂。如果限定一定的初始条件，例如限定摆动初始条件为:x(0)=A,y(0)=0，x`(0)=B,y`(0)=C,知道该方程条件后，可以得到对应解：
这个是较为准确的解析解，实际情况获得解析解过程很复杂，使用Matlab进行数值解计算。Matlab中dsolve函数可以用于求解符号微分方程组，并能给出其对应解析解：
y = dsolve('eq1,eq2,...' , 'cond1,cond2,...', 'x')
dsolve:该函数用来求解符号常微分方程或方程组;
'eql,eq2,..': 待求解微分方程（组);
'cond1,cond2,..….':初始（边界）条件;
'x':是独立变量，空缺时默认为't';
'y :求得的微分方程通解。`
[xcommon ycommon] = dsolve('D2x-2*w*Dy*sin(a) g/l*x = 0','D2y 2*w*Dx*sin(a) g/l*y = 0');dsolve函数还存在另外一种写法，直接将表达式定义为符号表达式，将上述方程组进行变形：
将该方程写成符号运算形式，然后带入dsolve函数中
syms y(t) g l a w x(t) f1=diff(y,2,t) == -2*w*diff(x,1,t)*sin(a)-g/l*y; f2=diff(x,2,t) == -2*w*diff(y,1,t)*sin(a)-g/l*x; [xcommon ycommon] = dsolve(f1,f2) xcommon = C1*exp(-(t*((-l*(g - l*w^2*sin(a)^2))^(1/2) l*w*sin(a)))/l)*(((-l*(g - l*w^2*sin(a)^2))^(1/2) l*w*sin(a))/g - (2*l*w*sin(a))/g) - C3*exp((t*((-l*(g - l*w^2*sin(a)^2))^(1/2) l*w*sin(a)))/l)*(((-l*(g - l*w^2*sin(a)^2))^(1/2) l*w*sin(a))/g - (2*l*w*sin(a))/g) - C2*exp((t*((-l*(g - l*w^2*sin(a)^2))^(1/2) - l*w*sin(a)))/l)*(((-l*(g - l*w^2*sin(a)^2))^(1/2) - l*w*sin(a))/g (2*l*w*sin(a))/g) C4*exp(-(t*((-l*(g - l*w^2*sin(a)^2))^(1/2) - l*w*sin(a)))/l)*(((-l*(g - l*w^2*sin(a)^2))^(1/2) - l*w*sin(a))/g (2*l*w*sin(a))/g) ycommon = C1*exp(-(t*((-l*(g - l*w^2*sin(a)^2))^(1/2) l*w*sin(a)))/l)*(((-l*(g - l*w^2*sin(a)^2))^(1/2) l*w*sin(a))/g - (2*l*w*sin(a))/g) C3*exp((t*((-l*(g - l*w^2*sin(a)^2))^(1/2) l*w*sin(a)))/l)*(((-l*(g - l*w^2*sin(a)^2))^(1/2) l*w*sin(a))/g - (2*l*w*sin(a))/g) - C2*exp((t*((-l*(g - l*w^2*sin(a)^2))^(1/2) - l*w*sin(a)))/l)*(((-l*(g - l*w^2*sin(a)^2))^(1/2) - l*w*sin(a))/g (2*l*w*sin(a))/g) - C4*exp(-(t*((-l*(g - l*w^2*sin(a)^2))^(1/2) - l*w*sin(a)))/l)*(((-l*(g - l*w^2*sin(a)^2))^(1/2) - l*w*sin(a))/g (2*l*w*sin(a))/g)
一般情况下，很难从上述解析解表达 直观看出摆锤运动形式。可借助Matlab可视化绘制出对应运动形式。 可视化最重要的两个步骤是画摆绳和摆锤，并让摆绳和摆锤按上述表达式在不同的时刻出现在不同的位置。画摆锤和摆绳可以用MATLAB的 line函数实现。
line(x,y,z,'PropertyName','propertyvalue’)
在三维坐标中绘制线条，如果 x 和 y ，z是矩阵则 line 将绘制多个线条。与 plot 函数不同，line 会向当前坐标区添加线条，而不删除其他图形对象或重置坐标区属性。
'PropertyName':需要设置的属性的名称。propertyvalue:设置属性具体值。
对应程序:
figure1 = figure;
colormap(gray);
axes1 = axes('Parent',figure1,'units','normalized','position',[0 0 1 1]);
uistack(axes1,'down')
image(blackboard,'Parent',axes1,'CDataMapping','scaled')
set(axes1,'handlevisibility','off','visible','off');
axes2 = axes('Parent',figure1);
pendulumline = line([0,10],[0,10],[0,10],'Parent',axes2,'color','y','linestyle', '-','linewidth',6); % 摆线
pendulumbob = line(10,10,10,'Parent',axes2,'color','r','Marker','.', 'Markersize',100); % 摆锤初始化
view( axes2,[15 30])
box(axes2,'on')
set(axes2,'Color','none','LineWidth',6,'TickLength',[0.01 0.005],'XColor',[1 0 1],'YColor',[1 0 1],'ZColor',[1 0 1]);在新的时刻，通过set函数设置摆线和摆锤的line对象的xdata、ydata、zdata属性，可以画一条新的摆线和摆锤,这就实现了动态显示。为更明显地显示摆平面变化，可通过line函数将摆锤的轨迹记录下来，进行实时更新。
以真实傅科摆(摆长67 m,摆锤质量28 kg)为例，对上述解析解进行可视化显示，真实傅科摆的摆平面转动周期非常长（为32 h)，为了使摆平面转动效果更明显，可在程序中增大地球自转角速度，使程序显示一种“夸大”的摆平面转动过程。我们模拟纬度值45°，x0=2，y0=2,u0=0, v0=0，且地球自转角速度为2π/100摆锤的运动轨迹。
其对应程序是：
a = 45; %纬度;
x0=2; % x 方向初始坐标;
y0=2;%y 方向初始坐标
u0=0 ;%x 方向初始速度
v0=0; %y 方向初始速度
c = a*pi/180;
w = pi/50;
l = 67;
t = 0:0.1:180;
g = 9.8;
[x y] = dsolve('D2x-2*w*Dy*sin(a) g/l*x = 0', 'D2y 2*w*Dx*sin(a) g/l*y = 0',['x(0)=',num2str(x0)],['y(0)=',num2str(y0)],['Dx(0)=',num2str(u0)],['Dy(0)=',num2str(v0)]);
funx=inline(x);funy=inline(y);
dx=real(funx(a,g,l,t,w));dy=real(funy(a,g,l,t,w));
dz =-(l^2-dx.^2-dy.^2).^0.5;
blackboard=imread('blackboard.jpg');
figure1 = figure;
colormap(gray);
axes1 = axes('Parent',figure1,'units','normalized','position',[0 0 1 1]);
uistack(axes1,'down')
image(blackboard,'Parent',axes1,'CDataMapping','scaled')
set(axes1,'handlevisibility','off','visible','off');
axes2 = axes('Parent',figure1);
view( axes2,[15 30])
box(axes2,'on')
set(axes2,'Color','none','LineWidth',6,'TickLength',[0.01 0.005],'XColor',[1 0 1],'YColor',[1 0 1],'ZColor',[1 0 1]);
verticalaxis = line([0,0],[0,0],[0,min(dz)],'color','k','linestyle','-','linewidth',1.5); % 竖轴
pendulumline = line([0,dx(1)],[0,dy(1)],[0,dz(1)],'color','y','linestyle', '-','linewidth',2); % 摆线初始化
pendulumbob = line(dx(1),dy(1),dz(1),'color','r','marker','.', 'markersize',25); % 摆锤初始化
pendulumtra = line([dx(1),dx(1)],[dy(1),dy(1)],[dz(1),dz(1)]); % 摆锤起点
axis ([min(dx) max(dx) min(dy) max(dy) min(dz) 0])
hold on
for i = 1:size(t,2)
set(pendulumline,'xdata',[0,dx(i)],'ydata',[0,dy(i)],'zdata',[0,dz(i)]); % 更新摆线位置
set(pendulumbob,'xdata',dx(i),'ydata',dy(i),'zdata',dz(i));% 更新摆锤位置
set(pendulumtra,'xdata',dx(1:i),'ydata',dy(1:i),'zdata',dz(1:i))
drawnow;
End微分数值方法
微分数值方法,不同于解析方法,其适合求解没有或很难求出微分表达式的实际问题计算。解析法中函数的导数是通过取极限来计算。当函数形式不明确,如以表格给出自变量和因变量关系时就不能通过解析法求导数,实际工作中通常是需要求列表函数在节点和非节点处的导数值,这正是数值微分所要解决的问题。数值微分方法可近似求出某点的导数值,或者将函数在某点的导数用该点附近节点上的函数近似表示。比如动力学中研究物体速度。在数值解法中,首先把区间[t0，tf]插入一系列分点ti,使t0<t1<…<ti<…<tf; 记hi;=ti 1-ti;，hi称为步长。所谓数值解就是求取方程满足定解条件下函数在节点ti,上的近似值ui。
欧拉法
将微分问题进行两端积分，
以上即为由y(ti)求y(ti 1)关系式,but ! f(t,y)是未知的,数值积分法通常是采用矩形法，把区间[ti, ti 1]上 f(t,y)近似看成是常数,这样
上述公式即给出了由y(ti)求y(ti 1)近似值方法,称为欧拉法。欧拉法对f(t，x)积分实际使用矩形公式,被积函数在求积区间的左端点被计算一次。当f(t，x)为常数该方法准确,当f(t，x)为线性时则不然,并且误差和h=(ti 1-ti)成正比,为获得一定位数的精度需要选很小步长。欧拉法最大缺点是没有提供误差估计方法,也就无法自动确定步长以得到期望的精度。
Runge-Kutta法
如果在欧拉法基础上再增加一次函数求值,则可能得到一个解决办法。类似于积分中的中点公式和梯形公式,这里也有两种选择。类比中点公式,先用欧拉法计算区间的一半，在中点处估算函数值v然后再用这里的斜率,进行实际的一步计算:
德国数学家C，Runge及M.W.Kutta提出对ti和ti 1之间几个不同t值估算f(t,x),然后由这些f(t,x)值的线性组合加上yi,得到需要的yi 1值。这套思想逐渐成为求解常微分方程问题所采用的主要算法，其对应的数学形式为：
根据上述思想可编出数值法求解微分方程程序
function dy=odefun1(x)
rho=1;
omega=pi/50;phi=45*pi/180;g=9.8;l=67;
dy(1)=x(2);
dy(2)=2*omega*sin(phi)*x(4)-g/l*x(1);
dy(3)=x(4);
dy(4)=-2*omega*sin(phi)*x(2)-g/l*x(3);
dy=[dy(1);dy(2);dy(3);dy(4)];
end
clear y
warning off all
f=@odefun1
h=0.1;l=67;
t=0:h:60;
x0=[2 0 2 0]
y(:,1)=x0;
for i=1:length(t)-1
L1=f(y(:,i));
L2=f(y(:,i) (h/2)*L1);
L3=f(y(:,i) (h/2)*L2);
L4=f(y(:,i) h*L3);
y(:,i 1)=y(:,i) (h/6)*(L1 2*L2 2*L3 L4);
end
z =-(l^2-y(1,:).^2-y(3,:).^2).^0.5;
figure1 = figure;
colormap(gray);
axes1 = axes('Parent',figure1,'units','normalized','position',[0 0 1 1]);
uistack(axes1,'down')
image(blackboard,'Parent',axes1,'CDataMapping','scaled')
set(axes1,'handlevisibility','off','visible','off');
axes2 = axes('Parent',figure1);
view( axes2,[15 30])
box(axes2,'on')
set(axes2,'Color','none','LineWidth',6,'TickLength',[0.01 0.005],'XColor',[1 0 1],'YColor',[1 0 1],'ZColor',[1 0 1]);
verticalaxis = line([0,0],[0,0],[0,min(z)],'color','k','linestyle','-','linewidth',1.5); % 竖轴
pendulumline = line([0,y(1,1)],[0,y(3,1)],[0,z(1,1)],'color','y','linestyle', '-','linewidth',2); % 摆线初始化
pendulumbob = line(y(1,1),y(3,1),z(1,1),'color','r','marker','.', 'markersize',25); % 摆锤初始化
pendulumtra = line([y(1,1),y(1,1)],[y(3,1),y(3,1)],[z(1,1),z(1,1)]); % 摆锤起点
view([15 30])
box on
%axis([-0.6 0.6 -1 0.2]);
axis ([min(y(1,:)) max(y(1,:)) min(y(3,:)) max(y(3,:)) min(z) 0])
hold on
gg=1
for i = 1:size(t,2);
set(pendulumline,'xdata',[0,y(1,i)],'ydata',[0,y(3,i)],'zdata',[0,z(1,i)]); % 更新摆线位置
set(pendulumbob,'xdata',y(1,i),'ydata',y(3,i),'zdata',z(1,i));% 更新摆锤位置
M(i,:)=[y(1,i) y(3,i) z(1,i)];
set(pendulumtra,'xdata',y(1,1:i),'ydata',y(3,1:i),'zdata',z(1,1:i))
fmat(:,i)=getframe; %拷贝祯到矩阵fmat中
im = frame2im(getframe);
[I,map] = rgb2ind(im,256);
if gg == 1
imwrite(I,map,'M1.gif','GIF', 'Loopcount',inf,'DelayTime',0.05);
gg = gg 1;
else
imwrite(I,map,'M1.gif','WriteMode','append','DelayTime',0.05);
end
end这样我们就获得了二种求解x,y方法，符号运算法，数值解法，数值法好处是一次性给出的x,y，以及对应偏导值dx、dy，通常表示物体运动速度，这在运功方程中是很重要的，如果你求出来是解析解，需要求dx,dy,就需要做一次求导，我们比较一下三种方法对应数值差别.
出乎意料的数值计算法和符号法得到的解竟然不完全一致。这就说明即使是 Runge-Kutta方法也是存在误差的，t:0~30 区间x重合度较高，y有偏差，t:30~60 y 重合度高，x差别大，过冷水需要进一步摸索，本期文章给了给出较为典型的求解微分方程组的方法，应用这种中方法求解常见微分问题基本是没有问题的，希望对有需要的读者快应用起来。
过冷水和你分享 matlab读取存储各种文件的方法 文末有独家金曲分享
