MATLAB/Simulink Applications in Solving Ordinary Differential Equations PDF

Preview

Summary

MATLAB/Simulink applications in solving Ordinary Differential Equations Sulaymon Eshkabilov Dynamics & Control Lab Tashkent Automobile Road Institute Str. Mavaraunnaxr 20, Tashkent -700060 Uzbekistan Email: [email protected] Abstract. This paper presents essential points with applications of ...

Description

MATLAB/Simulink applications in solving Ordinary Differential Equations Sulaymon Eshkabilov Dynamics & Control Lab Tashkent Automobile Road Institute Str. Mavaraunnaxr 20, Tashkent -700060 Uzbekistan Email: [email protected] Abstract. This paper presents essential points with applications of MATLAB/Simulink tools in solving initial value problems (IVP) of ordinary differential equations (ODEs) analytically and numerically. There are seven approaches of solving IVP of ODEs in MATLAB/Simulink are discussed via two different examples one which is simulation of thrown ball trajectories and the other is a non-homogeneous second order ODE. Comparative analyses of seven approaches are discussed by plotting their numerical results. Keywords: ODEs, numerical, analytical solutions, MATLAB/Simulink

1.

Introduction

Differential equations (DEs) are very important not only for mathematics or physics but also for broad areas of engineering and many other scientific fields and there are myriad of engineering and scientific problems formulated and simulated using differential equations. For instance, falling or flying object trajectories and equations of motion, radioactive decay, dynamics of processes, logistics problems, wave propagation equations, control problems, just to name a few, are formulated with differential equations. There are many definitions of differential equations and one of the simplest one is “A differential equation is any equation which contains derivatives, either ordinary derivatives or partial derivatives” by (Internet Source: Paul’s Online Math Notes). Despite the fact that modeling the behavior of a system or process or phenomenon with DEs can be relatively easy, but solving many DEs analytically are not plausible and only possibility is then to seek their numerical solutions. There are a number of methods to solve DEs analytically and numerically. For example, to define analytical solution(s) of the DEs separation of variables or introduction of new variables or by expressing them as linear equations or multiplying by some scalars type of approaches (Boyce and DiPrima, 2003), (Internet sources: SOS MATH and wikiHOW) are used depending on the type of the DEs, e.g. homogeneous or non-homogeneous, or linear or non-linear. To find numerical solutions of the DEs, there a number of methods viz. Euler, Runge-Kutta, Runge-Kutta-Gill, Adams, Adams-Bashforth method, Adams-Moulton method, Heun, Ralton, Milne, mid-point rule, Taylor series, and trapezoidal rule methods, developed. Nowadays, to find numerical solutions of differential equations computers and computer programs are mostly employed, and there are many commercial and open ware software packages to find numerical solutions of differential equations. One of such powerful software packages is MATLAB/Simulink that contains many easy to use tools and built-in functions to solve or simulate differential equations. We have studied a few tools and functions of the package in order to show how to employ it in solving initial value problems (IVP) of ordinary differential equations (ODEs). There are two relatively simple examples are taken to demonstrate our developed methodologies and approaches in solving IVP. One of our examples is to simulate thrown ball trajectories and the other is to find analytical and numerical solutions of a second order non-homogeneous ODE. 2. Method 2.1. Example 1. Modeling Thrown Ball Trajectories There are several factors influencing on flight trajectory of a hit or thrown ball. This can be taken as one of the good examples for applications of differential equations for educational purposes. All good tennis players know how to hit the ball to make it fly according to a certain angle to get over the net and at just the right spot of the court on the opponent’s side. Professional hits are the result of a few years of hard practice. If a new player hits the ball 1% faster than a professional tennis player does, it will land considerably far away (maybe about 0.5 m) beyond the baseline. The ball trajectory is a well formulated and explained with the theory of the ball trajectory that is all depicted hereafter.

z

x y Fig. 1. Forces acting on the spinning ball moving in air. From the figure, we can see that there are three forces acting upon the ball flying in air. They are: 1. Gravitational force (1) Where 2.

is vector of the gravitational acceleration. The drag force (2)

That is in opposite direction to the velocity 3.

of the ball.

Magnus force (3)

The magnitudes of the drag force ideal fluids (Richardson, 1961):

and the Magnus force

are usually supposed to have a form given by the theory of

(4) (5) Where is air density. The coefficients and are dependent in real fluids, for instance, air on the velocity , the ball revolution and the ball surface material. These coefficients are found experimentally and taken from the experimental results of the source (Stepanek, 1988), the tennis ball moving with velocity of , the coefficients of projection of the equatorial velocity coefficients were obtained:

and

depend on the ratio of

and with a ball revolution (spin) of only, where

of the spinning ball onto the velocity vector

in some cases the

. The following expressions for the

(6)

(7) For a tennis ball we can neglect the deceleration of the ball revolution, so defined via Newton’s Second Law equations for the position vector

is constant. The trajectory vector of the ball can be

