matlab 含阻力单摆微分方程可视化

最新推荐文章于 2024-08-30 20:19:02 发布

slandarer

最新推荐文章于 2024-08-30 20:19:02 发布

阅读量2.4k

点赞数 3

分类专栏： # MATLAB 笔记文章标签： MATLAB 单摆常微分可视化

本文链接：https://blog.csdn.net/slandarer/article/details/104429297

版权

MATLAB 笔记专栏收录该内容

46 篇文章

订阅专栏

文章目录:

1.简介

这是一期将单摆微分方程可视化的博文，我们都知道单摆常微分方程求解过程中会涉及到椭圆积分，难以用常见函数表示其结果，所以我们这篇博文想办法将其数值解可视化。

2.效果

在这里插入图片描述
其中横坐标为： $\theta(t),\theta(t)\in[-2.5\pi,2.5\pi]$
纵坐标为： $\dot{\theta}(t),\dot{\theta}(t)\in[-2\pi,2\pi]$

3.基本步骤

3.1 解方程

首先我们有微分方程：
$\ddot{\theta}(t)=-\mu\dot{\theta}(t)-\frac g Lsin(\theta(t))$
我们将其变形得到微分方程组：
$\frac d {dt} \begin{bmatrix} \theta(t)\\\dot{\theta}(t) \end{bmatrix}= \begin{bmatrix} \dot{\theta}(t)\\\ddot{\theta}(t) \end{bmatrix}= \begin{bmatrix} \dot{\theta}(t)\\-\mu\dot{\theta}(t)-\frac g Lsin(\theta(t)) \end{bmatrix}$
对此我们均匀的在平面上取点，并计算不同 $\begin{bmatrix} \theta(t)\\\dot{\theta}(t) \end{bmatrix}$ 对应的 $\frac d {dt} \begin{bmatrix} \theta(t)\\\dot{\theta}(t) \end{bmatrix}$ 就可完成微分方程的可视化。

3.2 背景绘制

如果你看过我之前写的小游戏，可会下意识的这样写：

axis([-2.5,2.5,-2,2].*pi)
set(gca,‘color’,[0 0 0.0078])
set(gca,‘xtick’,[],‘ytick’,[],‘xcolor’,‘w’,‘ycolor’,‘w’)

但这样写后的效果是这个样子的：
在这里插入图片描述
很明显不够高端大气，至少不如效果图里显得大气
我们不妨重新定义一个大的figure，并且重新设置axes大小：

penduium.fig=figure(‘units’,‘pixels’,…
‘position’,[350 100 800 500],…
‘Numbertitle’,‘off’,…
‘menubar’,‘none’,…
‘resize’,‘off’,…
‘name’,‘simple_penduium’,…
‘color’,[0.95 0.95 0.95]);
axes(‘position’,[0 0 1 1])
axis([-2.5,2.5,-2,2].*pi)
set(gca,‘color’,[0 0 0.0078])
set(gca,‘xtick’,[],‘ytick’,[],‘xcolor’,‘w’,‘ycolor’,‘w’)

