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

文章目录:1.简介2.效果3.基本步骤3.1 解方程3.2 背景绘制3.3计算并绘图3.4完整代码1.简介这是一期将单摆微分方程可视化的博文，我们都知道单摆常微分方程求解过程中会涉及到椭圆积分，难以用常见函数表示其结果，所以我们这篇博文想办法将其数值解可视化。2.效果其中横坐标为： θ(t),θ(t)∈[−2.5π,2.5π]\theta(t),\theta(t)\in[-2.5\pi,...

slandarer

2734人浏览 · 2020-02-21 18:35:55

slandarer · 2020-02-21 18:35:55 发布

文章目录:

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

魔乐社区

魔乐社区（Modelers.cn) 是一个中立、公益的人工智能社区，提供人工智能工具、模型、数据的托管、展示与应用协同服务，为人工智能开发及爱好者搭建开放的学习交流平台。社区通过理事会方式运作，由全产业链共同建设、共同运营、共同享有，推动国产AI生态繁荣发展。

更多推荐

全家桶集齐！Qwen3.5四款小模型上线魔乐社区，附昇腾全套实践教程

魔乐社区

Pont - 搭建前后端之桥：高效、灵活的接口管理工具

Pont 是一款强大的数据服务层解决方案，它能够帮助开发者快速搭建前后端之间的桥梁，实现接口的高效管理和代码自动生成。无论是新手还是有经验的开发者，都能通过 Pont 轻松处理接口文档、生成类型安全的 API 代码，从而显著提升开发效率。[![Pont 工具标志](https://raw.gitcode.com/gh_mirrors/po/pont/raw/3f1b7d4bbba3fd2dda

魔乐社区

如何快速上手 hvac：HashiCorp Vault Python 客户端零基础入门指南

**hvac** 是 HashiCorp Vault 的 Python 3.X 客户端库，专为开发者提供简单高效的 Vault 交互方式。无论你是需要管理密钥、配置身份验证，还是实现安全的秘密数据存储，hvac 都能帮助你轻松搞定 Vault 的各项操作。本文将带你零基础快速入门，从安装到基础操作，让你在几分钟内即可上手使用这个强大的工具。[![hvac 客户端 Logo](https://r