matlab动力学分析程序详解

合集下载

中心差分 动力学方程 matlab

中心差分 动力学方程 matlab

中心差分动力学方程matlab随着科学技术的不断发展,数值计算方法在各个领域得到了广泛应用。

中心差分法作为一种常用的数值计算方法,在科学研究和工程计算中具有重要意义。

本文将介绍中心差分法、动力学方程以及MATLAB在二者中的应用,并通过实例演示来展示具体操作。

一、中心差分法简介中心差分法是一种在离散点上计算连续函数的导数、微分方程数值解的数值方法。

它的基本思想是将函数在一个区间内离散化为若干个点,然后用差分的形式表示函数的导数,再通过求和得到整个区间上的积分。

中心差分法具有计算简便、精度较高、稳定性较好等特点。

二、动力学方程简介动力学方程是描述物体运动规律的数学方程,它可以用于分析各种物理、化学和生物过程。

动力学方程一般包括牛顿方程、拉格朗日方程、哈密顿方程等。

求解动力学方程可以得到物体的运动轨迹、速度、加速度等物理量,对于工程设计和控制系统分析具有重要意义。

三、MATLAB在中心差分法和动力学方程中的应用MATLAB是一款功能强大的数学软件,可以方便地实现中心差分法和动力学方程的数值求解。

以下将以一个简单的例子说明MATLAB在中心差分法和动力学方程中的应用。

假设我们要计算以下动力学方程的数值解:dx/dt = -ax + b初始条件:x(0) = 1,t(0) = 01.首先,编写MATLAB脚本,设置时间步长Δt和仿真时间t_end:```MATLABt_end = 10; % 仿真时间dt = 0.1; % 时间步长```2.接着,编写中心差分法求解动力学方程的代码:```MATLABx = zeros(1, t_end/dt + 1); % 初始化x arrayt = zeros(1, t_end/dt + 1); % 初始化t arrayx(1) = 1; % 设置初始条件for i = 2:t_end/dt + 1x(i) = x(i-1) + dt * (-a * x(i-1) + b);t(i) = i * dt;end```3.绘制数值解曲线:```MATLABfigure;plot(t, x);xlabel("t");ylabel("x");title("动力学方程的数值解");```四、结论与建议本文通过一个简单实例展示了MATLAB在中心差分法和动力学方程中的应用。

matlab 系统动力学模型

matlab 系统动力学模型

matlab 系统动力学模型Matlab 系统动力学模型引言:系统动力学是研究动态系统行为的数学方法,通过描述和分析系统在时间上的演化过程,揭示系统内部的关系和相互作用规律。

Matlab是一种强大的数值计算软件,广泛应用于系统动力学模型的建立和仿真。

本文将介绍Matlab在系统动力学建模中的应用,并结合实例进行说明。

一、系统动力学基本概念系统动力学是一种描述系统行为的数学工具,它将系统划分为不同的部分,并研究它们之间的相互作用。

系统动力学模型通常由一组关于系统部分之间关系的微分方程或差分方程组成。

在建立模型时,需要考虑系统的输入、输出以及系统内部的状态变量,并通过数学表达式描述它们之间的关系。

二、Matlab在系统动力学模型中的应用Matlab提供了丰富的数学函数和工具箱,使得系统动力学模型的建立和仿真变得更加简单和高效。

下面将以一个简单的例子来说明Matlab在系统动力学建模中的应用。

假设有一个简单的机械系统,由弹簧和质量构成。

假设弹簧的刚度为k,质量为m,阻尼系数为b。

我们想要建立一个系统动力学模型,来描述质点的运动过程。

我们需要确定系统的状态变量和输入输出。

在这个例子中,质点的位移x是系统的状态变量,外力F是系统的输入,质点的加速度a 是系统的输出。

根据牛顿第二定律,我们可以建立如下的微分方程:m * a = F - b * v - k * x其中,v是质点的速度。

为了建立系统的动力学模型,我们需要对该微分方程进行求解。

在Matlab中,可以使用ode45函数来解决常微分方程。

具体的Matlab代码如下:```matlabfunction dxdt = system_dynamics(t, x)m = 1;k = 10;b = 0.5;F = 5;v = x(2);a = (F -b * v - k * x(1)) / m;dxdt = [v; a];end[t, x] = ode45(@system_dynamics, [0, 10], [0, 0]);plot(t, x(:, 1));xlabel('Time');ylabel('Displacement');title('System Dynamics');```在上述代码中,system_dynamics函数定义了系统的微分方程,其中包括质点的加速度和速度的计算。

基于MATLAB的六杆机构动力学分析与仿真

基于MATLAB的六杆机构动力学分析与仿真

六杆机构的动力学分析仿真一系统模型建立为了对机构进行仿真分析,首先必须建立机构数学模型,即位置方程,然后利用MATLAB仿真分析工具箱Simulink对其进行仿真分析。

图3.24所示是由原动件(曲柄1)和RRR—RRP六杆机构。

各构件的尺寸为r1=400mm,r2=1200mm,r3=800mm,r4=1500mm,r5=1200mm;各构件的质心为rc1=200mm,rc2=600mm,rc3=400mm,rc5=600mm;质量为m1=1.2kg,m2=3kg,m3=2.2kg;m5=3.6kg,m6=6kg; 转动惯量为J1=0.016kg·m2,J2=0.25kg·m2;J3=0.09kg·m2,J5=0.45kg·m2;构件6的工作阻力F6=1000N,其他构件所受外力和外力矩均为零,构件1以等角速度10 rad/s逆时针方向回转,试求不计摩擦时,转动副A的约束反力、驱动力矩、移动副F的约束反力。

图1-1此机构模型可以分为曲柄的动力学、RRR II级杆组的动力学和RRP II级杆组的动力学,再分别对这三个模型进行相应参数的求解。

