使用符号数学工具箱金宝app验证simulink模型

打开实时脚本

This example shows how to model a bouncing ball, which is a classical hybrid dynamic system. This model includes both continuous dynamics and discrete transitions. It uses the Symbolic Math Toolbox™ to help explain some of the theory behind ODE solving in theSimulink® Model of a Bouncing Ball.

Assumptions

The ball rebounds with no angle
There is no drag
Height at time t=0 is 10 m
以15 m/s的速度抛出

Derivation

The coefficient of restitution is defined as

$C_{r} = v_{b} - v_{a} / u_{a} - u_{b}$

where v is the speed of object before impact and u is the speed after.

We split the second order differential equation

$\frac{d^{2} h}{d t^{2}} = - g$

into

$\frac{d h}{d t} = v$ 离散为 $\frac{h (t + δ t) - h (t)}{δ t} = v (t)$

和

$\frac{d v}{d t} = - g$ 离散为 $\frac{v (t + δ t) - v (t)}{δ t} = - g$

We will use basic 1st order numerical integration using forward Euler.

$h (t + δ t) = h (t) + v (t) δ t$

$v (t + δ t) = v (t) - g δ t$

分析解决问题

Using the Symbolic Math Toolbox, we can approach the problem analytically. This allows us to solve certain problems more easily, such as determining precisely when the ball first hits the ground (below).

Declare our symbolic variables.

symsgtH(t)h_0v_0

拆分第二阶微分方程 $\frac{d^{2} h}{d t^{2}} = - g$ into $\frac{d h}{d t} = v$ 和 $\frac{d v}{d t} = - g$ .

Dh = diff(H); D2h = diff(H, 2) == g

d2h（t）=
                  
                    $\frac{\partial^{2}}{\partial t^{2}} H (t) = g$

Solve the ODE usingDSOLVE:

eqn = dsolve(D2h, H(0) == h_0, Dh(0) == v_0)

eqn =
                  
                    $\frac{g t^{2}}{2} + v_{0} t + h_{0}$

参数探索使用运动的抛物线subs:

eqn = subs(eqn, [h_0, v_0, g], [10, 15, -9.81])

eqn =
                  
                    $- \frac{981 t^{2}}{200} + 15 t + 10$

通过求解零来找到球击中地面的时间：

假设（t> 0）thit = solve（eqn == 0）

tHit =
                  
                    $\frac{20 \sqrt{5} \sqrt{26}}{109} + \frac{500}{327}$

可视化解决方案：

fplot(eqn,[0 10]) ylim([0 25])

图包含一个轴对象。轴对象包含一个类型函数线的对象。

使用可变精度算术格式化您的确切结果vpa:

disp(['The ball with an initial height of 10m and velocity of 15m/s will hit the ground at 'char(vpa(tHit, 4))' seconds.'）））

The ball with an initial height of 10m and velocity of 15m/s will hit the ground at 3.621 seconds.

Solve the Problem Numerically

Setup Simulation Parameters

球的特性

c_bounce = .9;% Bouncing's coefficient of restitution; perfect restitution would be 1

Properties of the simulation

重力= 9.8;％重力加速度（M/S）height_0 = 10;% Initial height at time t=0 (m)Velocity_0 = 15;% Initial velocity at time t=0 (m/s)

Declaring the simulation timestep

dt = 0.05;％动画时间段（S）t_final = 25;％模拟周期t = 0:dt:t_final;% TimespanN = length(t);% Number of iterations

Initialize simulation quantities

h = [];% Height of the ball as a function of time (m)v = [];% Velocity of the ball (m/sec) as a function of time (m/s)h(1)=height_0; v(1)=velocity_0;

Simulate the bouncing ball (we will use basic 1st order numerical integration using forward Euler):

fori=1:N-1 v(i+1)=v(i)-gravity*dt; h(i+1)=h(i)+v(i)*dt;% When the ball bounces (the height is less than 0),% reverse the velocity and recalculate the position.% Using the coefficient of restitutionifh（i+1）<0 v（i）=  -  v（i）*c_bounce;v（i+1）= v（i）-gravity*dt;h（i+1）= 0+v（i）*dt;endend

Visualize and validate the simulation

情节（t，h，，'o'） 抓住在fplot(eqn,[0 10]) title('Height over time'25) ylim ([0])离开

图包含一个轴对象。The axes object with title Height over time contains 2 objects of type line, functionline.

plot(t,v) title('Velocity over time')

图包含一个轴对象。The axes object with title Velocity over time contains an object of type line.

通过分析验证数字

将您的分析结果与数字结果进行比较。

As a reminder the time of impact was:

disp(['The ball with an initial height of 10m and velocity of 15m/s will hit the ground at 'char(vpa(tHit, 4))' seconds.'）））

The ball with an initial height of 10m and velocity of 15m/s will hit the ground at 3.621 seconds.

From the numerical simulation we can find the closest value in the simulation when $h (t_{i}) \approx 0$

i = ceil（double（thit/dt））;t（[I-1 I I+1]）

ans =1×33.5500 3.6000 3.6500

情节（t，h，，'o'） 抓住在fplot(eqn,[0 10]) plot(t([i-1 i i+1 ]),h([i-1 i i+1 ]),'*R'） 标题（'Height over time') xlim([0 5]) ylim([0 25]) hold离开

图包含一个轴对象。带有标题高度的轴对象包含3个类型线，函数线的对象。