(8) With initial conditions and The equation is a nonlinear system of three differential equations and no analytical solution exists for it, so only way is to solve it numerically. Let us select the x and z axes in x-z plane and then derive the final formulation of the equation (8) by plugging in it and . (9a) (9b) and

Where initial conditions for

The parameter

describes direction of rotation, for instance, for topspin

As

, the following initial values are taken (10)

Where, is the magnitude of the initial velocity vector and is an angle between and x - axis. For our simulations we take two case-conditions which are: 1) ball flying in vacuum and 2) ball flying in air. 1) For the case of ball flying in vacuum, there is one force acting upon the ball is gravitational force and thus, the equations (9a, 9b) will gain a very simple form (11) : (11) 2)

For the case of ball flying in air, in order to find numerical solutions of the problem we rewrite the equations (9a, 9b) of the second order ODEs via a system of the first order ODEs shown in (12) by introducing new variable names.

(12)

. The initial conditions are

Where

For our simulations we use the following numerical values: diameter of the tennis all d

m its mass

and air density

conditions we choose and for the spin of the ball we take

for top spin.

Also, for initial

1)

Case # 1. Ball is hit in the vacuum. This case is very simple and its analytical solution can be easily found with MATLAB Sym olic MATH tool ox’s function dsolve(). Here is a script – MATLAB simulation model.

% Analytical solution of the problem: ball trajectory in the vacuum x_z_vac=dsolve('D2x=0', 'D2z=-g', 'x(0)=0', 'Dx(0)=v0*sin(theta*pi/180)', 'z(0)=h', 'Dz(0)=v0*cos(theta*pi/180)', 't'); x=x_z_vac.x; z=x_z_vac.z; By executing the above script, in the MATLAB command window, the next outputs – analytical solutions (13) of the ball trajectory in vacuum are found. ,

(13)

Here is a complete simulation model implemented via MATLAB (with Sym olic MATH tool ox’s function dsolve( )) script for the case when a tennis ball is hit in vacuum. % Analytical solution of the ball trajectory in vacuum % Case #1. Model_case1.m x_z_vac = dsolve('D2x=0', 'D2z=-g', 'x(0)=0', 'Dx(0)=v0*cos(theta*pi/180)', 'z(0)=h', 'Dz(0)=v0*sin(theta*pi/180)', 't'); xt = x_z_vac.x; zt = x_z_vac.z; g = 9.8; % gravitational acceleration % ball hit 1 m above ground h = 1; theta = 15; % ball hit under 15 degrees v0 = 25; % initial velocity of ball xt = vectorize(xt); zt = vectorize(zt); t = linspace(0,1.5, 200); xt = eval(xt); zt = eval(zt); z_i = find(zt==min(abs(zt))); touch_gr = xt(z_i); t_touch = t(z_i); 2)

Case #2. The numerical simulation models of the tennis ball trajectory are developed in MATLAB and Simulink. Simulations of the ball trajectory in the air under two sub-cases are considered when a ball is hit with topspin and without spin, and also, the numerical simulation model of the case # 1, when a ball is hit in the vacuum, is developed for comparison purposes. The first way, the MATLAB simulation model is created as an M-file with ODE solvers ode23, ode45 and ode113, that is given below.

% Ball Trajectory simulation model with ODE23, ODE45 and ODE113 solvers of MATLAB. h =1; % [m]; v0=25 ; % [m*sec^-1]; theta=15; % [deg]; ro=1.29; % [kg*m^-3]; g=9.81; % [kg*m] d=.063; % [m]; % [kg]; m=.05; w=20; % [m*sec^-1]; angular velocity of a spinning ball eta=1; % describes direction(+-) and presence of rotation; eta=1 is for topspin. alfa = ro*pi*d^2/(8*m); u = @(t,x)sqrt(x(3).^2+x(4).^2); CD = @(t,x)(0.508+(1./(22.053+4.196*(u(t,x)./w).^(5/2))).^(2/5)); CM = @(t,x)(1/(2.022+.981*(u(t,x)./w))); vacuum = @(t,x)([x(3); x(4); 0; -g]); nospin = @(t,x)([x(3); x(4); (-1)*CD(t,x)*alfa*u(t,x)*x(3); (-1)*g-CD(t,x)*alfa*u(t,x)*x(4)]); topspin = @(t,x)([x(3); x(4); (-1)*CD(t,x)*alfa*u(t,x)*x(3)+eta*CM(t,x)*alfa*u(t,x)*x(4); (-1)*g-CD(t,x)*alfa*u(t,x)*x(4)-eta*CM(t,x)*alfa*u(t,x)*x(3)]); time=0:.01:1.5; ICs=[0, 1, v0*cos(theta*pi/180), v0*sin(theta*pi/180)]; [tvac, XZvac] = ode23(vacuum, time, ICs, []); [tns, XZns] = ode45(nospin, time, ICs, []); [tts, XZts] = ode113(topspin, time, ICs, []); % It is also important to find when the ball hits the ground (along x-axis) and after how many seconds in: the vacuum & in the % air with no-spin and top-spin. % Case #1. Vacuum

z_i=find(abs(XZvac(:,2))...