像这样的到的背景长这样：
在这里插入图片描述
没有边框帅气很多，之后我们(残暴的）为其增添坐标系：

hold on
plot([1;1]*(-2.5:0.1:2.5).*pi,[-2;2]*ones(1,51).*pi,'color',[0.1569 0.2039 0.1882],'linewidth',0.01)
plot([-2.5;2.5]*ones(1,41).*pi,[1;1]*(-2:0.1:2).*pi,'color',[0.1569 0.2039 0.1882],'linewidth',0.01)
plot([-2.5;2.5].*pi,[0 0],'color',[0.9373 0.9569 0.9569],'linewidth',1)
plot([0 0],[-2;2].*pi,'color',[0.9373 0.9569 0.9569],'linewidth',1)
plot([-2.5;2.5]*ones(1,6).*pi,[1;1]*[1.5,1,0.5,-0.5,-1,-1.5].*pi,'color',[0.2353 0.4235 0.4745],'linewidth',0.8)
plot([1;1]*[2,1.5,1,0.5,-0.5,-1,-1.5 -2].*pi,[-2;2]*ones(1,8).*pi,'color',[0.2353 0.4235 0.4745],'linewidth',0.8)

就得到了效果图中的背景：
在这里插入图片描述

3.3计算并绘图

只需要写一个双层for循环计算方向向量，将其归一化后画在图中，并且根据不同的长度为其赋予不同的颜色就好了

for i=(-2.5+1/6:1/6:2.5-1/6).*pi
    for j=(-2+1/4:1/4:2-1/4).*pi
        a=f([i,j],M,g,L);
        len=norm(a);
        a=a/len;
        a=a.*vector_len;
        Spoint=pos_exchange([i,j]);
        Epoint=pos_exchange([i,j]+a);
        if ~any(isnan(a))
            annotation('arrow',[Spoint(1),Epoint(1)],[Spoint(2),Epoint(2)],...
            'color',color_exchange(len,max_len,color_map),...
            'linewidth',0.2)
        end
    end
end

其中f,pos_exchange,color_exchange三个匿名函数为：

%========================================================= 
 f=@(theta,M,g,L) [theta(2),-M*theta(2)-g/L*sin(theta(1))];
%theta(1) ->theta
%theta(2) ->theta'
%==============================================
 pos_exchange=@(pos) pos./[pi*5,pi*4]+[0.5,0.5];
%==============================================
 color_map=[0.9725    0.3804    0.3529
            0.9725    0.3804    0.3529
            0.9020    0.3843    0.3765
            0.9020    0.3843    0.3765
            0.9333    0.4118    0.3765
            0.9686    0.5804    0.2235
            0.9765    0.8392    0.1098
            0.9882    0.9804    0.0588
            0.7961    0.8353    0.2510
            0.6510    0.7373    0.2000
            0.5961    0.7059    0.4039];
 color_exchange=@(value,maxvalue,color_map) color_map(12-(floor(value/(maxvalue/11))+1),:);
%=========================================================================================

3.4完整代码

function simple_penduium
M=1.4;
g=9.8;
L=2;
vector_len=0.75;
%========================================================= 
 f=@(theta,M,g,L) [theta(2),-M*theta(2)-g/L*sin(theta(1))];
%theta(1) ->theta
%theta(2) ->theta'
%==============================================
 pos_exchange=@(pos) pos./[pi*5,pi*4]+[0.5,0.5];
%==============================================
 color_map=[0.9725    0.3804    0.3529
            0.9725    0.3804    0.3529
            0.9020    0.3843    0.3765
            0.9020    0.3843    0.3765
            0.9333    0.4118    0.3765
            0.9686    0.5804    0.2235
            0.9765    0.8392    0.1098
            0.9882    0.9804    0.0588
            0.7961    0.8353    0.2510
            0.6510    0.7373    0.2000
            0.5961    0.7059    0.4039];
 color_exchange=@(value,maxvalue,color_map) color_map(12-(floor(value/(maxvalue/11))+1),:);
%=========================================================================================
max_len=norm(f([pi/2,2*pi],M,g,L));
global penduium
penduium.fig=figure('units','pixels',...
    'position',[350 100 800 500],...
    'Numbertitle','off',...
    'menubar','none',...
    'resize','off',...
    'name','simple_penduium',...
    'color',[0.95 0.95 0.95]);
axes('position',[0 0 1 1])
axis([-2.5,2.5,-2,2].*pi)
set(gca,'color',[0 0 0.0078])
set(gca,'xtick',[],'ytick',[],'xcolor','w','ycolor','w')
hold on


plot([1;1]*(-2.5:0.1:2.5).*pi,[-2;2]*ones(1,51).*pi,'color',[0.1569 0.2039 0.1882],'linewidth',0.01)
plot([-2.5;2.5]*ones(1,41).*pi,[1;1]*(-2:0.1:2).*pi,'color',[0.1569 0.2039 0.1882],'linewidth',0.01)
plot([-2.5;2.5].*pi,[0 0],'color',[0.9373 0.9569 0.9569],'linewidth',1)
plot([0 0],[-2;2].*pi,'color',[0.9373 0.9569 0.9569],'linewidth',1)
plot([-2.5;2.5]*ones(1,6).*pi,[1;1]*[1.5,1,0.5,-0.5,-1,-1.5].*pi,'color',[0.2353 0.4235 0.4745],'linewidth',0.8)
plot([1;1]*[2,1.5,1,0.5,-0.5,-1,-1.5 -2].*pi,[-2;2]*ones(1,8).*pi,'color',[0.2353 0.4235 0.4745],'linewidth',0.8)



for i=(-2.5+1/6:1/6:2.5-1/6).*pi
    for j=(-2+1/4:1/4:2-1/4).*pi
        a=f([i,j],M,g,L);
        len=norm(a);
        a=a/len;
        a=a.*vector_len;
        Spoint=pos_exchange([i,j]);
        Epoint=pos_exchange([i,j]+a);
        if ~any(isnan(a))
            annotation('arrow',[Spoint(1),Epoint(1)],[Spoint(2),Epoint(2)],...
            'color',color_exchange(len,max_len,color_map),...
            'linewidth',0.2)
        end
    end
end
%text(pi,-0.1,'\pi','fontsize',20,'color','w','HorizontalAlignment', 'center')
%text(-pi,-0.1,['-','\pi'],'fontsize',20,'color','w','HorizontalAlignment', 'center')
%text(2*pi,-0.1,'  \pi','fontsize',20,'color','w','HorizontalAlignment', 'center')
%text(2*pi-0.1,-0.18,'2','fontsize',13,'color','w','HorizontalAlignment', 'center')
%text(-2*pi,-0.1,'- \pi','fontsize',20,'color','w','HorizontalAlignment', 'center')
%text(-2*pi-0.08,-0.18,'2','fontsize',13,'color','w','HorizontalAlignment', 'center')
%text(2.5*pi-0.3,0.5,'\theta','fontsize',25,'color','w','HorizontalAlignment', 'center')
%text(0.5,2*pi-0.3,['\theta','’'],'fontsize',25,'color','w','HorizontalAlignment', 'center')
end