HANDOUT ON SOLUTIONS OF DIFFERENTIAL EQUATIONS
MEEN 364
Parasuram July 27, 2001
HANDOUT M.4 - SOLUTIONS OF ORDINARY DIFFERENTIAL EQUATIONS Introduction The majority of physical systems of interest in this course can be represented in the form of ordinary differential equations. In order to study the response of the system modeling, the differential equations must be solved. The general solution of a linear differential equation is the sum of two components, the particular integral and the complementary function. Often the particular integral is the steady-state component of the solution of the differential equation, while the complimentary function, which is the solution of the corresponding homogenous equation, is the transient component of the solution. Often the steady state component of the solution has the same form as the driving or (external) function. Solution of a second order differential equation The model of a majority of mechanical systems can be approximated with a second order differential equation. So the solution of a second order differential equation is important, as it can be used to study the response of many physical systems. Case 1: Constant forcing function Consider, for example a mass spring damper system subjected to a constant force ‘F’, as shown in the following figure: F m
k
c
The governing differential equation of motion of the mass, ‘m’ is given by ..
.
m x(t ) + c x(t ) + kx(t ) = F ,
(1)
where dividing the entire equation by ‘m’, we have
1
MEEN 364
Parasuram July 27, 2001
Note that the time dependence is assumed and hence it is omitted for future manipulations. ..
x+
c . k F x+ x = . m m m
(2)
Let c = 2ζ ω n , m k = ω n2 . m So equation (2) reduces to ..
.
x + 2ζ ω n x + ω n2 x =
F . m
(3)
As explained earlier, the solution to the above differential equation has two parts, the particular integral and the complimentary function. To find the particular integral, which is nothing but the steady-state solution, equate all of the derivate terms to zero, i.e., ..
x = 0, .
x = 0. Substituting the above relation in equation (3), we have
ω n2 x =
F m
F ⇒ xs = 2 ωn m
(4)
The second of equations (4) represents the steady-state response of the system. To find the complimentary function or the transient response, we need to consider the homogenous equation, i.e., ..
.
x + 2ζ ω n x + ω n2 x = 0 .
(5)
Equation (5) is the left hand side of equation (3), with the right hand side set to zero. To solve the above equation, first assume a solution of the form
2
MEEN 364
Parasuram July 27, 2001
x = Ae λt .
(6)
Substituting the above relation in equation (5), we get 2
Aλ e λt + 2 Aζ ω n λ e λt + Aω n2 e λt = 0 ⇒ ( Aλ
2
+ 2 Aζ ω n λ + Aω n2 )e λt = 0
(7) (8)
e λt = 0 is a trivial solution, as equation (6) will then be identically equal to zero, i.e., x ≡ 0 . But we need to find the non-trivial solution, so e λt ≠ 0 . So from equation (8), we have 2
λ + 2ζ ω n λ + ω n2 = 0
(9)
Solving the above quadratic equation, we find the following roots 2
λ 1, 2 =
− 2ζ ω n ± 4ζ ω n2 − 4ω n2 2
⇒λ
1, 2
= −ζ ω n ± ω n ζ
⇒λ
1, 2
= −ζ ω n ± ω d
2
−1
where
ωd = ωn ζ
2
− 1.
Therefore the transient solution of the system is given as xt = Ae
( −ζ ω n +ω d ) t
+ Be
( −ζ ω n −ω d ) t
(10)
where A and B can be found based on the initial conditions. So the final solution of the differential equation given by equation (1) is x = x s + xt ⇒x=
F ( −ζ ω n +ω d ) t ( −ζ ω n −ω d ) t + Ae + Be 2 mω n
For an unforced system, the differential equation is given by ..
.
m x + c x + kx = 0 .
(11)
3
MEEN 364
Parasuram July 27, 2001
Since the forcing function is equal to zero, here the steady state response of the system will be equal to zero, while the transient response will be the same as that of the forced system. In other words the solution of the equation (11) is given as x = x s + xt ⇒ x = 0 + Ae ⇒ x = Ae
( −ζ ω n +ω d ) t
( −ζ ω n +ω d ) t
+ Be
+ Be
( −ζ ω n −ω d ) t
( −ζ ω n −ω d ) t
Case 2: Time varying forcing function In this case, instead of a constant forcing function driving the system as in case 1, let us consider a harmonic function, Fejwt, which by Euler’s expansion is equal to F (cos ω t + j sin ω t ) , driving the system. ’w’ is the frequency with which the driving function is applied on the system. The governing differential equation will be ..
.
m x + c x + kx = Fe jwt ..
.
⇒ x + 2ζ ω n x + ω n2 x =
F jwt e = Fo e jwt . m
(12)
Usually the form of the solution is similar to that of the driving function. In this case, to obtain the steady state component, we assume a solution of the form x = Ce jwt . Substituting the above relation in equation (12), we have 2
− Cω e jwt + 2Cζ ω n jω e jwt + Cω n2 e jwt = F0 e jwt 2
⇒ {(ω n2 − ω ) + 2 jζ ω nω }C = F0
(13)
Diving the second equation of equation (13) by ω n2 , results in F0
C=
F ω n2 k = , 2 2 {(1 − r ) + 2ζ r} {(1 − r ) + 2ζ r}
where ω r= . ωn Once the value of ‘C’ is obtained, the steady state solution of the system is x s = Ce jwt .
4
MEEN 364
Parasuram July 27, 2001
The complimentary function or the transient response of this system is the same as that of the previous case, as the homogenous equation is the same. So the entire solution to the system subjected to a harmonic driving function is given by x = xt + x s ⇒ x = Ae where A
( −ζ ω n +ω d ) t
and
+ Be
( −ζ ω n −ω d ) t
B − − − obtained
+ Ce jwt ,
from initial
conditions
F
C=
k , {(1 − r ) + 2ζ r} 2
where ω r= . ωn Using MATLAB to solve the differential equations The ‘ODE’ command in MATLAB can be used to solve various differential equations. The use of this command is explained in detail with the help of the following example. Example 1 Consider the system below. m
k
c
This system is similar to the one in case 1, except that the forcing function is equal to zero. Therefore the governing differential equation of motion of the system is ..
.
m x + c x + kx = 0 .
(14)
5
MEEN 364
Parasuram July 27, 2001
Problem statement Solve for the response of the unforced system defined by the differential equation given by equation (14). Assume ξ = 0.1; m = 1Kg; k = 100 N/m. Let the initial conditions be chosen as x(0) = 0.02 m .
x(0) = 0 m / s Solution Equation (14) is a second order constant-coefficient differential equation. To solve this equation we have to reduce it into two first order differential equations. This step is taken because MATLAB uses a Runge-Kutta method to solve differential equations, which is valid only for first order differential equations. Let .
x = v.
(15)
From the above expression we see that .
m v + cv + kx = 0 ,
(16)
so equation (16) reduces to .
v = [(
−c k )v − ( ) x]. m m
(17)
We can see that the second order differential equation (14) has been reduced to two first order differential equations (15) and (17). For our convenience, put x = y (1); .
x = v = y (2); Note that, the variables ‘x’ and ‘v’ are the states of the system, and are each equated to another variable y(1) and y(2) respectively to be compatible with MATLAB. These variables are not the outputs of the system.
6
MEEN 364
Parasuram July 27, 2001
Equations (15) and (17) reduce to .
y (1) = y (2);
(18)
.
y (2) = [(-c/m)*y (2) – (k/m)*y (1)];
(19)
To calculate the value of ‘c’, compare equation (14) with the following generalized equation. ..
.
2
x + 2ξω n x + ω n = 0 . Equating the coefficients of the similar terms we have c = 2ξω n m k 2 ωn = . m
(20) (21)
Using the values of ‘m’ and ‘k’, calculate the value of ‘c’ corresponding to the value of ξ. Once the value of ‘c’ is known, equations (18) and (19) can be solved using MATLAB. Some of the properties of the system are defined as follows. Natural frequency ωn = For ξ = 0.1
(k / m) = 10 rad/sec.
Damped natural frequency ωd = ωn 1 − ξ
2
= 9.95 rad/sec.
Damped time period Td = 2π/ωd = 0.63 sec. Therefore for to plot five time cycles, the time interval should be 5 times the damped time period, i.e., 3.16 sec. MATLAB Code In order to apply the ODE45 or any other numerical integration procedure, a separate function file must be generated to define equations (18) and (19). Actually, the right hand side of equations (18) and (19) are what is stored in the file. The equations are written in the form of a vector. The MATLAB code is given below. function yp = unforced(t,y) yp = [y(2);(-((c/m)*y(2))-((k/m)*y(1)))];
(22)
In the above code ‘yp’ is a two dimensional vector representing the right hand side of the equations (18) and (19), respectively. Open a new M-file and write down the above two
7
MEEN 364
Parasuram July 27, 2001
lines. The first line of the function file must start with the word “function” and the file must be saved corresponding to the function call; i.e., in this case, the file is saved as unforced.m. The derivatives are stored in the form of a vector. This example problem has been solved for ξ = 0.1. We need to find the value of ‘c/m’ and ‘k/m’ so that the values can be substituted in equation (22). Substituting the values of ξ and ωn in equations (20) and (21) the values of ‘c/m’ and ‘k/m’ can be found out. After finding the values, substitute them into equation (22). In this example c = 2 × 0.1 × 10 = 2, m k 100 = = 100. m 1 Now we need to write a code, which calls the above function and solves the differential equation and plots the required result. First open another M-file and type the following code. tspan=[0 4]; y0=[0.02;0]; [t,y]=ode45('unforced',tspan,y0); figure(1); plot(t,y(:,1)); grid on xlabel(‘time’) ylabel(‘Displacement’) title(‘Displacement Vs Time’)
The ode45 command in the main body calls the function ‘unforced’, which defines the systems first order derivatives. The response is then plotted using the plot command. ‘tspan’ represents the time interval and ‘y0’ represents the initial conditions for y(1) and y(2), which in turn represent the displacement ‘x’ and the first derivative of ‘x’, i.e., the velocity. In this example, the initial conditions are taken as 0.02 m for ‘x’ and 0 m/sec for the first derivative of ‘x’ or the velocity of the mass ‘m’. MATLAB uses a default step value for the vector ‘tspan’. In order to change the step size use the following code tspan=0:0.001:4;
This tells MATLAB that the vector ‘tspan’ is from 0 to 4 with a step size of 0.001. This example uses an arbitrary step size. If the step size is too large, the plot obtained will not be a smooth curve. So it is always better to use relatively smaller step size. But it also takes longer to simulate programs with smaller step size. So you have to decide which is the best step size you can use for a given problem.
8
MEEN 364
Parasuram July 27, 2001
In the above code ‘y(:,1) represents the displacement ‘x’. To plot the velocity, change the variable in the plot command line to ‘y(:,2)’. The plot is attached below
Note: Any nth order differential equation can be solved using MATLAB by following the procedure outlined in example 1. In other words, the nth order differential equation can be reduced to ‘n’ first order differential equations and the ‘ode’ command in MATLAB can be used to solve the system of ‘n’ first order differential equations.
9
MEEN 364
Parasuram July 27, 2001
Example 2. Double Pendulum Derive the governing differential equation of motion of a double pendulum. Solve the derived differential equation and plot the values of θ1 and θ2 with respect to time. Choose l1 = 1m; l2 = 1m; m1 = m2 = 5Kg. The initial conditions for θ1 and θ2 are 0.5233 and 0.5233 radians respectively.
θ1 l1, m1 θ2
l2 , m2
Solution When represented in a matrix form, the differential equations of motion for the above system is . 2 .. m 2 l 2 cos(θ 2 − θ 1 ) θ 1 − ( w1 + w2 ) sin θ 1 + m 2 l 2 θ 2 sin(θ 2 − θ 1) (m 2 + m1 )l1 m l cos(θ − θ ) .. = . m2 l 2 2 1 21 θ 2 − w sin θ − m l θ 1 2 sin(θ − θ 2 2 1 1 2 1)
MATLAB Code These coupled second order differential equations must be converted into a vector of first order differential equation. This step is done by the following substitutions.
θ 1 = y (1); .
θ 1 = y (2); θ 2 = y (3); .
θ 2 = y (4); Substituting the above relations in the original nonlinear differential equation, we get the following first order differential equations, . y (2) 0 0 0 1 y (1) 2 . 0 . ( m + m ) l 0 m l cos( θ − θ ) − ( w + w ) sin θ + m l θ sin( θ − θ ) y ( 2 ) 2 1 1 2 2 2 1 1 2 1 2 2 2 2 1 = . 0 y (3) 0 1 0 y (4) . . 10 2 m2 l 2 0 m2 l1 cos(θ 2 − θ 1 ) 0 y (4) − w2 sin θ 2 − m1l1 θ 1 sin(θ 2 − θ 1 )
MEEN 364
Parasuram July 27, 2001
In this type of a problem where the coefficient matrix is a function of the states or the variables, a separate M-file must be written which incorporates a switch/case programming with a flag case of ‘mass’. For example if the differential equation is of the form, .
M (t , y ) y (t ) = F (t , y ) then the right hand side of the above equation has to be stored in a separate m-file called ‘F.m’. Similarly the coefficient matrix should be stored in a separate m-file named ‘M.m’. So, when the flag is set to the default, the function ‘F.m’ is called and later when the flag is set to ‘mass’ the function ‘M.m’ is called. The code with the switch/case is given below. Note that it is a function file and should be saved as ‘indmot_ode.m’ in the current directory. function varargout=indmot_ode(t,y,flag) switch flag case '' %no input flag varargout{1}=pend(t,y); case 'mass' %flag of mass calls mass.m varargout{1}=mass(t,y); otherwise error(['unknown flag ''' flag '''.']); end
To store the right hand side of the state variable matrix form of the model, a separate function file must be written as shown below. Note that the name of the function is ‘pend’, so this file must be saved as ‘pend.m’. %the following function contains the right hand side of the %differential equation of the form %M(t,y)*y'=F(t,y) %i.e. it contains F(t,y).it is also stored in a separate file named, pend.m. function yp= pend(t,y) M1=5; M2=5; g=9.81; l1=1; l2=1; w2=M2*9.81; w1=M1*9.81; yp=zeros(4,1); yp(1)=y(2); yp(2)=-(w1+w2)*sin(y(1))+M2*l2*(y(4)^2)*sin(y(3)-y(1)); yp(3)=y(4); yp(4)=-w2*sin(y(3))-M2*l1*(y(2)^2)*sin(y(3)-y(1));
11
MEEN 364
Parasuram July 27, 2001
Similarly, to store the coefficient matrix a separate function file is written which is stored as ‘mass.m’. % the following function contains the mass matrix. %it is separately stored in a file named, mass.m function m = mass(t,y) M1=5; M2=5; g=9.81; l1=1; l2=1; m1=[1 0 0 0]; m2=[0 (M1+M2)*l1 0 M2*l2*cos(y(3)-y(1))]; m3=[0 0 1 0]; m4=[0 M2*l1*cos(y(3)-y(1)) 0 M2*l2]; m=[m1;m2;m3;m4];
To plot the response, the main file should call the function ‘indmot_ode.m’, which has the switch/case programming which in turn calls the corresponding functions depending on the value of the flag. For the main file to recognize the coefficient matrix, the MATLAB command ODESET is used to set the mass to ‘M (t, y)’. The MATLAB code for the main file is given below % this is the main file, which calls the equations and solves using ode113. %it then plots the first variable. tspan=[0 10] y0=[0.5233;0;0.5233;0] options=odeset('mass','M(t,y)') [t,y]=ode113('indmot_ode',tspan,y0,options) subplot(2,1,1) plot(t,y(:,1)) grid xlabel('Time') ylabel('Theta1') subplot(2,1,2) plot(t,y(:,3)) grid xlabel('Time') ylabel('Theta2')
For help regarding ‘subplot’ use the online help in MATLAB by typing the following statement at the command prompt. help subplot subplot(m,n,p), breaks the Figure window into an m-by-n matrix of small axes and selects the p-th axes for the current plot. The axes are counted along the top row of the Figure window, then the second row, etc. For example, 12
MEEN 364
Parasuram July 27, 2001
subplot(2,1,1), plot(income) subplot(2,1,2), plot(outgo) plots “income” on the top half of the window and “outgo” on the bottom half. The plot has been attached.
Assignment 1) For the system shown below,
θ 13
MEEN 364
Parasuram July 27, 2001
l
m m = 10 g; l = 5 cms; .
The initial conditions are θ (0) = 90° and θ (0) = 0 a) Derive the governing differential equation. b) Solve the above differential equation obtained analytically and plot the value of theta with respect to time. c) Solve the differential equation using MATLAB and again plot the value of theta with respect to time. Compare the results obtained from (2) and (3). 2) For the system shown below F m
k
c
Using MATLAB, plot the response of a forced system given by the equation ..
.
m x + c x + kx = f sin ω t For ξ = 0.1. .
Take m = 5 kg; k = 1000 N/m; f = 50 N; ω = 4ωn, x(0) = 5 cms; x (0) = 0.
14
MEEN 364
Parasuram July 27, 2001
15
Parasuram July 27, 2001
HANDOUT M.4 - SOLUTIONS OF ORDINARY DIFFERENTIAL EQUATIONS Introduction The majority of physical systems of interest in this course can be represented in the form of ordinary differential equations. In order to study the response of the system modeling, the differential equations must be solved. The general solution of a linear differential equation is the sum of two components, the particular integral and the complementary function. Often the particular integral is the steady-state component of the solution of the differential equation, while the complimentary function, which is the solution of the corresponding homogenous equation, is the transient component of the solution. Often the steady state component of the solution has the same form as the driving or (external) function. Solution of a second order differential equation The model of a majority of mechanical systems can be approximated with a second order differential equation. So the solution of a second order differential equation is important, as it can be used to study the response of many physical systems. Case 1: Constant forcing function Consider, for example a mass spring damper system subjected to a constant force ‘F’, as shown in the following figure: F m
k
c
The governing differential equation of motion of the mass, ‘m’ is given by ..
.
m x(t ) + c x(t ) + kx(t ) = F ,
(1)
where dividing the entire equation by ‘m’, we have
1
MEEN 364
Parasuram July 27, 2001
Note that the time dependence is assumed and hence it is omitted for future manipulations. ..
x+
c . k F x+ x = . m m m
(2)
Let c = 2ζ ω n , m k = ω n2 . m So equation (2) reduces to ..
.
x + 2ζ ω n x + ω n2 x =
F . m
(3)
As explained earlier, the solution to the above differential equation has two parts, the particular integral and the complimentary function. To find the particular integral, which is nothing but the steady-state solution, equate all of the derivate terms to zero, i.e., ..
x = 0, .
x = 0. Substituting the above relation in equation (3), we have
ω n2 x =
F m
F ⇒ xs = 2 ωn m
(4)
The second of equations (4) represents the steady-state response of the system. To find the complimentary function or the transient response, we need to consider the homogenous equation, i.e., ..
.
x + 2ζ ω n x + ω n2 x = 0 .
(5)
Equation (5) is the left hand side of equation (3), with the right hand side set to zero. To solve the above equation, first assume a solution of the form
2
MEEN 364
Parasuram July 27, 2001
x = Ae λt .
(6)
Substituting the above relation in equation (5), we get 2
Aλ e λt + 2 Aζ ω n λ e λt + Aω n2 e λt = 0 ⇒ ( Aλ
2
+ 2 Aζ ω n λ + Aω n2 )e λt = 0
(7) (8)
e λt = 0 is a trivial solution, as equation (6) will then be identically equal to zero, i.e., x ≡ 0 . But we need to find the non-trivial solution, so e λt ≠ 0 . So from equation (8), we have 2
λ + 2ζ ω n λ + ω n2 = 0
(9)
Solving the above quadratic equation, we find the following roots 2
λ 1, 2 =
− 2ζ ω n ± 4ζ ω n2 − 4ω n2 2
⇒λ
1, 2
= −ζ ω n ± ω n ζ
⇒λ
1, 2
= −ζ ω n ± ω d
2
−1
where
ωd = ωn ζ
2
− 1.
Therefore the transient solution of the system is given as xt = Ae
( −ζ ω n +ω d ) t
+ Be
( −ζ ω n −ω d ) t
(10)
where A and B can be found based on the initial conditions. So the final solution of the differential equation given by equation (1) is x = x s + xt ⇒x=
F ( −ζ ω n +ω d ) t ( −ζ ω n −ω d ) t + Ae + Be 2 mω n
For an unforced system, the differential equation is given by ..
.
m x + c x + kx = 0 .
(11)
3
MEEN 364
Parasuram July 27, 2001
Since the forcing function is equal to zero, here the steady state response of the system will be equal to zero, while the transient response will be the same as that of the forced system. In other words the solution of the equation (11) is given as x = x s + xt ⇒ x = 0 + Ae ⇒ x = Ae
( −ζ ω n +ω d ) t
( −ζ ω n +ω d ) t
+ Be
+ Be
( −ζ ω n −ω d ) t
( −ζ ω n −ω d ) t
Case 2: Time varying forcing function In this case, instead of a constant forcing function driving the system as in case 1, let us consider a harmonic function, Fejwt, which by Euler’s expansion is equal to F (cos ω t + j sin ω t ) , driving the system. ’w’ is the frequency with which the driving function is applied on the system. The governing differential equation will be ..
.
m x + c x + kx = Fe jwt ..
.
⇒ x + 2ζ ω n x + ω n2 x =
F jwt e = Fo e jwt . m
(12)
Usually the form of the solution is similar to that of the driving function. In this case, to obtain the steady state component, we assume a solution of the form x = Ce jwt . Substituting the above relation in equation (12), we have 2
− Cω e jwt + 2Cζ ω n jω e jwt + Cω n2 e jwt = F0 e jwt 2
⇒ {(ω n2 − ω ) + 2 jζ ω nω }C = F0
(13)
Diving the second equation of equation (13) by ω n2 , results in F0
C=
F ω n2 k = , 2 2 {(1 − r ) + 2ζ r} {(1 − r ) + 2ζ r}
where ω r= . ωn Once the value of ‘C’ is obtained, the steady state solution of the system is x s = Ce jwt .
4
MEEN 364
Parasuram July 27, 2001
The complimentary function or the transient response of this system is the same as that of the previous case, as the homogenous equation is the same. So the entire solution to the system subjected to a harmonic driving function is given by x = xt + x s ⇒ x = Ae where A
( −ζ ω n +ω d ) t
and
+ Be
( −ζ ω n −ω d ) t
B − − − obtained
+ Ce jwt ,
from initial
conditions
F
C=
k , {(1 − r ) + 2ζ r} 2
where ω r= . ωn Using MATLAB to solve the differential equations The ‘ODE’ command in MATLAB can be used to solve various differential equations. The use of this command is explained in detail with the help of the following example. Example 1 Consider the system below. m
k
c
This system is similar to the one in case 1, except that the forcing function is equal to zero. Therefore the governing differential equation of motion of the system is ..
.
m x + c x + kx = 0 .
(14)
5
MEEN 364
Parasuram July 27, 2001
Problem statement Solve for the response of the unforced system defined by the differential equation given by equation (14). Assume ξ = 0.1; m = 1Kg; k = 100 N/m. Let the initial conditions be chosen as x(0) = 0.02 m .
x(0) = 0 m / s Solution Equation (14) is a second order constant-coefficient differential equation. To solve this equation we have to reduce it into two first order differential equations. This step is taken because MATLAB uses a Runge-Kutta method to solve differential equations, which is valid only for first order differential equations. Let .
x = v.
(15)
From the above expression we see that .
m v + cv + kx = 0 ,
(16)
so equation (16) reduces to .
v = [(
−c k )v − ( ) x]. m m
(17)
We can see that the second order differential equation (14) has been reduced to two first order differential equations (15) and (17). For our convenience, put x = y (1); .
x = v = y (2); Note that, the variables ‘x’ and ‘v’ are the states of the system, and are each equated to another variable y(1) and y(2) respectively to be compatible with MATLAB. These variables are not the outputs of the system.
6
MEEN 364
Parasuram July 27, 2001
Equations (15) and (17) reduce to .
y (1) = y (2);
(18)
.
y (2) = [(-c/m)*y (2) – (k/m)*y (1)];
(19)
To calculate the value of ‘c’, compare equation (14) with the following generalized equation. ..
.
2
x + 2ξω n x + ω n = 0 . Equating the coefficients of the similar terms we have c = 2ξω n m k 2 ωn = . m
(20) (21)
Using the values of ‘m’ and ‘k’, calculate the value of ‘c’ corresponding to the value of ξ. Once the value of ‘c’ is known, equations (18) and (19) can be solved using MATLAB. Some of the properties of the system are defined as follows. Natural frequency ωn = For ξ = 0.1
(k / m) = 10 rad/sec.
Damped natural frequency ωd = ωn 1 − ξ
2
= 9.95 rad/sec.
Damped time period Td = 2π/ωd = 0.63 sec. Therefore for to plot five time cycles, the time interval should be 5 times the damped time period, i.e., 3.16 sec. MATLAB Code In order to apply the ODE45 or any other numerical integration procedure, a separate function file must be generated to define equations (18) and (19). Actually, the right hand side of equations (18) and (19) are what is stored in the file. The equations are written in the form of a vector. The MATLAB code is given below. function yp = unforced(t,y) yp = [y(2);(-((c/m)*y(2))-((k/m)*y(1)))];
(22)
In the above code ‘yp’ is a two dimensional vector representing the right hand side of the equations (18) and (19), respectively. Open a new M-file and write down the above two
7
MEEN 364
Parasuram July 27, 2001
lines. The first line of the function file must start with the word “function” and the file must be saved corresponding to the function call; i.e., in this case, the file is saved as unforced.m. The derivatives are stored in the form of a vector. This example problem has been solved for ξ = 0.1. We need to find the value of ‘c/m’ and ‘k/m’ so that the values can be substituted in equation (22). Substituting the values of ξ and ωn in equations (20) and (21) the values of ‘c/m’ and ‘k/m’ can be found out. After finding the values, substitute them into equation (22). In this example c = 2 × 0.1 × 10 = 2, m k 100 = = 100. m 1 Now we need to write a code, which calls the above function and solves the differential equation and plots the required result. First open another M-file and type the following code. tspan=[0 4]; y0=[0.02;0]; [t,y]=ode45('unforced',tspan,y0); figure(1); plot(t,y(:,1)); grid on xlabel(‘time’) ylabel(‘Displacement’) title(‘Displacement Vs Time’)
The ode45 command in the main body calls the function ‘unforced’, which defines the systems first order derivatives. The response is then plotted using the plot command. ‘tspan’ represents the time interval and ‘y0’ represents the initial conditions for y(1) and y(2), which in turn represent the displacement ‘x’ and the first derivative of ‘x’, i.e., the velocity. In this example, the initial conditions are taken as 0.02 m for ‘x’ and 0 m/sec for the first derivative of ‘x’ or the velocity of the mass ‘m’. MATLAB uses a default step value for the vector ‘tspan’. In order to change the step size use the following code tspan=0:0.001:4;
This tells MATLAB that the vector ‘tspan’ is from 0 to 4 with a step size of 0.001. This example uses an arbitrary step size. If the step size is too large, the plot obtained will not be a smooth curve. So it is always better to use relatively smaller step size. But it also takes longer to simulate programs with smaller step size. So you have to decide which is the best step size you can use for a given problem.
8
MEEN 364
Parasuram July 27, 2001
In the above code ‘y(:,1) represents the displacement ‘x’. To plot the velocity, change the variable in the plot command line to ‘y(:,2)’. The plot is attached below
Note: Any nth order differential equation can be solved using MATLAB by following the procedure outlined in example 1. In other words, the nth order differential equation can be reduced to ‘n’ first order differential equations and the ‘ode’ command in MATLAB can be used to solve the system of ‘n’ first order differential equations.
9
MEEN 364
Parasuram July 27, 2001
Example 2. Double Pendulum Derive the governing differential equation of motion of a double pendulum. Solve the derived differential equation and plot the values of θ1 and θ2 with respect to time. Choose l1 = 1m; l2 = 1m; m1 = m2 = 5Kg. The initial conditions for θ1 and θ2 are 0.5233 and 0.5233 radians respectively.
θ1 l1, m1 θ2
l2 , m2
Solution When represented in a matrix form, the differential equations of motion for the above system is . 2 .. m 2 l 2 cos(θ 2 − θ 1 ) θ 1 − ( w1 + w2 ) sin θ 1 + m 2 l 2 θ 2 sin(θ 2 − θ 1) (m 2 + m1 )l1 m l cos(θ − θ ) .. = . m2 l 2 2 1 21 θ 2 − w sin θ − m l θ 1 2 sin(θ − θ 2 2 1 1 2 1)
MATLAB Code These coupled second order differential equations must be converted into a vector of first order differential equation. This step is done by the following substitutions.
θ 1 = y (1); .
θ 1 = y (2); θ 2 = y (3); .
θ 2 = y (4); Substituting the above relations in the original nonlinear differential equation, we get the following first order differential equations, . y (2) 0 0 0 1 y (1) 2 . 0 . ( m + m ) l 0 m l cos( θ − θ ) − ( w + w ) sin θ + m l θ sin( θ − θ ) y ( 2 ) 2 1 1 2 2 2 1 1 2 1 2 2 2 2 1 = . 0 y (3) 0 1 0 y (4) . . 10 2 m2 l 2 0 m2 l1 cos(θ 2 − θ 1 ) 0 y (4) − w2 sin θ 2 − m1l1 θ 1 sin(θ 2 − θ 1 )
MEEN 364
Parasuram July 27, 2001
In this type of a problem where the coefficient matrix is a function of the states or the variables, a separate M-file must be written which incorporates a switch/case programming with a flag case of ‘mass’. For example if the differential equation is of the form, .
M (t , y ) y (t ) = F (t , y ) then the right hand side of the above equation has to be stored in a separate m-file called ‘F.m’. Similarly the coefficient matrix should be stored in a separate m-file named ‘M.m’. So, when the flag is set to the default, the function ‘F.m’ is called and later when the flag is set to ‘mass’ the function ‘M.m’ is called. The code with the switch/case is given below. Note that it is a function file and should be saved as ‘indmot_ode.m’ in the current directory. function varargout=indmot_ode(t,y,flag) switch flag case '' %no input flag varargout{1}=pend(t,y); case 'mass' %flag of mass calls mass.m varargout{1}=mass(t,y); otherwise error(['unknown flag ''' flag '''.']); end
To store the right hand side of the state variable matrix form of the model, a separate function file must be written as shown below. Note that the name of the function is ‘pend’, so this file must be saved as ‘pend.m’. %the following function contains the right hand side of the %differential equation of the form %M(t,y)*y'=F(t,y) %i.e. it contains F(t,y).it is also stored in a separate file named, pend.m. function yp= pend(t,y) M1=5; M2=5; g=9.81; l1=1; l2=1; w2=M2*9.81; w1=M1*9.81; yp=zeros(4,1); yp(1)=y(2); yp(2)=-(w1+w2)*sin(y(1))+M2*l2*(y(4)^2)*sin(y(3)-y(1)); yp(3)=y(4); yp(4)=-w2*sin(y(3))-M2*l1*(y(2)^2)*sin(y(3)-y(1));
11
MEEN 364
Parasuram July 27, 2001
Similarly, to store the coefficient matrix a separate function file is written which is stored as ‘mass.m’. % the following function contains the mass matrix. %it is separately stored in a file named, mass.m function m = mass(t,y) M1=5; M2=5; g=9.81; l1=1; l2=1; m1=[1 0 0 0]; m2=[0 (M1+M2)*l1 0 M2*l2*cos(y(3)-y(1))]; m3=[0 0 1 0]; m4=[0 M2*l1*cos(y(3)-y(1)) 0 M2*l2]; m=[m1;m2;m3;m4];
To plot the response, the main file should call the function ‘indmot_ode.m’, which has the switch/case programming which in turn calls the corresponding functions depending on the value of the flag. For the main file to recognize the coefficient matrix, the MATLAB command ODESET is used to set the mass to ‘M (t, y)’. The MATLAB code for the main file is given below % this is the main file, which calls the equations and solves using ode113. %it then plots the first variable. tspan=[0 10] y0=[0.5233;0;0.5233;0] options=odeset('mass','M(t,y)') [t,y]=ode113('indmot_ode',tspan,y0,options) subplot(2,1,1) plot(t,y(:,1)) grid xlabel('Time') ylabel('Theta1') subplot(2,1,2) plot(t,y(:,3)) grid xlabel('Time') ylabel('Theta2')
For help regarding ‘subplot’ use the online help in MATLAB by typing the following statement at the command prompt. help subplot subplot(m,n,p), breaks the Figure window into an m-by-n matrix of small axes and selects the p-th axes for the current plot. The axes are counted along the top row of the Figure window, then the second row, etc. For example, 12
MEEN 364
Parasuram July 27, 2001
subplot(2,1,1), plot(income) subplot(2,1,2), plot(outgo) plots “income” on the top half of the window and “outgo” on the bottom half. The plot has been attached.
Assignment 1) For the system shown below,
θ 13
MEEN 364
Parasuram July 27, 2001
l
m m = 10 g; l = 5 cms; .
The initial conditions are θ (0) = 90° and θ (0) = 0 a) Derive the governing differential equation. b) Solve the above differential equation obtained analytically and plot the value of theta with respect to time. c) Solve the differential equation using MATLAB and again plot the value of theta with respect to time. Compare the results obtained from (2) and (3). 2) For the system shown below F m
k
c
Using MATLAB, plot the response of a forced system given by the equation ..
.
m x + c x + kx = f sin ω t For ξ = 0.1. .
Take m = 5 kg; k = 1000 N/m; f = 50 N; ω = 4ωn, x(0) = 5 cms; x (0) = 0.
14
MEEN 364
Parasuram July 27, 2001
15
