ode45求解多自由度动力学方程实例
ode45求解多自由度动力学方程实例
在matlab命令窗口里输入一下命令：
y0=[1 1 1 1 1 1];
tspan=[0 30];
[t,y]=ode45(@func,tspan,y0);
figure(1)
plot(t,y(:,1),t,y(:,3),t,y(:,5));
legend('x1','x2','x3');
xlabel('时间(s)','FontSize',10);
ylabel('振动位移曲线','FontSize',10);
figure(2)
plot(t,y(:,2),t,y(:,4),t,y(:,6));
legend('v1','v2','v3');
xlabel('时间(s)','FontSize',10);
ylabel(‘振动速度曲线’,’FontSize’,10);
运行结果：
Ode45函数调用形式如下：
[T,Y]=ode45(odefun,tspan,y0)un
用于存放待求解的方程的m文件名，方程必须用y’=f(t,y)的形式存放
tspan
指定自变量范围的向量，通常用[t0，tf]指定
y0
函数的边界条件，即y0=y(t0)，对于方程组，y0也可以是向量
例：若一三自由度多体动力学系统方程如下：
初始条件：
由于方程必须用y’=f(t,y)的形式存放，因此需要对方程组进行降阶处理。
令
则方程组可化为：
因此建立M函数文件来定义此方程组如下：
function dy=func(t,y)
弹簧-质量-阻尼实验指导书汇总
质量-弹簧-阻尼系统实验教学指导书北京理工大学机械与车辆学院2016.3实验一：单自由度系统数学建模及仿真 1 实验目的（1）熟悉单自由度质量-弹簧-阻尼系统并进行数学建模； （2）了解MATLAB 软件编程，学习编写系统的仿真代码； （3）进行单自由度系统的仿真动态响应分析。
2 实验原理单自由度质量-弹簧-阻尼系统，如上图所示。
由一个质量为m 的滑块、一个刚度系数为k 的弹簧和一个阻尼系数为c 的阻尼器组成。
系统输入：作用在滑块上的力f (t )。
系统输出：滑块的位移x (t )。
建立力学平衡方程：m x c x kx f ∙∙∙++=变化为二阶系统标准形式：22f x x x mζωω∙∙∙++=其中：ω是固有频率，ζ是阻尼比。
ω=2c m ζω== 2.1 欠阻尼(ζ<1)情况下，输入f (t )和非零初始状态的响应：()()sin()))]t t x t t d e ζωττζωττ+∞--=∙-=-+-+⎰2.2 欠阻尼(ζ<1)情况下，输入f(t)=f0*cos(ω0*t) 和非零初始状态的的响应：02230022222002222222()cos(arctan())2f[(0)]cos()[()(2)]sin(ttx t tx ekeζωζωζωωωωωζωωωωζωω-∙-=--++-++)输出振幅和输入振幅的比值：A=3 动力学仿真根据数学模型，使用龙格库塔方法ODE45求解，任意输入下响应结果。
仿真代码见附件4 实验4.1 固有频率和阻尼实验（1）将实验台设置为单自由度质量-弹簧-阻尼系统。
（2）关闭电控箱开关。
点击setup菜单，选择Control Algorithm，设置选择Continuous Time Control，Ts=0.0042，然后OK。
（3）点击Command菜单，选择Trajectory，选取step,进入set-up,选取Open Loop Step 设置(0)counts, dwell time=3000ms,(1)rep, 然后OK。
微分方程的数值解法matlab(四阶龙格—库塔法)
解析解： x x x1 3 2(((ttt))) 0 .0 8 1 1 2 P k 8 0siw n t) (2 .6 3 0 3 3 P k 0siw n t) (0 .2 12 2 2 P k 0siw n t)(
第一个质量的位移响应时程
Y (t)A(Y t)P(t)
(2)
Y (t)A(Y t)P(t)
3. Matlab 程序(主程序：ZCX)
t0;Y0;h;N;P0,w; %输入初始值、步长、迭代次数、初始激励力；
for i = 1 : N
t1 = t0 + h
P=[P0*sin(w*t0);0.0;0.0]
%输入t0时刻的外部激励力
Van der Pol方程
% 子程序 （程序名： dYdt.m ） function Ydot = dYdt (t, Y) Ydot=[Y(2);-Y(2)*(Y(1)^2-1)-Y(1)];
或写为
function Ydot = dYdt (t, Y) Ydot=zeros(size(Y)); Ydot(1)=Y(2); Ydot(2)=-Y(2)*(Y(1).^2-1)-Y(1)];
Solver解算指令的使用格式
说明：
t0：初始时刻；tN：终点时刻 Y0：初值； tol：计算精度
[t, Y]=solver (‘ODE函数文件名’, t0, tN, Y0, tol);
ode45
输出宗量形式
y1 (t0 )
Y
y1
(t1
)
y
1
(t
2
)
y2 (t0 )
y
2
(
t1
)
y
2
(
t
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在求解动力学微分方程方面具有非常大的优势。
其丰富的函数库和高效的求解器为研究人员提供了便利，使他们能够更深入地探究动力学系统的行为。
ode45函数用法
ode45函数用法ode45函数是MATLAB中的一个常用函数，用于求解常微分方程初值问题。
它实质上是利用了4阶和5阶Runge-Kutta方法的一种变种形式，可以对一阶和高阶常微分方程进行数值求解。
ode45函数的主要语法格式为：[t, y] = ode45(@fun, tspan, y0)其中，@fun是一个函数句柄，用于定义微分方程dy/dt=f(t,y)的右侧函数f(t,y)。
tspan是求解的时间范围，通常为一个二元向量[tstart, tend]。
y0是微分方程的初始条件。
当调用ode45函数时，它会返回两个向量t和y，其中t是时间向量，y是对应时间点的解向量。
可以通过绘图或其他后续处理分析这些解。
ode45函数的使用方法具有一定的灵活性。
首先，可以通过在函数f(t,y)中定义额外的参数来解决一些特定问题。
例如，可以将常数或其他变量作为输入参数传递给f函数。
其次，可以在tspan中指定更复杂的时间范围。
除了指定[tstart, tend]之外，还可以指定等间隔的时间点或自定义时间点。
这样可以对非等间隔的时间步长进行模拟。
另外，ode45函数还可以通过指定其他选项来进行更精确的求解或控制数值计算。
例如，可以通过指定'AbsTol'和'RelTol'选项来调整计算的绝对误差和相对误差。
还可以通过指定'Events'选项来检测事件，并在事件发生时中断求解。
除了以上基本用法外，ode45函数还可以与其他函数和工具箱进行配合使用。
例如，可以使用odeplot函数在求解过程中实时绘制解曲线。
还可以使用odephas2函数对微分方程的相平面进行绘制。
总之，ode45函数是MATLAB中一个强大的求解常微分方程的工具，它不仅具有灵活的使用方法，而且可以与其他函数和工具箱进行配合使用。
熟练掌握ode45函数的用法对于数学建模和科学计算等应用非常有价值。
多自由度系统动力学方程
01 22 2 n 2
i ：第 i 阶固有频率
1 ：基频
仅取决于系统本身的刚度、质量等物理参数
多自由度系统振动 / 多自由度系统的自由振动
采用位移方程求解固有频率
位移方程： FM X XF(P t) XRn 自由振动的位移方程： FM X X0
FK1柔度矩阵
主振动： Xφ si nt () φ [1 2 n]T
作用力方程： M X KX P(t) 固有振动方程： M X KX 0
XRn M、 KRnn
P(t)Rn
在考虑系统的固有振动时，首先考察系统的同步振动，即系
统在各个坐标上除了运动幅值不相同外，随时间变化的规律
都相同的运动。
假设系统的运动为： Xφf(t)
XRn φ Rn
常数列向量 运动规律的时间函数 f (t)R1
FD (T T )1 FC
将 X D TXC 代入D点的方程，并左乘 T T ： T T MDTXC T T K DTXC T T FD FC
得： T T MDT MC
T T KDT KC
验证：
m MD me
me
IC
me2
m 0
MC
0
I
C
1 0 m me 1 e m 0
k212m21
k222m22
k2n 2m2n
0
kn12mn1 kn2 2mn2 knn2mnn
多自由度系统振动 / 多自由度系统的自由振动
k112m11 k122m12 k1n 2m1n
k212m21
k222m22
k2n 2m2n
0
kn12mn1 kn2 2mn2 knn2mnn
matlab利用ode45求解二元二阶微分方程
题目：探究matlab利用ode45求解二元二阶微分方程的方法与应用在数学和工程领域，微分方程是一类重要的数学工具，它可以描述自然界中众多的现象和规律。
而求解微分方程的问题一直是科学家和工程师们所关注的重要问题之一。
在计算机辅助数学建模领域，matlab作为一种强大的数值计算工具，可以通过内置的函数ode45来求解常微分方程初值问题。
本文将探讨matlab利用ode45求解二元二阶微分方程的方法与应用。
一、二元二阶微分方程的基本概念二元二阶微分方程是指含有两个自变量、二阶导数和一阶导数的微分方程。
一般形式如下：\[ F(x, y, \frac{dy}{dx}, \frac{d^2y}{dx^2}) = 0 \]其中x为自变量，y为因变量，\(\frac{dy}{dx}\)为y关于x的一阶导数，\(\frac{d^2y}{dx^2}\)为y关于x的二阶导数。
二、matlab中ode45函数的基本原理在matlab中，ode45是求解常微分方程初值问题的函数，它使用了一种自适应步长的Runge-Kutta方法来求解微分方程。
ode45可以求解一阶或高阶的常微分方程组，是matlab中最常用的求解微分方程的函数之一。
对于二元二阶微分方程，可以通过一些简单的变换和处理，转化为一组一阶微分方程的形式，然后利用ode45进行求解。
三、matlab利用ode45求解二元二阶微分方程的具体步骤1. 将二元二阶微分方程转化为一组一阶微分方程。
对于形如\(\frac{d^2y}{dx^2} = f(x, y, \frac{dy}{dx})\)的二阶微分方程，可以引入新的变量z = \(\frac{dy}{dx}\)，转化为一组一阶微分方程：\[\frac{dy}{dx} = z\]\[\frac{dz}{dx} = f(x, y, z)\]2. 编写matlab脚本文件。
在matlab中，编写脚本文件来定义微分方程的函数形式，并调用ode45函数来求解微分方程。
ode45的原理
ode45的原理引言在数值计算中，求解常微分方程是一个重要的问题。
常微分方程是描述自然界中很多现象的数学模型，例如物理学中的运动学、动力学问题，生物学中的生物传播模型等等。
由于常微分方程往往难以求得解析解，因此需要借助数值方法进行求解。
ode45是一种常用的求解常微分方程的数值方法，本文将对ode45的原理进行全面、详细、完整且深入地探讨。
什么是ode45ode45是MATLAB中的一个函数，用于求解常微分方程的初值问题。
它采用了一种称为Runge-Kutta-Fehlberg（RKF）的方法。
RKF方法是一种经典的自适应步长的数值积分方法，它能够在保证一定的精度的前提下自动调整步长，从而提高计算效率。
RKF方法的基本原理RKF方法是一种多步骤的迭代方法，其基本思想是将步长h分为两个部分：主步长h和辅助步长h/2。
在每个主步长上，RKF方法使用两个不同的步长来计算解的近似值，然后通过比较两个近似值的差异来估计误差，并根据误差的大小来调整步长。
具体的步骤如下：1.使用主步长h计算解的近似值y1和y2。
其中y1是通过一个四阶的Runge-Kutta方法计算得到的近似值，y2是通过一个五阶的Runge-Kutta方法计算得到的近似值。
2.计算y1和y2之间的差异delta。
如果差异太大，说明步长h可能过大，需要减小步长；如果差异较小，说明步长h可能适合，可以增大步长。
3.根据差异delta来调整步长h。
如果delta小于某个阈值，可以增大步长；如果delta大于某个阈值，需要减小步长。
步长的调整可以使用一些启发式的方法，例如根据delta的大小来乘以一个系数来调整步长。
4.重复步骤1~3，直到达到指定的终止条件，例如达到指定的时间点或者误差限度。
RKF方法的优点是可以在保证一定精度的前提下自动调整步长，从而提高计算效率。
然而，RKF方法也有一些缺点，例如计算量较大，不适合求解高阶的常微分方程等。
ode45函数的使用方法在MATLAB中，可以使用ode45函数来求解常微分方程的初值问题。
实验二MATLAB数值计算常微分方程的求解
实验二MATLAB数值计算常微分方程的求解引言:微分方程是描述物理问题和工程问题的重要工具，也是数学和工程学科中的重要课题。
解析方法可以求得一些简单微分方程的解析解，但对于复杂问题，常常无法得到解析解。
此时，数值求解方法成为一种有效的工具。
MATLAB是一个功能强大的数值计算软件，对于求解常微分方程组也提供了多种数值方法。
本实验将介绍MATLAB中求解常微分方程(组)的方法，并通过实例来演示这些方法的应用。
一、MATLAB的常微分方程(组)求解函数MATLAB提供了多种函数用于求解常微分方程(组)，其中最常用的函数是ode45、ode23和ode15s。
这些函数采用不同的数值方法，精度和效率也不同。
下面分别对这些函数进行简单介绍：1. ode45函数:ode45函数采用了一种自适应的步长控制算法，可以在一个时间段内自动调整步长。
这个函数的语法是：[t,y] = ode45(odefun,tspan,y0,options)其中，odefun是一个函数句柄，表示待解的常微分方程组，tspan是求解区间，y0是初始条件，options是一个结构体，包含其他的参数和选项。
2. ode23函数:ode23函数也采用了自适应的步长控制算法，但相对于ode45，它是一个简化版本，只允许使用比较简单的问题。
这个函数的语法与ode45相似：[t,y] = ode23(odefun,tspan,y0,options)3. ode15s函数:ode15s函数采用了一种隐式数值方法，适用于比较刚性（stiff）的常微分方程。
这个函数的语法与前两个函数也相似：[t,y] = ode15s(odefun,tspan,y0,options)以上是MATLAB提供的三种解常微分方程(组)的函数，根据问题的特点和求解的需求选择相应的函数。
二、实例演示在本实验中，我们将通过一个实例来演示如何使用MATLAB解常微分方程组。
假设有一个简单的线性常微分方程组：dy1/dt = -2y1 + y2dy2/dt = y1 - 2y2给定初始条件为y1(0)=1，y2(0)=0，求在t=[0,5]时，y1和y2的解。
ode45解微分方程组
ode45解微分方程组要解决微分方程组，我们可以使用MATLAB中的ODE45（常微分方程数值求解器）函数。
ode45可以求解形如dy/dt = f(t, y)的一阶常微分方程，其中t是自变量，y是因变量，f(t, y)是描述y在t上的变化率的函数。
对于多个变量的微分方程组，可以将其处理为向量形式并使用ode45以下是一个例子来说明如何使用ode45函数来解微分方程组。
假设我们要求解以下微分方程组：dy1/dt = y2dy2/dt = -2*y1 + 3*y2其中，y1和y2都是关于t的函数。
首先，我们需要将这个微分方程组转化为向量形式。
定义一个向量Y，其中Y(1)=y1，Y(2)=y2，我们可以将系统改写为：dY/dt = [Y(2); -2*Y(1) + 3*Y(2)]然后，在MATLAB中，我们可以定义一个函数来表示这个系统。
我们将其命名为"system"，输入参数t和Y，输出参数是dY。
function dY = system(t, Y)dY(1,1)=Y(2);dY(2,1)=-2*Y(1)+3*Y(2);end接下来，我们可以使用ode45函数来求解该微分方程组。
首先，我们需要指定初始条件。
令t0为初始时间，Y0为初始向量，我们可以从初始条件中获得这些值。
在这个例子中，假设我们要从t=0开始，并且y1(0)=1，y2(0)=0。
因此，t0=0，Y0=[1;0]。
然后，我们可以调用ode45，将解密度'R'设置为1e-3（即0.001），并指定时间间隔[0 10]：tspan = [0 10];Y0=[1;0];R=1e-3;最后，我们得到了在时间间隔[010]上的数值解，通过't'和'Y'来表示。
t是时间向量，Y是与时间相对应的解向量。
我们可以绘制结果以获得更好的理解。
假设我们想绘制y1和y2的图像。
y1=Y(:,1);y2=Y(:,2);figure;plot(t, y1, 'r', t, y2, 'b')xlabel('t')ylabel('y')legend('y1', 'y2')title('Solution of the differential equations')以上就是求解微分方程组的一个例子。
07-多自由度系统的运动方程
多自由度系统的运动微分方程牛顿第二定律矢量建模方法 影响系数法•刚度影响系数法•柔度影响系数法Lagrange方程方法•约束、自由度与广义坐标•Lagrange方程建模方法多自由度系统的运动微分方程牛顿第二定律建模多自由度系统的运动微分方程--牛顿第二定律建模以1m 为研究对象，有()())(1122111221111t F x x c x c x x k x k xm +−+−−+−=&&&&& （1） 以2m 为研究对象，有()())(2232321221212t F x c x k x x c x x k xm +−−−+−=&&&&& （2）将方程(1)、（2）整理可得多自由度系统的运动微分方程--牛顿第二定律建模()[]()[])(1221212212111t F x k x k k x c x c c xm =−++−++&&&& （3） ()[]()[])(2232122321212t F x k k x k x c c x c xm =++−+++−+&&&& （4） 将方程(3)、(4)写成矩阵形式⎭⎬⎫⎩⎨⎧=⎭⎬⎫⎩⎨⎧⎥⎦⎤⎢⎣⎡+−−++⎭⎬⎫⎩⎨⎧⎥⎦⎤⎢⎣⎡+−−++⎭⎬⎫⎩⎨⎧⎥⎦⎤⎢⎣⎡)()(002121322221213222212121t F t F x x k k k k k k x xc c c c c c x x m m &&&&&&即)(t F kx x C xM =++&&& 多自由度系统的运动微分方程--牛顿第二定律建模这种用矩阵写出的运动微分方程与单自由度系统的运动微分方程非常相似。
象例题中在各个离散质量上建立的坐标系为描述系统的物理坐标系，在此坐标下的系统质量矩阵、阻尼矩阵和刚度矩阵为系统的物理参数。
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 映射。
ode45 matlab 例子
ode45 matlab 例子ode45 MATLAB 例子1. 简介ode45 是 MATLAB 中一个常用的求解常微分方程（ODE）的函数。
它使用了一个基于 Runge-Kutta 方法的算法来解决ODE问题。
2. 线性ODE的例子下面我们通过一个线性ODE的例子来说明 ode45 的基本用法。
线性ODE的定义：假设有一个一阶线性ODE，形式如下：dy/dt = a*y + b其中 a 和 b 是常数，y 是未知函数。
例子：考虑一个简单的线性ODE：dy/dt = 2*t*y + t^2为了使用 ode45 解决这个问题，我们需要定义一个 MATLAB 函数来表示这个方程：function dydt = myODE(t, y)dydt = 2*t*y + t^2;end然后，我们可以使用 ode45 来求解该问题：[t, y] = ode45(@myODE, [0 1], 1);上述代码中： - @myODE表示将myODE函数作为变量传递给ode45 函数。
- [0 1]表示求解的时间区间是 [0, 1]。
- 1表示初始条件。
3. 非线性ODE的例子除了线性ODE，ode45 也可以求解非线性ODE。
下面我们通过一个非线性ODE的例子来说明。
非线性ODE的定义：假设有一个二阶非线性ODE，形式如下：d^2y/dt^2 + a*dy/dt + b*y + c*y^2 = 0其中 a、b 和 c 是常数，y 是未知函数。
例子：考虑一个简单的非线性ODE：d^2y/dt^2 + 2*dy/dt + 3*y + 4*y^2 = 0为了使用 ode45 解决这个问题，我们需要将这个二阶ODE转化为两个一阶ODE。
令 z = dy/dt，我们可以得到以下一阶ODE系统：dz/dt = -2*z - 3*y - 4*y^2dy/dt = zfunction dydt = myODE(t, y)dydt = zeros(2,1);dydt(1) = -2*y(2) - 3*y(1) - 4*y(1)^2;dydt(2) = y(1);end然后，我们可以使用 ode45 来求解该问题：[t, y] = ode45(@myODE, [0 1], [1 0]);上述代码中，初始条件是[1 0]，因为我们需要同时指定 y 和z 的初始值。
第十三讲—多自由度系统的运动方程
−
ka2θ2
+
⎛ ⎜⎝
ka 2
+
1 2
mgl
⎞⎟⎠θ3
=
0
⎡ ⎢ ⎢
1 3
ml
2
⎢ ⎢
0
⎢
0 1 ml2 3
⎢ ⎢⎣
0
0
质量矩阵
0
⎤ ⎥ ⎥
⎡⎢ka2 ⎢
+
1 2
mgl
0
⎥⎢ ⎥⎢
−ka2
1
ml
2
⎥ ⎥
⎢ ⎢
0
3 ⎥⎦ ⎢⎣
−ka2 2ka2 + 1 mgl
2 −ka2
刚度矩阵
0
⎤ ⎥
⎥
−ka2
⎥ ⎥
θ12
+
θ
2 2
+ θ32
O1
O2
O3
a
k
k
势能
θ1
θ2
θ3
U
=
1 2
k
( aθ1
−
aθ2
)2
+
1 2
k
( aθ 2
−
aθ3 )2
+
mg
l 2
(1 −
cos θ1 )
+
mg
l 2
(1 −
cosθ2
)
+
mg
l 2
(1 −
cos θ1 )
( ) 1
2
k
( aθ1
−
aθ2
)2
+
1 2
k
( aθ 2
−
aθ3
2i
j
mij
ode45 matlab 用法
在本文中，我将为您介细介绍MATLAB中的ode45函数的用法。
ode45函数是MATLAB中常用的求解常微分方程(ODE)的函数，可以对一阶或高阶ODE进行求解，并且能够处理刚体动力学、化学反应动力学、电路等各种领域的问题。
我将从ode45函数的基本用法、参数设定、示例应用以及个人理解这几个方面展开讨论。
1. ode45函数的基本用法ode45函数是MATLAB中常用的求解ODE的函数。
其基本用法为：[t, y] = ode45(odefun, tspan, y0)其中，odefun为自定义的求解函数，tspan为时间范围，y0为初值条件。
ode45函数通过自适应步长Runge-Kutta法对ODE进行求解，得到时间t与对应的解y。
2. 参数设定在使用ode45函数时，需要设定求解的ODE函数odefun、时间范围tspan和初值条件y0。
还可以设定相对误差容限RelTol和绝对误差容限AbsTol，以控制求解的精度。
还可以通过设定事件函数、输出函数等对求解过程进行控制和监测。
3. 示例应用下面通过一个简单的例子来说明ode45函数的应用。
考虑一个简单的一阶ODE问题：dy/dt = -y，y(0) = 1其MATLAB代码如下：```matlabodefun = @(t, y) -y;tspan = [0, 5];y0 = 1;[t, y] = ode45(odefun, tspan, y0);plot(t, y, '-o');```这段代码首先定义了ODE函数odefun，然后设定时间范围tspan和初值条件y0，最后利用ode45函数求解ODE并绘制解y随时间t的图像。
4. 个人观点和理解在我看来，ode45函数是MATLAB中非常强大且常用的求解ODE的工具。
它能够高效地对各种类型的ODE进行求解，并且通过设定参数和监控函数，可以实现对求解过程的精确控制。
在工程领域和科学研究中，ode45函数的灵活性和高效性使其成为了不可或缺的工具。
第4章 多自由度系统的振动题解
62 / 2962习 题4-1 在题3-10中，设m 1=m 2=m ，l 1=l 2=l ，k 1=k 2=0，求系统的固有频率和主振型。
解：由题3-10的结果22121111)(l g m l g m m k k +++=，2221l gm k -=，2212l g m k -=，22222l gm k k += 代入m m m ==21，021==k k ，l l l ==21可求出刚度矩阵K 和质量矩阵M⎥⎦⎤⎢⎣⎡=m m M 00；⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=l mg lmg l mg l mg K 3 由频率方程02=-M p K ，得0322=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----=mp l mg l mg lmgmp l mg B 0242222242=+-∴l g m p l g m p ml g p )22(1-=∴ ，lgp )22(2+= 为求系统主振型，先求出adjB 的第一列⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=l mg mp lmg adjB 2分别将频率值21p p 和代入，得系统的主振型矩阵为题4-1图63 / 2963⎥⎦⎤⎢⎣⎡-=112)1(A ⎥⎦⎤⎢⎣⎡+=112)2(A4-2 题4-2图所示的均匀刚性杆质量为m 1，求系统的频率方程。
解：设杆的转角θ和物块位移x 为广义坐标。
利用刚度影响系数法求刚度矩阵k 。
设0,1==x θ，画出受力图，并施加物体力偶与力2111,k k ，由平衡条件得到，222111a k b k k +=， a k k 221-=设1,0==x θ，画出受力图，并施加物体力偶与力2212,k k ，由平衡条件得到，12k a k 2-=， a k k 222=得作用力方程为⎥⎦⎤⎢⎣⎡=⎭⎬⎫⎩⎨⎧⎥⎦⎤⎢⎣⎡--++⎭⎬⎫⎩⎨⎧⎥⎥⎦⎤⎢⎢⎣⎡0000312222221221x a k a k a k a k b k x m a m θθ由频率方程02=-M K p ，得031222222212221=----+p m a k ak a k p a m a k b k4-3 题4-3图所示的系统中，两根长度题4-3图题4-2图64 / 2964为l 的均匀刚性杆的质量为m 1及m 2，求系统的刚度矩阵和柔度矩阵，并求出当m 1=m 2=m 和k 1=k 2=k 时系统的固有频率。
关于ode45函数的说明
y0 y 1 x y2 ... yn 1
Dx是一个数组，其第一个元素Dx(1)表示dy0/dt 第二个元素Dx(2)表示dy1/dt 如此类推。
x也是一个数组，其第一个元素x(1)表示y0 第二个元素x(2)表示y1 如此类推。 根据这个关系，很容易可以把一阶微分方程组用 数组Dx与数组x的元素表示出来。例如（见下页）
dy0 f 0 t , y0 , y1 , y2 ,..., yn 1 dt dy1 f1 t , y0 , y1 , y2 ,..., yn 1 dt dy2 f 2 t , y0 , y1 , y2 ,..., yn 1 dt dyn 1 f n 1 t , y0 , y1 , y2 ,..., yn 1 dt
[T,Y]=ode45(‘odefun’,tspan,y0)
tspan为求解区间 y0为初始条件
（1）根据问题所属学科中的规律、定律、公式， 用微分方程与初始条件进行描述。
F y , y , y ,..., y ,t 0 微分方程
n
初始条件 y 0 c0 , y 0 c1 ,..., y
（3）根据（1）和（2）的结果，编写计算导数的 M函数文件。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Ode45函数调用形式如下：
[T,Y]=ode45(odefun,tspan,y0)
相关参数介绍如下：参数名称
参数说明odefun
用于存放待求解的方程的m 文件名，方程必须用y’=f(t,y)的形式存放tspan
指定自变量范围的向量，通常用[t0，tf]指定y0函数的边界条件，即y0=y(t0)，对于方程组，y0也可以是向量
例：若一三自由度多体动力学系统方程如下：
1121221231233232323 1.510050 2.0sin(3.754t)
2 1.5
3 1.55010050 2.0cos(2.2t)
2 1.5350100 1.0sin(2.8t)x x x x x x x x x x x x x x x x x +-+-=-+--+-=-++-+=
初始条件：
1020301020301
1x x x x x x ======
由于方程必须用y’=f(t,y)的形式存放，因此需要对方程组进行降阶处理。
令11
3253214263y x y x y x y x y x y x ======
则方程组可化为：
12
2241334
424613556
646350.5*(3 1.510050 2.0sin(3.754t))
0.5*(1.53 1.55010050 2.0cos(2.2t))
0.5*(1.5350100 1.0sin(2.8t))y y y y y y y y y y y y y y y y y y y y y y y ==-+-++==-++-+-==--+-+
因此建立M函数文件来定义此方程组如下：
function dy=func(t,y)
dy=zeros(6,1);
dy(1)=y(2);
dy(2)=0.5*(-3*y(2)+1.5*y(4)-100*y(1)+50*y(3)+2.0*sin(3.754*t));
dy(3)=y(4);
dy(4)=0.5*(1.5*y(2)-3*y(4)+1.5*y(6)+50*y(1)-100*y(3)+50*y(5)-2.0*cos(2.2*t)); dy(5)=y(6);
dy(6)=0.5*(-1.5*y(4)-3*y(6)+50*y(3)-100*y(5)+1.0*sin(2.8*t));
end
在matlab命令窗口里输入一下命令：
y0=[111111];
tspan=[030];
[t,y]=ode45(@func,tspan,y0);
figure(1)
plot(t,y(:,1),t,y(:,3),t,y(:,5));
legend('x1','x2','x3');
xlabel('时间(s)','FontSize',10);
ylabel('振动位移曲线','FontSize',10);
figure(2)
plot(t,y(:,2),t,y(:,4),t,y(:,6));
legend('v1','v2','v3');
xlabel('时间(s)','FontSize',10);
ylabel(‘振动速度曲线’,’FontSize’,10);
运行结果：
0510********-1-0.5
0.51
1.5
时间（s ）振动位移曲线
0510********-6-5
-4
-3
-2-1012
3
时间（s ）振动速度曲
线。