Matlab tutorial
12-13/09/2019
Ing. Marianna Magni
TOPICS:
• Numerical differentiation
• Numerical integration
• Ordinary differential equations
2
METHOD FOR PARTIAL DIFFERENTIAL EQUATION
To solve a partial derivative in several variables, the function diff
allows to compute the derivation with respect to one variable.
𝑄𝑠
𝑓𝑢𝑛𝑐𝑡𝑖𝑜𝑛: 𝑘 = 2
𝑙 (𝑇2 − 𝑇1 )
diff function, variable
𝑑𝑖𝑓𝑓 𝑘, 𝑄
𝛿𝑘 𝑠
𝑃𝑎𝑟𝑡𝑖𝑎𝑙 𝑑𝑒𝑟𝑖𝑣𝑎𝑡𝑖𝑣𝑒: =
𝛿𝑄 𝑙 2 (𝑇2 − 𝑇1 )
3
CASE STUDY 1:
COMPUTATION OF THE MEASUREMENT UNCERTAINTY
For measuring the thermal conductivity of a plate the relation 𝑘 =
𝑄𝑠
2 is used. In steady conditions the power measurement is 𝑄 =
𝑙 (𝑇2 −𝑇1 )
500 𝑚𝑊 ± 0.5% 𝐿𝐶 99.7% . The plate has a 𝑠 = 9.82 𝑚𝑚 thickness,
measured with a centesimal gauge, and the heat exchange area A is a
square whose side measured 30 times with a twentieths caliber, has
provided an average value of 0.12567 𝑚 and a standard deviation of
1.2 𝑚𝑚 . The measurements of the temperatures 𝑇1 and 𝑇2 are
respectively 20°C and 40°C obtained with a thermometer having a
resolution 0.01 K. Compute the measurement uncertainty of the thermal
conductivity.
4
CASE STUDY 1: COMPUTATION OF THE UNCERTAINTY
0.5%𝑄
𝑢𝑄 = = 8.33 10−4 𝑊
3
𝑄𝑠 0.01
𝑘= 2 𝑢𝑠 = = 2.88 10−3 𝑚𝑚
𝑙 (𝑇2 − 𝑇1 ) 2 3
𝑠𝑡
𝑢𝑙 = = 0.219 𝑚𝑚
30
0.01
𝑢𝑇 = = 2.88 10−3 𝐾
2 3
Uncertainty propagation:
2 2 2 2 2
𝛿𝑘 𝛿𝑘 𝛿𝑘 𝛿𝑘 𝛿𝑘
𝑢𝑘 = 𝑢 + 𝑢 + 𝑢 + 𝑢 + 𝑢
𝛿𝑄 𝑄 𝛿𝑠 𝑠 𝛿𝑙 𝑙 𝛿𝑇2 𝑇2 𝛿𝑇1 𝑇1
5
CASE STUDY 1: COMPUTATION OF THE UNCERTAINTY
𝑢𝑄 = 8.33 10−4 𝑊
𝑢𝑠 = 2.88 10−3 𝑚𝑚
𝑢𝑙 = 0.219 𝑚𝑚
𝑢 𝑇 = 2.88 10−3 𝐾
𝑊
𝑘 = 1.55449 10−5
𝑚𝑚 °𝐶
2 2 2 2 2
𝛿𝑘 𝛿𝑘 𝛿𝑘 𝛿𝑘 𝛿𝑘
𝑢𝑘 = 𝑢 + 𝑢 + 𝑢 + 𝑢 + 𝑢
𝛿𝑄 𝑄 𝛿𝑠 𝑠 𝛿𝑙 𝑙 𝛿𝑇2 𝑇2 𝛿𝑇1 𝑇1
2 2
= 3.109 10−5 ∙ 8.33 10−4 2 + 1.58 10−6 ∙ 2.88 10−3 2 + −2.474 10−7 ∙ 0.219 2 + −7.77 10−7 ∙ 2.88 10−3 + 7.77 10−7 ∙ 2.88 10−3
𝑊
= 6.03 10−8
𝑚𝑚 °𝐶
6
CASE STUDY 1: MATLAB SCRIPT
%% 1-DEFINE THE FUNCTION
syms Q s l T1 T2 % syms: declare the symbolic variable and function
….
pretty(k) % pretty: prints the symbolic expression in a format
% that resembles type-set mathematics.
%% 2-COMPUTE THE PARTIAL DERIVATIVE OF THE FUNCTION k
…..
%% 3-RESULTS SUBSTITUTING THE VALUES INTO THE PARTIAL DERIVATIVES
….
dk/dQ=double(subs(dk/dQ)) % subs: substitute the value into the partial derivative
% double: conversion in a decimal value
7
CASE STUDY 1: REMARK
8
METHOD OF NUMERICAL INTEGRATION
To compute the integration of a general function, Matlab uses the int
function:
𝑖𝑛𝑡(𝑓, 𝑥, 𝑎, 𝑏)
It gives the definite integral of f with respect to x from a to b.
The integration interval can also be specified using a row or a column
vector with two elements, i.e., 𝑖𝑛𝑡 𝑓, 𝑥, 𝑎 𝑏 and 𝑖𝑛𝑡 𝑓, 𝑥, 𝑎; 𝑏 .
9
CASE STUDY 2: THE PROJECTILE MOTION
𝑚
Consider a projectile, launched upward with an initial velocity 𝑣0 = 2
𝑠
which forms an angle 𝜗 = 45° with the horizontal direction. The speed in
the Cartesian axes is defined as:
𝑣0𝑥 = 𝑣0 cos 𝜗
𝑣0𝑦 = 𝑣0 sin 𝜗
Along the x axis the speed remains constant, while along the y axis it
varies as:
𝑣𝑦 = 𝑣0 sin 𝜗 − 𝑔𝑡
Define the components
of the projectile motion and
plot the trajectory.
10
CASE STUDY 2: THE PROJECTILE MOTION
Integrating the speed components, the displacement along the x and y
axes is computed as:
𝑣𝑥 𝑑𝑡 → 𝑥 = 𝑡 𝑣0 cos 𝜗
𝑔𝑡 2
𝑣𝑦 𝑑𝑡 → 𝑦 = 𝑡 𝑣0 sin 𝜗 −
2
And finally substituting t, from the first equation, into the second, the
trajectory of the projectile is given by:
sin 𝜗 𝑔 2
𝑦= 𝑥− 2 𝑥
cos 𝜗 2𝑣0 𝑐𝑜𝑠 2 𝜗
11
CASE STUDY 2: MATLAB SCRIPT
%% 1-DEFINE THE FUNCTION
…..
%% 2-NUMERICAL INTEGRATION
y_th=@(xp)(...) % @: is used to create a function handle.
x_vect=….
plot….
…..
%% 3-INTEGRATION OF A SYMBOLIC FUNCTION
….
x=int(vx,t);
y=int(vy,t);
…..
%% 4-PLOT THE TRAJECTORY
…..
12
CASE STUDY 2: THE PROJECTILE MOTION RESULTS
Analytical computation
Matlab computation
13
CASE STUDY 2: REMARK FUNTOOL
14
NUMERICAL METHOD FOR ORDINARY DIFFERENTIAL
EQUATIONS
An Ordinary Differential Equation (ODE) contains
• one or more derivatives of a dependent variables with respect to a
single independent variable usually referred to as time.
The form:
𝑑
𝑦 𝑡 = 𝑓 𝑡, 𝑦 𝑡 , 𝑦 𝑡0 = 𝑦0
𝑑𝑡
can be solved by using the ODE function
15
[T,Y]=solver (odefun,tspan,y0)
with:
• T column vector of time points
• Y solution array. Each row in Y correspond to the solution at a time returned
in the corresponding row of T.
• tspan integrates the system differential equations y’ from time 𝑡0 to 𝑡𝑓 with
initial conditions 𝑦0 .
• odefun is a function handle.
16
EXAMPLE: SOLVING VAN DER POL EQUATION
This example illustrates the steps for solving an initial value ODE problem:
1. Rewrite the problem as a system of first-order ODEs. Rewrite the van der Pol
equation (second-order):
𝑦1′′ − 𝜇 1 − 𝑦12 𝑦1′ + 𝑦1 = 0
The resulting system of first-order ODEs becomes:
𝑦1′ = 𝑦2
𝑦2′ = 𝜇 1 − 𝑦12 𝑦2 − 𝑦1
2. Code the system of first-order ODEs. Code the system as a function that an ODE
solver can use. The function must be of the form:
𝑓𝑢𝑛𝑐𝑡𝑖𝑜𝑛 𝑑𝑦𝑑𝑡 = 𝑓(𝑡, 𝑦)
𝑑𝑦𝑑𝑡 = 𝑦2 ; 𝜇 1 − 𝑦12 𝑦2 − 𝑦1 ;
3. Apply a solver to the problem. Call the solver and pass it the function you created
to describe the first-order system of ODEs, the time interval on which you want to
solve the problem, and an initial condition vector.
𝑡, 𝑦 = 𝑜𝑑𝑒45(@𝑓, 0 20 , 2; 0 )
17
EXAMPLE: SOLVING VAN DER POL EQUATION
4. View the solver output. Use the plot command:
𝑝𝑙𝑜𝑡 (𝑡, 𝑦 : , 1 , ′ −′ , 𝑡, 𝑦 : , 2 ,′ −−′ )
18
EXAMPLE: REMARK
19
CASE STUDY 3: 1 DOF SYSTEM SPRING MASS DAMPER SYSTEM
x F(t) 𝑚 = 1 𝑘𝑔
m
𝑘 = 10 𝑁/𝑚
𝑟 = 0.1 𝑁𝑠/𝑚
r
k 𝑚𝑥 + 𝑟𝑥 + 𝑘𝑥 = 𝐹(𝑡)
Tasks:
1. Free Vibration of the system (initial position and null speed) varying
system damping
2. System response with sinusoidal input 𝑓0 = 0.05 Hz
20
CASE STUDY 3: 1 DOF SYSTEM
The equation 𝑚𝑥 + 𝑟𝑥 + 𝑘𝑥 = 𝐹 𝑡 can be written as a
system:
𝑑𝑥
=𝑥
𝑑𝑡
𝑑 𝑥 𝐹(𝑡) 𝑟𝑥 𝑘𝑥
= − −
𝑑𝑡 𝑚 𝑚 𝑚
21
CASE STUDY 3: 1 DOF SYSTEM
To solve the tasks, the numerical integration procedure is applied:
• a separate function file must be generated to define the previous
system of equations. The equations are written in the form of a
vector.
• a code is written, which calls the above function and solves the
differential equation and plots the required result
22
CASE STUDY 3: 1 DOF SYSTEM MATLAB SCRIPT
• function:
function dy=….
dy=[…..]
• script:
%% VARIABLE DEFINITION
…
%% FREE VIBRATION
…
%% SOLVE THE SINGLE DOF SISTEM
…
[…]=ode45(…)
…
%% FIGURES
…
23
CASE STUDY 3: TASK 1
Free Vibration Response r=0.1 Ns/m
24
CASE STUDY 3: TASK 1
Free Vibration Response r=0.01 Ns/m
25
CASE STUDY 3: TASK 2
Transient
Steady-state
26
CASE STUDY 3: REMARK
27
CASE STUDY 4: 2 DOF SYSTEM MECHANICAL SYSTEM
𝐹0 cos(2𝜋𝑓𝑡)
𝑥1 𝑥2 𝑚1 = 1 𝑘𝑔
𝑘 𝑘 𝑘 𝑚2 = 1 𝑘𝑔
𝑚1 𝑚2 𝑘 = 1 𝑁/𝑚
𝐹0 = 10 𝑁
𝑓 = 1 𝐻𝑧
Tasks:
1. Find the dynamic equations
2. Plot the displacement and the speed as function of time
28
CASE STUDY 4: TASK 1
The equations of motion of the system:
𝑑 2 𝑥1
𝑚1 2
+ 𝑘1 + 𝑘2 𝑥1 − 𝑘2 𝑥2 = 𝐹0 cos (2π𝑓𝑡)
𝑑𝑡
𝑑2 𝑥2
𝑚2 2
+ 𝑘2 + 𝑘3 𝑥2 − 𝑘2 𝑥1 = 0
𝑑𝑡
Imposing:
𝑦 1 = 𝑥1
𝑦 2 = 𝑥1
𝑦 3 = 𝑥2
𝑦 4 = 𝑥2
29
CASE STUDY 4: TASK 1
The system becomes:
𝑦 2
𝐹0 cos (2π𝑓𝑡)− 𝑘1 +𝑘2 𝑦 1 +𝑘2 𝑦(3)
𝑑𝑦 𝑚1
=
𝑑𝑡 𝑦 4
− 𝑘2 +𝑘3 𝑦 3 +𝑘2 𝑦(1)
𝑚2
30
CASE STUDY 4: TASK 1
global m1 m2 k1 k2 k3 F0
%% Variables definition
…….
%% Initial condition
……
%% Integration
t=0:0.01:100;
y0=[x10,v10,x20,v20]';
[ti,y]=ode45('linear_system',t,y0);
s1=y(:,1);
v1=y(:,2); function dy=linear_system (t,y)
s2=y(:,3); global m1 m2 k1 k2 k3 F0
v2=y(:,4);
dy=[y(2);
(F0*cos(2*pi*f*t)-(k1+k2)*y(1)+k2*y(3))/m1;
y(4);
(-(k2+k3)*y(3)+k2*y(1))/m2];
31
CASE STUDY 4: TASK 2
CASE 1: imposing 𝐹0 = 0 and the initial conditions equal to 0
32
CASE STUDY 4: TASK 2
𝑥01 =1
𝑥 =0
CASE 2: imposing 𝐹0 = 0 and the initial conditions: 𝑥01 =1
02
𝑥02 =0
33
CASE STUDY 4: TASK 2
𝑥01 = 1
𝑥 =0
CASE 3: imposing 𝐹0 = 0 and the initial conditions: 𝑥 01= −1
02
𝑥02 = 0
34
CASE STUDY 4: TASK 2
CASE 4: imposing 𝑓 = 0.01 𝐻𝑧 (before the natural frequency f=0.16 Hz)
and the initial conditions equal to 0
35
CASE STUDY 4: TASK 2
CASE 5: imposing 𝑓 = 0.16 𝐻𝑧 (the natural frequency f=0.16 Hz)
and the initial conditions equal to 0
36
CASE STUDY 4: TASK 2
CASE 6: imposing 𝑓 = 0.185 𝐻𝑧 (after the natural frequency f=0.16 Hz)
and the initial conditions equal to 0
37
CASE STUDY 4: TASK 2
CASE 7: imposing 𝑓 = 0.274 𝐻𝑧 (the second natural frequency f=0.274 Hz)
and the initial conditions equal to 0
38
CASE STUDY 4: REMARK
39
