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 中求解拉格朗日动力学方程的方法，并探讨其在物体运动研究中的应用。
【目录】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中，可以利用符号计算功能和数值求解方法来求解拉格朗日方程。
符号计算功能可以将拉格朗日方程中的变量和函数表示为符号表达式，从而进行符号计算和求解。
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仿真分析工具箱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程序-概述说明以及解释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作为一种强大的数学计算软件，提供了丰富的工具和函数用于求解微分方程。
其拥有的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在求解动力学微分方程方面具有非常大的优势。
其丰富的函数库和高效的求解器为研究人员提供了便利，使他们能够更深入地探究动力学系统的行为。
机械原理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应用之动力学与振动
一.运动微分方程 当
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
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
··
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;%显示网格线
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');
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));
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');
8.分岔图的绘制
随F变化的分岔图。
3.03=
+
-
+
x2.1
F
t
x
x
x
cos
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;
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');
··
F=0.20
··
·
·
F=0.27
F=0.275
F=0.2875
·
·
F=0.32
F=0.36
F=0.4
··
F=0.652
F=0.8。