图1-2 AB 构件受力模型如上图1-2对于曲柄AB 由理论力学可以列出表达式:111XA Re R ••=+-s m F R X XB 111y A Im R ••=+-s m F R y yB1111111111111cos )(sin )(cos sin ••=---+-++θθθθθJ r r R r r R r R r R M M c yB c XB c yA c XA F由运动学知识可以推得:)cos()2/cos(Re Re 12111111πθθπθθ++++=•••••••c c r r A s )sin()2/sin(Im Im 12111111πθθπθθ++++=•••••••c c r r A s将上述各式合并成矩阵形式有,⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡++-+++++-++++=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡••••••••••g m R F r m r m A m R F r m r m A m M R R yB y c c XBX c c yA XA 111211111111112111111111)sin()2/sin(Im )cos()2/cos(Re πθθπθθπθθπθθ(1-21)如图1-3,对构件BC 的约束反力推导如下,图1-3 BC 构件受力模型222Re ••=++s m R F R XC X XB 2222Im ••=-++s m g m R F R yC y yB2222222222222cos )(sin )(cos sin ••=-----+θθθθθJ r r R r r R r R r R M c yC c XC c yB c XB如图1-4,对构件BC 的约束反力推导如下,图 1-4 CD 构件受力模型333Re ••=-+s m R F R XC X XD 3333Im ••=-++s m g m R F R yC y yD3333333333333cos )(sin )(cos sin ••=-----+θθθθθJ r r R r r R r R r R M c yC c XC c yD c XD由运动学可以推导得,)sin()2/sin(Im Im 22222222πθθπθθ++++=•••••••c c r r B s )cos()2/cos(Re Re 22222222πθθπθθ++++=•••••••c c r r B s )cos()2/cos(Re Re 32333333πθθπθθ++++=•••••••c c r r D s )sin()2/sin(Im Im 32333333πθθπθθ++++=•••••••c c r r D s将上述BC 构件,CD 构件各式合并成矩阵形式有,⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----------33333333332222222222cos sin cos )(sin )(0010100001010000cos )(sin )(cos sin 001010000101θθθθθθθθc c c c c c c c r r r r r r r r r r r r ⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡yD XD yC XC yB XB R R R R R R =⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-+-++++-++++-+-++++-++++••••••••••••••••••••••••3333332333333333323333333322222222222222222222222222)sin()2/sin(Im )cos()2/cos(Re )sin()2/sin(Im )cos()2/cos(Re M J g m F r m r m D m F r m r m D m M J g m F r m r m B m F r m r m B m y c c X c c y c c X c c θπθθπθθπθθπθθθπθθπθθπθθπθθ (1-22)如图1-5 对构件5进行约束反力的推导如下,图1-5 CE 杆件受力模型••=++s m R F R xE x xC Re 55 ••=-++s m g m R F R yE y yC Im 5555555555555555cos )(sin )(cos sin ••=-+-+--θθθθθJ r r R r r R r R r R M c yE c xE c yC c xC如图1-6 对滑块进行受力分析如下,滑块受力模型••=--E m R R F F xE x Re sin 666θ ••=-+-E m g m R R F F yE y Im cos 6666θ由运动学可推,)cos()2/cos(C Re Re 5255555πθθπθθ•••••••+++=c c r r s )sin()2/sin(C Re Im 5255555πθθπθθ•••••••+++=c c r r s66cos Re θ••••=s E 66sin Im θ••••=s E⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡+---+-++++-++++=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-------••••••••••••••••g m F s m F s m M J g m F r m r m m F r m r m m R R R R R r r r r r r y x y c c x c c F yE xE yC xC c c c c 66666666655555525555555555255555555665555555555sin cos )sin()2/sin(C Re )cos()2/cos(C Re cos 1000sin 01000cos )(sin )(cos sin 010*******θθθπθθπθθπθθπθθθθθθθθ(1-23)二编程与仿真利用MATLAB进行仿真分析,主要包括两个步骤:首先是编制计算所需要的函数模块,然后利用其仿真工具箱Simulink建立仿真系统框图,设定初始参数进行仿真分析。

sindy的matlab程序-概述说明以及解释

sindy的matlab程序-概述说明以及解释

sindy的matlab程序-概述说明以及解释1.引言1.1 概述Sindy是一种基于数据驱动的系统辨识方法,通过对系统的动态行为进行分析和建模,可以帮助我们更好地理解系统的运行机制和规律。

Matlab作为一种强大的科学计算工具,能够提供丰富的功能和工具,帮助我们进行数据处理、模型建立和结果分析。

本文将详细介绍Sindy在Matlab环境下的应用,探讨其在不同领域中的作用和价值。

通过对Sindy程序的优势和局限性进行分析,可以更全面地了解其在系统辨识方面的特点和适用范围。

最后,我们将总结Sindy 的Matlab程序的重要性,展望其未来在系统辨识领域的发展,并希望能为相关研究提供一定的参考和启发。

1.2 文章结构本篇文章主要分为三个部分:引言、正文和结论。

在引言部分,将会对文章的主题进行一定的概述,介绍Sindy的Matlab程序的背景和意义,以及对文章的结构进行简要的介绍。

在正文部分,将详细介绍Sindy的Matlab程序的相关内容,包括程序的介绍、应用领域、优势和局限性等方面。

最后,在结论部分,将总结Sindy的Matlab程序的重要性,展望其在未来的发展,并给出一些结束语,为全文画上一个完美的句号。

1.3 目的本文的目的是介绍Sindy的Matlab程序,探讨其在科学研究和工程领域中的应用情况。

通过对Sindy程序的介绍和分析,读者可以更深入地了解Sindy程序的原理和特点,以及其在系统辨识、动力系统建模等方面的重要性和价值。

同时,本文也将讨论Sindy程序的优势和局限性,对于读者在选择合适的程序工具时提供参考。

通过本文的阐述,旨在激发读者对于Sindy程序的兴趣,促进该程序在未来的发展和应用。

2.正文2.1 Sindy的Matlab程序介绍Sindy是一个用于系统辨识和模型推断的Matlab程序。

该程序的全称为Sparse Identification of Nonlinear Dynamics,意为稀疏非线性动力学识别。

matlab求解动力学微分方程

matlab求解动力学微分方程

【文章标题】:深度探究:使用Matlab求解动力学微分方程一、探索动力学微分方程在物理学、工程学和生物学等领域中,动力学微分方程广泛应用于描述系统随时间演化的规律。

动力学系统的行为可以通过微分方程来捕捉,而解析求解这些微分方程通常是困难的。

数值方法成为了研究这些系统行为的重要工具。

二、Matlab在动力学微分方程求解中的应用Matlab作为一种强大的数学计算软件,提供了丰富的工具和函数用于求解微分方程。

其拥有的ode45、ode23等求解器能够高效地求解多种微分方程,并且具有较好的稳定性和精度。

三、从简到繁:使用Matlab求解一阶动力学微分方程我们可以从一阶动力学微分方程开始探讨。

考虑如下的一阶常微分方程:\[ \frac{dy}{dt} = f(t, y) \]我们可以使用Matlab的ode45函数来求解此类微分方程。

通过传入方程的右侧函数f(t, y)、初始条件和求解的时间范围,ode45可以给出相应的数值解。

四、由浅入深:使用Matlab求解高阶动力学微分方程我们可以深入探讨高阶动力学微分方程的求解。

考虑如下的二阶常微分方程:\[ \frac{d^2y}{dt^2} = g(t, y, \frac{dy}{dt}) \]同样地,Matlab的ode45函数也可以用于求解高阶微分方程。

我们需要将这个二阶微分方程化为一个一阶方程组的形式,然后传入ode45函数进行求解。

五、总结回顾:深刻理解动力学微分方程的数值求解通过本文的讨论,我们深入了解了如何使用Matlab求解动力学微分方程。

从简单的一阶微分方程到复杂的高阶微分方程,Matlab都提供了强大的求解工具。

通过数值求解,我们可以得到系统随时间演化的行为规律,从而加深我们对动力学系统的理解。

六、个人观点与理解作为文章写手,我个人认为Matlab在求解动力学微分方程方面具有非常大的优势。

其丰富的函数库和高效的求解器为研究人员提供了便利,使他们能够更深入地探究动力学系统的行为。

matlab newmark法

matlab newmark法

matlab newmark法Matlab Newmark法是一种非线性动力学分析方法,主要用于求解动力学系统的时间响应。

该方法由Newmark在20世纪50年代提出,在工程结构领域得到了广泛应用。

本文将分步骤回答关于Matlab Newmark法的问题,包括算法原理、计算步骤、优缺点以及实际案例的应用。

一、算法原理1.1 基本原理Matlab Newmark法是一种基于离散时间步长的计算方法。

其基本原理是通过将系统的运动方程转化为等效的一阶微分方程组,然后使用步进法进行数值求解。

该方法采用了二阶精度的数值积分公式,具有较高的计算精度和稳定性。

1.2 新马克法公式Matlab Newmark法的核心公式为:δu(t+Δt) = u(t) + Δt * v(t) + Δt^2 * (0.5 - β) * a(t)δv(t+Δt) = v(t) + Δt * (1 - γ) * a(t)δa(t+Δt) = (1 - γ) * a(t) + γ* a(t+Δt)其中,δ表示增量,u(t)、v(t)和a(t)分别表示位移、速度和加速度在时间t的值,β和γ为Newmark法的两个参数。

二、计算步骤2.1 确定系统参数首先,需要确定系统的质量矩阵、刚度矩阵和阻尼矩阵,以及外部激励载荷等参数。

2.2 确定时间步长根据求解精度和计算效率的要求,选择合适的时间步长Δt。

2.3 初始化位移、速度和加速度给定初始位移、速度和加速度的值。

2.4 进行时间循环使用Newmark法的公式,根据当前时刻的位移、速度和加速度的值,计算下一时刻的位移、速度和加速度。

2.5 判断收敛条件在每个时间步长内,判断计算结果是否满足收敛要求。

如果满足要求,则继续计算下一个时间步长;如果不满足要求,则重新选择适当的步长,并重新进行计算。

2.6 输出结果将每个时间步长内计算得到的位移、速度和加速度的值保存起来,以获取系统的时间响应曲线。

三、优缺点3.1 优点Matlab Newmark法具有以下优点:- 可以处理复杂的非线性动力学系统。

第6章Matlab应用之动力学与振动

第6章Matlab应用之动力学与振动
6.2.2 线性系统的自由振动
一.运动微分方程 当
0, F ( ) 0 时,得到线性振动系统的自由振动方程。
d x dx 2 x0 2 d d
2
上一页
目录
返回
下一页
13
6.2 单自由度系统
二.MATLAB求解 编写方程对应的函数文件FreeOscillation.m function xdot=FreeOscillation(t,x,zeta,Alpha) xdot=[x(2);-2.0*zeta*x(2)-x(1)-Alpha*x(1)^3]; end
上一页
目录
返回
下一页
21
6.2 单自由度系统
续上: figure(1) xlabel('\tau'); ylabel('x(\tau)'); axis([0.0,30.0,-3.0,3.0]); legend(d(1,:),d(2,:),d(3,:)); figure(2) xlabel('x(\tau)'); ylabel('dx/d\tau'); axis([-2.0,3.0,-2.0,3.0]); legend(d(1,:),d(2,:),d(3,:));
dx1 x2 d dx2 x1 d signum( x2 ) d
d x dx x d signum( ) 2 d d
上一页 目录 返回 下一页
25
6.2 单自由度系统
2、Matlab求解
编写常微分方程对应的函数文件FrictionOscillation.m
function xdot=FrictionOscillation(t,x,d) % 非线性阻尼系统ode文件 if abs(x(1))<=d && x(2)==0.0; xdot=[0;0]; else xdot=[x(2);-d*sign(x(2))-x(1)]; end

matlab使用拉格朗日法计算动力学方程参数

matlab使用拉格朗日法计算动力学方程参数

matlab使用拉格朗日法计算动力学方程参数标题:用MATLAB使用拉格朗日法计算动力学方程参数导言:在工程和科学领域中,我们经常需要对系统进行动力学建模和分析。

动力学方程是描述系统运动的数学表达式,它们可以帮助我们理解和预测系统的行为。

而计算动力学方程的参数对于系统的设计、优化和控制具有重要意义。

本文将介绍如何使用MATLAB编程语言和拉格朗日法来计算动力学方程参数,通过数值求解和优化方法,为工程师和科学家提供一个有力的工具。

1. 动力学方程与参数计算在动力学中,系统可以通过一组微分方程来描述。

这些方程表示系统的行为和演变,通常包括质量、速度、加速度和其他相关因素。

对于复杂系统,计算动力学参数可能是一项繁琐且复杂的任务。

幸运的是,拉格朗日法可以简化这一过程,通过定义能量或运动方程来推导系统动力学方程。

2. 拉格朗日法介绍拉格朗日法是一种基于能量和运动方程的变分法。

通过将系统的动能和势能组合成拉格朗日函数,并使用欧拉-拉格朗日方程,可以得到系统的运动方程。

拉格朗日法不仅适用于具有广义坐标的连续系统,还适用于离散系统和有束缚条件的系统。

3. MATLAB编程实现在MATLAB中,我们可以使用符号计算工具包来处理复杂的数学表达式和符号计算。

我们需要定义系统的广义坐标、广义速度和拉格朗日函数。

通过欧拉-拉格朗日方程生成系统的运动方程。

我们可以使用数值求解和优化方法来计算动力学参数。

4. 示例:双摆系统为了更好地理解如何使用MATLAB实现拉格朗日法,我们将以一个双摆系统为例。

双摆系统由两个连杆组成,每个连杆上挂有一个质点。

我们需要计算系统的动力学参数,包括质点的位置、速度和加速度。

通过拉格朗日法,我们可以推导出系统的运动方程,并使用MATLAB 进行求解和优化。

5. 讨论与总结本文介绍了使用MATLAB编程语言和拉格朗日法计算动力学方程参数的方法。

通过拉格朗日法,我们可以简化动力学参数的计算,并为系统的设计和优化提供可行性的解决方案。

matlab动力学分析程序详解

matlab动力学分析程序详解

matlab动力学分析程序详解··1.微分方程的定义对于duffing 方程032=++x x xω ,先将方程写作--==3112221x x x x x ω function dy=duffing(t,x) omega=1;%定义参数f1=x(2);f2=-omega^2*x(1)-x(1)^3; dy=[f1;f2];2.微分方程的求解function solve (tstop) tstop=500;%定义时间长度 y0=[0.01;0];%定义初始条件[t,y]=ode45('duffing',tstop,y0,[]);function solve (tstop) step=0.01;%定义步长y0=rand(1,2);%随机初始条件tspan=[0:step:500];%定义时间范围[t,y]=ode45('duffing',tspan,y0);3.时间历程的绘制时间历程横轴为t ,纵轴为y ,绘制时只取稳态部分。

plot(t,y(:,1));%绘制y 的时间历程 xlabel('t')%横轴为t ylabel('y')%纵轴为y grid;%显示网格线axis([460 500 -Inf Inf])%图形显示范围设置4.相图的绘制相图的横轴为y ,纵轴为dy/dt ,绘制时也只取稳态部分。

红色部··分表示只取最后1000个点。

plot(y(end-1000:end,1),y(end-1000:end,2));%绘制y 的时间历程xlabel('y')%横轴为yylabel('dy/dt')%纵轴为dy/dt grid;%显示网格线5.Poincare 映射的绘制对于不同的系统,Poincare 截面的选取方法也不同对于自治系统一般每过其对应线性系统的固有周期,截取一次对于非自治系统,一般每过其激励的周期,截取一次例程:duffing 方程032=++x x x ω的poincare 映射function poincare(tstop) global omega; omega=1;T=2*pi/omega;%线性系统的周期或激励的周期step=T/100;%定义步长为T/100 y0=[0.01;0];%初始条件tspan=[0:step:100*T];%定义时间范围[t,y]=ode45('duffing',tspan,y0);for i=5000:100:10000%稳态过程每个周期取一个点plot(y(i,1),y(i,2),'b.'); hold on;% 保留上一次的图形 endxlabel('y');ylabel('dy/dt');Poincare 映射也可以通过取极值点得到 function poincare(tstop) y0=[0.01;0];tspan=[0:0.01:500];[t,y]=ode45('duffing',tspan,y0); count=find(t>100);%截取稳态过程 y=y(count,:);n=length(y(:,1));%计算点的总数··for i=2:n-1if y(i-1,1)+epsy(i+1,1)+eps % 简单的取出局部最大值plot(y(i,1),y(i,2),'.'); hold on end endxlabel('y');ylabel('dy/dt');6.频谱yy=fft(y(end-1000:end,1)); N=length(yy); power=abs(yy);freq=(1:N-1)*1/step/N;plot(freq(1:N/2),power(1:N/2)); xlabel('f(y)') ylabel('y')7.算例duffing 方程03=++x x x的时间历程,相图,频谱和poincare 映射。

机械原理4-23MATLAB平面连杆机构动力学分析

机械原理4-23MATLAB平面连杆机构动力学分析

基于MATLAB/Solidworks COSMOSMotion的平面连杆机构动力学分析07208517王锡霖4-23在图示的正弦机构中,已知l AB =100 mm,h1=120 mm,h2 =80 mm,W1 =10 rad/s(常数),滑块2和构件3的重量分别为G2 =40 N和G3 =100 N,质心S2 和S3 的位置如图所示,加于构件3上的生产阻力Fr=400 N,构件1的重力和惯性力略去不计。

试用解析法求机构在Φ1=60°、150°、220°位置时各运动副反力和需加于。

构件1上的平衡力偶Mb分别对三个构件进行受力分析如图:构件3受力图构件2受力图构件1受力图(1)滑块2:V S2 =L AB W1 ①a s2 = L AB W12②构件3:S=L AB sinΦ1 ③V3=L AB W1 COSΦ1 ④a3=-L AB W12 sinΦ1 ⑤(2)确定惯性力:F12=m2as2=(G2/g)LABW12 ⑥F13=m3a3=(G3/g)LABW12sinΦ1 ⑦(3)各构件的平衡方程:构件3:∑Fy=0,FR23 =Fr-F13∑Fx=0,FR4’=FR4∑MS3 =0,FR4=FR23LAcosΦ1/h2构件2:∑Fx=0,FR12x=F12cosΦ1∑Fy=0,FR12y=FR32-F12sinΦ1构件1:∑Fx=0,FR41x=FR12x∑Fy=0,FR41y=FR12y∑MA =0,Mb=FR32LABcosΦ1总共有八个方程,八个未知数。

归纳出一元八次方程矩阵:1 0 0 0 0 0 0 0 FR23 Fr-F130 1 -1 0 0 0 0 0 FR4’ 0-LAB COSΦ1/h20 1 0 0 0 0 0 FR40 0 0 1 0 0 0 0 FR12x = F12cosΦ1-1 0 0 0 1 0 0 0 FR12y -F12sinΦ10 0 0 -1 0 1 0 0 FR41x 00 0 0 0 -1 0 1 0 FR41y 0-LABCOSΦ1 0 0 0 0 0 0 1 Mb 0 AX=B进而可得:X=A\B。

第6章 Matlab应用之动力学与振动概要

第6章  Matlab应用之动力学与振动概要

习题
上一页
目录
返回
下一页
3
6.1 轨迹
举例说明:重力场中有两个物体,其中质量为m2的物体固定,而 质量为m1的物体绕m2做平面圆周运动.做圆周运动的m1物体
的轨道半径用变量r表示,角度用变量a表示.
m1 r a m2
上一页
目录
返回
下一页
4
6.1 轨迹
例6.1:卫星绕地球转动时,m2等于地球的质量,m1等于卫 星的质量,r为卫星球心与地球球心间的距离。其运动轨迹由 下列方程组决定:
2
与例6.2相同,只是 改变了Alpha的值, 可以直接借用例6.2 的函数文件
上一页
目录
返回
下一页
19
6.2 单自由度系统
由初始条件建立执行文件(execute_63.m)
程序如下 zeta=0.2;Alpha=[0.00,-0.25,-0.25];
x0=[-2.00,-2.00,-2.00];v0=[2.00,2.00,2.31]; tspan=linspace(0.0,30.0,401); lintyp=char('-k','--k','-.k'); options=odeset('RelTol',1e-8,'AbsTol',[1e-8 1e-8]); d=char('Linear:x_0=-2 v_0=2 \alpha=0',... 'Nonlinear:x_0=-2 v_0=2 \alpha=-0.25',... 'Nonlinear:x_0=-2 v_0=2.31 \alpha=-0.25');
15

matlab求解动力学微分方程

matlab求解动力学微分方程

实用文档matlab求解动力学微分方程如今随着科学技术的持续发展和进步,动力学微分方程的求解成为了科研工作和工程应用中的一项基本任务。

作为一种广泛应用的计算工具,MATLAB可以通过其强大的数值计算和仿真功能来解决这一问题。

本文将深入探讨MATLAB在求解动力学微分方程方面的应用,包括其基本原理、解决方法以及一些实例分析,旨在帮助读者更全面地理解这一主题。

1. 动力学微分方程简介动力学微分方程是描述物质或系统中的运动过程的数学模型。

它们通常通过描述物体的运动、变化或响应来研究和分析不同领域的问题,例如物理、化学、生物学和工程。

2. MATLAB在求解动力学微分方程中的基本原理MATLAB提供了许多用于求解微分方程的函数和工具箱。

其中最常用且强大的函数是ode45,它基于龙格-库塔方法实现了自适应步长控制和高阶插值技术,可以有效地求解一般形式的动力学微分方程。

实用文档3. 使用MATLAB求解动力学微分方程的实例为了更好地理解MATLAB在求解动力学微分方程中的应用,我们将通过一些具体的例子来演示其使用方法。

我们可以考虑一个简单的弹簧振动方程,其中有一个质点通过弹簧受到外力作用。

通过建立该系统的微分方程模型,并利用MATLAB进行求解,我们可以得到质点的运动轨迹和其他相关信息。

4. 对MATLAB求解动力学微分方程的个人观点和理解作为一个计算工具,MATLAB无疑为求解动力学微分方程提供了便利和高效的方式。

其强大的数值计算和仿真功能能够帮助研究人员和工程师更好地理解和分析系统的运动行为。

然而,我们也应该注意,对于一些复杂的非线性动力学问题,可能需要更高级的数值方法和算法才能得到准确的解。

MATLAB作为一种常用的计算工具,在求解动力学微分方程方面具有广泛的应用。

通过掌握其基本原理和使用方法,我们可以有效地解决各类动力学问题,并更好地理解系统的运动行为。

当然,对于更复杂的问题,我们也应该不断地学习和探索更高级的数值方法,以求得更准确的解。

第6章 Matlab平面连杆机构的动力学分析

第6章  Matlab平面连杆机构的动力学分析

§6-1 曲柄的动力学仿真模块
由运动学知识可推得:
Re i Re A rcii cos i 2 rcii2 cos i s Im i Im A rcii sin i 2 rcii2 sin i s
§6-1 曲柄的动力学仿真模块
1.曲柄的动力学矩阵表达式 曲柄AB复向量的模 ri 为常数、幅角 i 为变量。 质心到转动副A的距离为 rci ,质量为 mi ,绕质心的转动惯量为 Ji , 作用于质心上的外力为 Fxi 和 Fyi 、 外力矩为M i ,曲柄与机架联接, 转动副A的约束反力为 RxA 和 RyA , 驱动力矩为 M 1 。
由理论力学可得:
RxA RxB Fxi mi Re i s
RyA RyB Fyi mi g mi Im i s
M1 M i RxArci sin i RyArci cos i RxB ri rci sin i RyB ri rci cos i J ii
§6-2 RRR II级杆组的动力学仿真模块
2.RRR II级杆组动力学分析M函数
g=9.8; %重力加速度 ri=1; rj=07; %两杆的长度 rci=0.5;rcj=0.35; %质心到铰链B的距离 %质心到铰链D的距离 mi=3; mj=2.2; %两杆的质量 Ji=0.25;Jj=0.09;%两杆的转动惯量 ReddD=0;ImddD=0; Fxi=0;Fyi=0;Mi=O; %i杆的外力和外力矩 a=zeros(6); a(1,1)=1;a(1,3)=1; a(2,2)=1; a(2,4)=1; a(3,1)=rci*sin(x(1)); a(3,2)=-rci*cos(x(1)); a(3,3)=-(ri-rci)*sin(x(1)); a(3,4)=(ri-rci)*cos(x(1)); a(4,3)=-1; a(4,5)=1; a(5,4)=-1; a(5,6)=1; a(6,3)=(rj-rcj)*sin(x(2)); a(6,4)=-(rj-rcj)*cos(x(2)); a(6,5)=rcj*sin(x(2)); a(6,6)=-rcj*cos(x(2));

基于MATLAB_SIMULINK的牛头刨床导杆机构运动学及动力学分析

基于MATLAB_SIMULINK的牛头刨床导杆机构运动学及动力学分析

基于MATLAB/SIMULINK的牛头刨床导杆机构运动学及动力学分析Kinematics and dynamics analysis of guide-bar mechanism in shaping machine based on MATLAB/SIMULINK钱文婷1,徐承妍2,李滨城1QIAN Wen-ting1, XU Cheng-yan2, LI Bin-cheng1(1. 江苏科技大学机械工程学院,镇江 212003;2. 华东师范大学软件学院,上海 200062)摘 要:本文运用MATLAB/SIMULINK 软件对牛头刨床导杆机构进行了运动学和动力学分析,并通过仿真直观地揭示了机构的运动规律和受力情况,为对其深入研究提供了基础。

另外,在应用SIMULINK进行运动学仿真过程中采用了微分的方法,从而使数学建模过程更为简便。

关键词:牛头刨床导杆机构;运动学分析;动力学分析; MATLAB/SIMULINK中图分类号:TP321 文献标识码:A 文章编号:1009-0134(2011)2(上)-0104-03 Doi: 10.3969/j.issn.1009-0134.2011.2(上).330 引言随着计算机技术的发展,在机构综合中计算机辅助分析得到了迅猛发展,特别为其建模仿真提供了极大的方便[1],为后续综合出机械最优化设计参数提供了可能。

牛头刨床是一种常见的金属切削机床,其导杆机构作为牛头刨床的主要执行机构,实现将回转运动转变为直线往复运动的重要功能,如图1所示。

本文运用矢量法和矩阵法建立了牛头刨床导杆机构分析的数学模型,在对其进行运动学和动力学分析的基础上,运用MATLAB/SIMULINK进行了运动学和动力学仿真,从而获得了机构的运动曲线图及运动副反力曲线图,为对其进一步深入研究提供了基础。

1 运动学分析1.1 数学模型的建立建立如图1所示的直角坐标系,将各构件视为杆矢量。

车辆工程基于MATLAB的动力性仿真分析及优化设计程序

车辆工程基于MATLAB的动力性仿真分析及优化设计程序

n=linspace(600,4000,100);%均分计算指令,600最低转速,4000最高转速,均分为100等分r=0.367;i0=5.83;nt=0.85;G=3880*9.8;f=0.013;CDA=2.77;If=0.218;Iw1=1.798;Iw2=3.598;m=3880;L=3.2;a=1.947;hg=0.9;ig=[6.09,3.09,1.71,1.00];%输入参数ua1=0.377*r*n/i0/ig(1);ua2=0.377*r*n/i0/ig(2);ua3=0.377*r*n/i0/ig(3);ua4=0.377*r*n/i0/ig(4);%各转速各挡位下的速度Tq=-19.313+295.27*(n/1000)-165.44*(n/1000).^2+40.874*(n/1000).^3-3.8445*(n/1000).^4;%从600~4000rpm油拟合公式计算发动机转距Ft1=Tq*i0*ig(1)*nt/r;Ft2=Tq*i0*ig(2)*nt/r;Ft3=Tq*i0*ig(3)*nt/r;Ft4=Tq*i0*ig(4)*nt/r;%从600~4000rpm各挡位的驱动力Ff=G*f;ua=linspace(0,200,100);Fw=CDA*ua.*ua/21.15;%空气阻力plot(ua1,Ft1,ua2,Ft2,ua3,Ft3,ua4,Ft4,ua,Ff+Fw);%画出各挡位的Ua-Ft,及Ua-Ff+Ftxlabel('ua/ km/h');ylabel('F/N');%标注横纵轴title('汽车驱动力-行驶阻力平衡图');%标注图形题目gtext('Ft1'),gtext('Ft2'),gtext('Ft3'),gtext('Ft4'),gtext('Ff+Fw');%给每根线条添加符号legend('Ft1','Ft2','Ft3','Ft4','Ff+Fw');%标注图例umax=max(ua4);disp('汽车最高车速=');disp(umax);disp('km/h');imax=tan(asin(max((Ft1-(Ff+Fw))/G)));%最大爬坡度的公式disp('汽车最大爬坡度=');disp(imax);%输出最高车速,与最大爬坡度的结果n=600:1:4000;%600最低转速,4000最高转速,相邻数组间隔1 r=0.367;i0=5.83;eff=0.85;f=0.013;CdA=2.77;m=3880;g=9.8; %输入参数G=m*g;Ttq=-19.313+295.27*n/1000-165.44*(n/1000).^2+40.874*(n/1000).^3-3.8445*(n/1000).^4; %从600~4000rpm油拟合公式计算发动机转距for ig=[6.09,3.09,1.71,1.00]Ua=0.377*r*n/ig/i0; %各转速各挡位下的速度Pe=Ttq.*n/9550; %各转速下的功率plot(Ua,Pe);hold on; %使当前轴及图形保持而不被刷新,准备承受此后将绘制的图形,多图共存endUa=0:0.1:max(Ua);Pf=G*f*Ua/3600; %滚动阻力Pw=CdA*Ua.^3/76140; %空气阻力plot(Ua,(Pf+Pw)/eff);title('汽车的功率平衡图'),xlabel('Ua/(km/h)'),ylabel('P/kw');%画出汽车的功率平衡图gtext('Ft1'),gtext('Ft2'),gtext('Ft3'),gtext('Ft4'),gtext('(Pf+Pw)/nt');legend('Ⅰ','Ⅱ','Ⅲ','Ⅳ','Pf+Pw/nt');n=600:1:4000;%600最低转速,4000最高转速,相邻数组间隔r=0.367;i0=5.83;nt=0.85;f=0.013;CdA=2.77;m=3880;g=9.8; %输入参数G=m*g;Ttq=-19.313+295.27*n/1000-165.44*(n/1000).^2+40.874*(n/1000).^3-3.8445*(n/1000).^4; %从600~4000rpm油拟合公式计算发动机转距for ig=[6.09,3.09,1.71,1.00]Ua=0.377*r*n/ig/i0;Ft=Ttq*i0*ig*nt/r;Fw=CdA*Ua.^2/21.15;D=(Ft-Fw)/G %汽车动力因子公式plot(Ua,D); %画出汽车动力特性图hold on; %使当前轴及图形保持而不被刷新,准备承受此后将绘制的图形,多图共存endf=0.0076+0.000056*Ua%滚动阻力与速度之间的关系plot(Ua,f); %画出速度与滚动阻力图title('汽车动力特性图'),%给图加题目xlabel('Ua/(km/h)'),ylabel('D');gtext('Ⅰ'),gtext('Ⅱ'),gtext('Ⅲ'),gtext('Ⅳ'),gtext('f');legend('Ⅰ','Ⅱ','Ⅲ','Ⅳ','f');n=600:10:4000; %600最低转速,4000最高转速,相邻数组间隔10m=3880;g=9.8;nmin=600;nmax=4000;G=m*g;ig=[6,09 3.09 1.71 1.00];nT=0.85;r=0.367;f=0.013;CDA=2.77;i0=5.83;L=3.2;a=1.947;hg=0.9;If=0.218;Iw1=1.798;Iw2=3.598;%输入参数Tq=-19.313+295.27*(n/1000)-165.44*(n/1000).^2+40.874*(n/1000).^3-3.8445*(n/1000).^4;%从600~4000rpm油拟合公式计算发动机转距Ft1=Tq*ig(1)*i0*nT/r;Ft2=Tq*ig(2)*i0*nT/r;Ft3=Tq*ig(3)*i0*nT/r;Ft4=Tq*ig(4)*i0*nT/r; %各转速各挡位下的驱动力ua1=0.377*r*n/ig(1)/i0;ua2=0.377*r*n/ig(2)/i0;ua3=0.377*r*n/ig(3)/i0;ua4=0.377*r*n/ig(4)/i0; %各挡位各转速下的速度Fw1=CDA*ua1.^2/21.15;Fw2=CDA*ua2.^2/21.15;Fw3=CDA*ua3.^2/21.15;Fw4=CDA*ua4.^2/21.15; %不同速度下的空气阻力Ff=G*f;deta1=1+(Iw1+Iw2)/(m*r^2)+(If*ig(1)^2*i0^2*nT)/(m*r^2);deta2=1+(Iw1+Iw2)/(m*r^2)+(If*ig(2)^2*i0^2*nT)/(m*r^2);deta3=1+(Iw1+Iw2)/(m*r^2)+(If*ig(3)^2*i0^2*nT)/(m*r^2);deta4=1+(Iw1+Iw2)/(m*r^2)+(If*ig(4)^2*i0^2*nT)/(m*r^2); %不同挡位下的汽车旋转质量换算系数a1=(Ft1-Ff-Fw1)/(deta1*m);ad1=1./a1;a2=(Ft2-Ff-Fw2)/(deta2*m);ad2=1./a2;a3=(Ft3-Ff-Fw3)/(deta3*m);ad3=1./a3;a4=(Ft4-Ff-Fw4)/(deta4*m);ad4=1./a4; %各挡位下的加速度plot(ua1,ad1,ua2,ad2,ua3,ad3,ua4,ad4);title('汽车的加速度倒数曲线');xlabel('ua(km/h)'); ylabel('1/a〕'); %作汽车加速度倒数曲线gtext('1/a1'),gtext('1/a2'),gtext('1/a3'),gtext('1/a4');legend('1/a1','1/a2','1/a3','1/a4');n=600:10:4000;m=3880;g=9.8;nmin=600;nmax=4000;G=m*g;ig=[6.09 3.09 1.71 1.00];nT=0.85;r=0.367;f=0.013;CDA=2.77;i0=5.83;L=3.2;a=1.947;hg=0.9;If=0.218;Iw1=1.798;Iw2=3.598; %输入参数Tq=-19.313+295.27*(n/1000)-165.44*(n/1000).^2+40.874*(n/1000).^3-3.8445*(n/1000).^4; %从600~4000rpm油拟合公式计算发动机转距Ft1=Tq*ig(1)*i0*nT/r;Ft2=Tq*ig(2)*i0*nT/r;Ft3=Tq*ig(3)*i0*nT/r;Ft4=Tq*ig(4)*i0*nT/r;%各转速各挡位下的驱动力ua1=0.377*r*n/ig(1)/i0;ua2=0.377*r*n/ig(2)/i0;ua3=0.377*r*n/ig(3)/i0;ua4=0.377*r*n/ig(4)/i0;%各挡位各转速下的速度Fw1=CDA*ua1.^2/21.15;Fw2=CDA*ua2.^2/21.15;Fw3=CDA*ua3.^2/21.15;Fw4=CDA*ua4.^2/21.15;%不同速度下的空气阻力Ff=G*f;i1=asin((Ft1-Ff-Fw1)/G);i2=asin((Ft2-Ff-Fw2)/G);i3=asin((Ft3-Ff-Fw3)/G);i4=asin((Ft4-Ff-Fw4)/G);%不同档位下的坡度plot(ua1,i1,ua2,i2,ua3,i3,ua4,i4);title('汽车的爬坡度图');xlabel('ua/(km*h^-1)');ylabel('i/%');%作汽车的坡度图gtext('Ⅰ'),gtext('Ⅱ'),gtext('Ⅲ'),gtext('Ⅳ');m=3880;g=9.8;r=0.367;nt=0.85;f=0.013;CdA=2.77;i0=5.83;pg=7.1;%汽油的重度取7.1N/Lig=[6.09 3.09 1.71 1];n=600:1:4000;n0=[815 1207 1614 2012 2603 3006 3403 3804];B00=[1326.8 1354.7 1284.4 1122.9 1141.0 1051.2 1233.9 1129.7];B10=[-416.46 -303.98 -189.75 -121.59 -98.893 -73.714 -84.478 -45.291];B20=[72.379 36.657 14.524 7.0035 4.4763 2.8593 2.9788 0.71113];B30=[-5.8629 -2.0553 -0.51184 -0.18517 -0.091077 -0.05138 -0.047449 -0.00075215];B40=[0.17768 0.043072 0.0068164 0.0018555 0.00068906 0.00035032 0.00028230 -0.000038568]; %输入参数B0=spline(n0,B00,n);B1=spline(n0,B10,n);B2=spline(n0,B20,n);B3=spline(n0,B30,n);B4=spline(n0,B40,n);%使用三次样条插值,保证曲线的光滑连续ua3=0.377*r*n/ig(3)/i0;ua4=0.377*r*n/ig(4)/i0; %求出发动机转速围对应的3、4档车速Pe3=(m*g*f*ua3/3600+CdA*ua3.^3/76140)/0.85;Pe4=(m*g*f*ua4/3600+CdA*ua4.^3/76140)/0.85; %发动机功率for i=1:1:3401 %用拟合公式求出各个燃油消耗率b3(i)=B0(i)+B1(i)*Pe3(i)+B2(i)*Pe3(i).^2+B3(i)*Pe3(i).^3+B4(i)*Pe3(i).^4;b4(i)=B0(i)+B1(i)*Pe4(i)+B2(i)*Pe4(i).^2+B3(i)*Pe4(i).^3+B4(i)*Pe4(i).^4; %插值得出对应速度的燃油消耗率endQ3=Pe3.*b3./(1.02.*ua3.*pg);Q4=Pe4.*b4./(1.02.*ua4.*pg); %3.4挡等速百公里燃油消耗量plot(ua3,Q3,ua4,Q4);title('最高档与次高档等速百公里油耗曲线'); %画出最高档与次高档等速百公里油耗曲线xlabel('ua(km/h)'); ylabel('百公里油耗〔L/100km〕');gtext('3档'),gtext('4档');。

汽车系统动力学Matlab

汽车系统动力学Matlab

汽车系统动力学Matlab作业报告小组成员:一、组内任务分配二、Matlab程序与图形1、不同转向特性车辆在不同车速下的系统特征根m=1000;I=1500;a1=1。

15;b1=1。

35;Caf=53000;Car=53000;i=1;R=[];for uc=10:5:100;D=(I*(Caf+Car)+m*(a1^2*Caf+b1^2*Car))/(m*I*uc);S=(a1+b1)^2*Caf*Car/(m*I*uc^2)+(b1*Car—a1*Caf)/I;P=[1 D S];r=roots(P);R(i,1)=r(1,1);R(i,2)=r(2,1);i=i+1;endplot(real(R(:,1)),imag(R(:,1)),’bo’);holda2=1.25;b2=1。

25;t=1;S=[];for uc=10:5:100P=[m 0;0 I];Q=[(Caf+Car)/uc,m*uc+(a2*Caf—b2*Car)/uc;(a2*Caf—b2*Car)/uc,(a2^2*Caf+b2^2*Car)/uc];R=[Caf;a2*Caf];A=-P^(—1)*Q;d=eig(A);i=imag(d);r=real(d);S(t,1)=r(1);S(t,2)=i(1);t=t+1;endplot(S(:,1),S(:,2),’*')a3=1.35;b3=1。

15;for uc=10:5:100P=[m 0;0 I];Q=[(Caf+Car)/uc,m*uc+(a3*Caf—b3*Car)/uc;(a3*Caf-b3*Car)/uc,(a3^2*Caf+b3^2*Car)/uc];R=[Caf;a3*Caf];A=-P^(—1)*Q;d=eig(A);i=imag(d);r=real(d);S(t,1)=r(1);S(t,2)=i(1);t=t+1;endgrid onplot(S(:,1),S(:,2),’d’);axis([-14 2 0 3]);xlabel('实轴(Re)');ylabel('虚轴(Im)');text(-8,2.8,'不足转向');text(0,0.2,'过多转向');text(-3,0。

matlab求解动力学微分方程

matlab求解动力学微分方程

matlab求解动力学微分方程摘要:I.引言A.动力学微分方程简介B.Matlab在求解动力学微分方程中的应用II.Matlab基础知识回顾A.Matlab简介B.Matlab的基本操作C.Matlab的数值计算功能III.动力学微分方程的Matlab求解方法A.常微分方程的数值求解1.欧拉方法2.改进欧拉方法3.龙格-库塔方法B.动力学微分方程的符号求解1.符号运算的基本概念2.符号求解动力学微分方程IV.案例分析A.弹簧振子的动力学模型1.建立弹簧振子的微分方程2.使用Matlab求解弹簧振子的运动轨迹B.简谐振子的动力学模型1.建立简谐振子的微分方程2.使用Matlab求解简谐振子的运动轨迹V.结论A.Matlab在求解动力学微分方程中的优势B.未来发展方向正文:I.引言动力学微分方程是描述物体运动规律的数学模型,广泛应用于物理学、工程学等领域。

Matlab作为一种功能强大的数学软件,可以方便地求解动力学微分方程,为相关领域的研究提供有力支持。

本文将介绍Matlab求解动力学微分方程的基本方法,并通过案例分析展示其应用效果。

II.Matlab基础知识回顾Matlab是一种数值计算软件,具有丰富的函数库和强大的绘图功能。

在求解动力学微分方程之前,我们先简要回顾一下Matlab的基本知识。

A.Matlab简介Matlab(Matrix Laboratory)是一款由美国MathWorks公司开发的数学软件,主要用于数值计算、数据分析、可视化以及算法开发等。

Matlab具有丰富的函数库,用户可以通过编写简单的指令对数据进行处理和分析。

B.Matlab的基本操作Matlab的基本操作包括:变量赋值、矩阵运算、图形绘制等。

用户可以通过这些操作对数据进行处理和分析。

C.Matlab的数值计算功能Matlab提供了丰富的数值计算功能,包括求和、求积、求导、积分等。

此外,Matlab还支持符号计算,可以处理复杂的数学表达式。

matlab 拉格朗日 动力学方程

matlab 拉格朗日 动力学方程

matlab 拉格朗日动力学方程【指定主题】MATLAB拉格朗日动力学方程【引言】MATLAB是一款广泛应用于科学计算和数据分析的工具,它具备强大的数学计算能力,包括用于求解动力学问题的拉格朗日方程。

拉格朗日动力学方程是描述物体运动的基本方程之一,通过该方程可以推导出物体在给定条件下的运动轨迹和动力学行为。

本文将介绍MATLAB 中求解拉格朗日动力学方程的方法,并探讨其在物体运动研究中的应用。

【目录】1. 第一部分:拉格朗日动力学方程简介1.1 拉格朗日动力学方程的定义1.2 拉格朗日方程的推导方法1.3 MATLAB中的拉格朗日方程求解2. 第二部分:MATLAB实例演示2.1 简谐振子的运动模拟2.2 刚体的运动模拟2.3 多体系统的模拟3. 第三部分:拉格朗日动力学方程的应用3.1 机械臂的建模与控制3.2 汽车悬挂系统的分析3.3 飞行器的动力学建模与仿真4. 总结与展望4.1 文章总结4.2 个人观点与经验分享4.3 对未来发展的展望1. 第一部分:拉格朗日动力学方程简介1.1 拉格朗日动力学方程的定义拉格朗日动力学方程是描述物体运动的基本方程之一。

它基于钦定了某种能量函数——拉格朗日量,并通过对该能量函数进行变分运算得到。

拉格朗日动力学方程可以从哈密顿原理或变分原理推导出来,它可以用于建立物体的动力学模型,研究物体的运动规律和相互作用。

1.2 拉格朗日方程的推导方法拉格朗日方程的推导可以通过哈密顿原理或变分原理进行。

哈密顿原理是一种最小作用量原理,它假设物体遵循最小作用量路径进行运动。

通过对作用量进行变分运算,并使用欧拉-拉格朗日方程,就可以得到拉格朗日方程。

变分原理是通过对拉格朗日量进行变分运算来得到拉格朗日方程。

1.3 MATLAB中的拉格朗日方程求解在MATLAB中,可以利用符号计算功能和数值求解方法来求解拉格朗日方程。

符号计算功能可以将拉格朗日方程中的变量和函数表示为符号表达式,从而进行符号计算和求解。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

1
1.微分方程的定义
对于duffing 方程03
2
=++x x x ω
,先将方程写作⎩⎨⎧
--==3
1122
21x x x x x ω function dy=duffing(t,x) omega=1;%定义参数 f1=x(2);
f2=-omega^2*x(1)-x(1)^3; dy=[f1;f2];
2.微分方程的求解
function solve (tstop) tstop=500;%定义时间长度 y0=[0.01;0];%定义初始条件
[t,y]=ode45('duffing',tstop,y0,[]);
function solve (tstop) step=0.01;%定义步长
y0=rand(1,2);%随机初始条件
tspan=[0:step:500];%定义时间范围 [t,y]=ode45('duffing',tspan,y0);
3.时间历程的绘制
时间历程横轴为t ,纵轴为y ,绘制时只取稳态部分。

plot(t,y(:,1));%绘制y 的时间历程 xlabel('t')%横轴为t ylabel('y')%纵轴为y grid;%显示网格线
2
axis([460 500 -Inf Inf])%图形显示范围设置
4.相图的绘制
相图的横轴为y ,纵轴为dy/dt ,绘制时也只取稳态部分。

红色部分表示只取最后1000个点。

plot(y(end-1000:end ,1),y(end-1000:end ,2));%绘制y 的时间历程
xlabel('y')%横轴为y
ylabel('dy/dt')%纵轴为dy/dt grid;%显示网格线
5.Poincare 映射的绘制
对于不同的系统,Poincare 截面的选取方法也不同
对于自治系统一般每过其对应线性系统的固有周期,截取一次 对于非自治系统,一般每过其激励的周期,截取一次
例程:duffing 方程03
2=++x x x ω
的poincare 映射 function poincare(tstop)
global omega; omega=1;
T=2*pi/omega;%线性系统的周期或激励的周期 step=T/100;%定义步长为T/100 y0=[0.01;0];%初始条件
tspan=[0:step:100*T];%定义时间范围 [t,y]=ode45('duffing',tspan,y0);
for i=5000:100:10000%稳态过程每个周期取一个点 plot(y(i,1),y(i,2),'b.'); hold on;% 保留上一次的图形 end
xlabel('y');ylabel('dy/dt');
3
Poincare 映射也可以通过取极值点得到 function poincare(tstop) y0=[0.01;0];
tspan=[0:0.01:500];
[t,y]=ode45('duffing',tspan,y0); count=find(t>100);%截取稳态过程 y=y(count,:);
n=length(y(:,1));%计算点的总数 for i=2:n-1
if y(i-1,1)+eps<y(i,1) && y(i,1)>y(i+1,1)+eps % 简单的取出局部最大值
plot(y(i,1),y(i,2),'.'); hold on end end
xlabel('y');ylabel('dy/dt');
6.频谱
yy=fft(y(end-1000:end,1)); N=length(yy); power=abs(yy);
freq=(1:N-1)*1/step/N;
plot(freq(1:N/2),power(1:N/2)); xlabel('f(y)') ylabel('y')
7.算例
duffing 方程03
=++x x x
的时间历程,相图,频谱和poincare 映射。

function dy=duffing(t,x)
omega=1;%定义参数
f1=x(2);
f2=-omega^2*x(1)-x(1)^3;
dy=[f1;f2]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function duffsim(tstop)
step=0.01
y0=[0.1;0];
tspan=[0:step:500];
[t,y]=ode45('duffing',tspan,y0); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,2,1)
plot(t,y(:,1));%绘制y的时间历程
xlabel('t')%横轴为t
ylabel('y')%纵轴为y
grid;%显示网格线
axis([460 500 -Inf Inf])%显示范围设置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,2,2)
plot(y(end-1000:end,1),y(end-1000:end,2));%绘制y的时间历程
xlabel('y')%横轴为y
ylabel('dy/dt')%纵轴为dy/dt
grid;%显示网格线%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,2,3)
yy=fft(y(end-1000:end,1));
4
N=length(yy);
power=abs(yy);
freq=(1:N-1)*1/step/N;
plot(freq(1:N/2),power(1:N/2));
xlabel('f(y)')
ylabel('y') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,2,4)
count=find(t>100);%截取稳态过程
y=y(count,:);
n=length(y(:,1));%计算点的总数
for i=2:n-1
if y(i-1,1)+eps<y(i,1) && y(i,1)>y(i+1,1)+eps % 简单的取出局部最大值
plot(y(i,1),y(i,2),'.');hold on;
end
end
xlabel('y');ylabel('dy/dt');
5
8.分岔图的绘制
随F变化的分岔图。

+
-
+
t
3.03=
cos
x
x
x
x2.1
F
function dy=duffing(t,x)
global c;
omega=1;%定义参数
f1=x(2);
f2=omega^2*x(1)-x(1)^3-0.3*x(2)+c*cos(1.2*t); dy=[f1;f2]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
global c; %定义全局变量
range=[0.1:0.002:0.9];%定义参数变化范围
k=0;
6
YY=[];%定义空数组
for c=range
y0=[0.1;0];%初始条件
k=k+1;
tspan=[0:0.01:400];
[t,Y]=ode45('duffing',tspan,y0);
count=find(t>200);
Y=Y(count,:);
j=1;
n=length(Y(:,1));
for i=2:n-1
if Y(i-1,1)+eps<Y(i,1) && Y(i,1)>Y(i+1,1)+eps % 简单的取出局部最大值。

YY(k,j)=Y(i,1);
j=j+1;
end
end
if j>1
plot(c,YY(k,[1:j-1]),'k.','markersize',3);
end
hold on;
index(k)=j-1;
end
xlabel('c');
ylabel('y');
7
随F变化的分岔图
F=0.20
8
F=0.27
F=0.275
F=0.2875
9
F=0.32
F=0.36
F=0.4
10
F=0.652 F=0.8
11。

相关文档
最新文档