Engineering and Scientific Computation using MATLAB

1
Engineering and Scientific Computation using
MATLAB

• Assoc. Prof. Ph.D. eng. Zagan Remus
• OVIDIUS University of Constanta
• ROMANIA

2
Introduction

• Define the subject matter
• State what the audience will learn in this
  session
• Find out any relevant background and interest of
  the audience

3
OPTIMIZATION AND NUMERICAL SIMULATIONS OF
DIFFERENTIAL EQUATIONS

• List the topics to be covered
• 4.1. Limits of Indeterminate Forms
• 4.2. Derivative a function
• 4.3. Numerical Integration
• 4.4. Numerical Solutions of ODE
• 4.5. MATLAB ODE solvers
• 4.6. Symbolic in MATLAB
• 4.7. Optimization in MATLAB case study

4
4.1. Limits of Indeterminate Forms

• In this section, we discuss a simple algorithm to
  obtain limits of
• using MATLAB. The method consists of the
  following steps
• Construct a sequence of points whose limits is
  x0. In the example below, consider the sequence
• Recall in this regard that as , the
  nth power of any number whose magnitude is
  smaller than one goes to zero.

5

• Construct the sequence of function values
  corresponding to the x-sequence, and find its
  limit.

6
Example 1

• Compute numerically the
• Solution
• Enter the following instructions in your MATLAB
  command window
• N20 n1N
• x00
• dxn-(1/2).n
• xnx0dxn
• ynsin(xn)./xn
• plot(xn,yn)

7

• Result of plot command is presented in figure
  bellow

8

• The limit of the yn sequence is clearly equal to
  1. The deviation of the sequence of the yn from
  the value of the limit can be obtained by
  entering
• dynyn-1
• semilogy(n,dyn)
• The last command plots the curve with the
  ordinate y expressed logarithmically (see fig 2).
• This mode of display is the most convenient in
  this case because the ordinate spans many decades
  of values.

9

• Figure 2

10
4.2 Derivative of a Function

• The derivative of certain function at a
  particular point is defined as

11

• Numerically, the derivative is computed at the
  point x0 as follows
• Construct an x-sequence that approaches x0.
• Compute a sequence of the function values
  corresponding to the x-sequence.
• Evaluate the sequence of the ratio, appearing in
  the definition of derivative.
• Read off the limit of this ratio sequence. This
  will be the value of derivative at the point x0.

12
Example 2

• Find numerically the derivative of the function
• ln(1 x) at x 0

13

• Solution
• Edit and execute the following script M-file
• N20n1N
• x00
• dxn(1/2).1N
• xnx0dxn
• ynlog(1xn)
• dynyn-log(1x0)
• deryndyn./dxn
• plot(n,deryn)
• The result of plot command is presented in figure
  3.

14

• Figure 3

15

• The limit of the deryns sequence is clearly
  equal to 1
• and
• the value of this function derivative at 0.

16
Example 3

• Plot the derivative of the function
• x2 sin(x)
• over the interval 0 x 2 p

17

• Solution
• Edit and execute the following script M-file
• dx10(-4)
• x0dx2pidx
• dfdiff(sin(x).x.2)/dx
• plot(0dx2pi,df)

18

• where diff is a MATLAB command, which when acting
  on an array X, gives the new array
• X(2) X(1), X(3) X(2) ,, X(n) X(n 1),
• whose length is one unit shorter than the array
  X.
• The accuracy of the above algorithm depends on
  the choice of dx. Ideally, the smaller it is, the
  more accurate the result.
• The result of plot command is presented in figure
  4

19

• Figure 4

20
4.3 Numerical Integration

• The algorithm for integration discussed in this
  section is the second simplest available (the
  trapezoid rule being the simplest, beyond the
  trivial, is given at the end of this section as a
  problem).
• It has been generalized to become more accurate
  and efficient through other approximations,
  including Simpsons rule, the Newton-Cotes rule,
  the Gaussian-Laguerre rule, etc.

21

• Here, we perform numerical integration through
  the means of a Riemann sum we subdivide the
  interval of integration into many subintervals.
• Then we take the area of each strip to be the
  value of the function at the midpoint of the
  subinterval multiplied by the length of the
  subinterval, and we add the strip areas to obtain
  the value of the integral.
• This technique is referred to as the midpoint
  rule.

22

• We can justify the above algorithm by recalling
  the Mean Value Theorem of Calculus, which states
  that
• where c is between a, b. Thus, if we divide the
  interval of integration into narrow subintervals,
  then the total integral can be written as the sum
  of the integrals over the subintervals, and we
  approximate the location of c in a particular
  subinterval by the midpoint between its
  boundaries .

23
Example 4

• Use the above algorithm to compute the value of
  the definite integral of the function
• sin(x)
• from 0 to p.

24

• Solution
• Edit and execute the following program
• dxpi/200
• x0dxpi-dx
• xshiftxdx/2
• yshiftsin(xshift)
• Intdxsum(yshift)
• Give the answer Int 2.0000

25
4.4 Numerical Solutions of Ordinary Differential
Equations

• Ordinary linear differential equations are of the
  form

26

• The as are called the coefficients and u(t) is
  called the source (or input) term.
• Ordinary differential equations (ODEs) show up in
  many problems of electrical engineering,
  particularly in circuit problems where, depending
  on the circuit element, the potential across it
  may depend on the deposited charge, the current
  (which is the time derivative of the charge), or
  the derivative of the current (i.e., the second
  time derivative of the charge) that is, in the
  same equation, we may have a function and its
  first- and second-order derivatives

27
Example 5

• Solve, over the interval 0 lt t lt 1, the following
  second order differential equation
• with the initial conditions
• y(t0)3/8
• and
• y(t0)0

28

• First we make some notation
• And the system before become

29

• Also is useful to consider that
• and the improved expression for the numerical
  differentiator, D(k)

30

• Substituting the above first-order expression of
  the iterators for the first-order and
  second-order numerical differentiators, we deduce
  the following iterative equation for y(k)

31

• With the notation maked before, now we are ready
  to give the solution of the our problem in MATLAB
  instruction

32

• Solution Edit and execute the following script
  M-file
• tin0
• tfin1
• tlinspace(tin,tfin,2000)
• Nlength(t)
• a1-t.2
• b-2t
• c20ones(1,N)
• yzeros(1,N)
• Dzeros(1,N)
• dt(tfin-tin)/(N-1)
• uzeros(1,N)
• y(1)3/8
• D(1)0
• D2(1)(1/a(1))(-b(1)D(1)-c(1)y(1)u(1))

33

• for k2N
• y(k)((4a(k)/dt22b(k)/dtc(k))(-1))...
• (y(k-1)(4a(k)/dt22b(k)/dt)D(k-1)...
• (4a(k)/dtb(k))a(k)D2(k-1)u(k))
• D(k)(2/dt)(y(k)-y(k-1))-D(k-1)
• D2(k)(4/dt2)(y(k)-y(k-1))-(4/dt)D(k-1)-...
• D2(k-1)
• end
• yanal(35t.4-30t.23)/8
• plot(t,y,t,yanal,'--')

34

• As you will observe upon running this program,
  the numerical solution and the analytical
  solution agree very well.
• The plot result is in figure 5

35

• Figure 5

36
4.5 Higher-Order Iterators The Runge-Kutta
Method

• In this subsection, we want to explore the
  possibility that if we sampled the function
  n-times per step, we will obtain a more accurate
  solution to the ODE (Ordinary Differential
  Equation) than that obtained from the .rst-order
  iterator for the same value of t.

37

• To focus the discussion, consider the ODE
• y'(t) f(t, y(t))
• Higher-order ODEs can be reduced, to a system of
  equations having the same functional form as
  previous equation.
• The derivation of a technique using higher-order
  iterators will be shown below in detail for two
  evaluations per step. Higher-order recipes can be
  found in most books on numerical methods for ODE.

38
4.6. MATLAB ODE Solvers

• MATLAB has many ODE solvers, ODE23 and ODE45
  being most commonly used. ODE23 is based on a
  pair of second-order and third-order Runge-Kutta
  methods running simultaneously to solve the ODE.
  The program automatically corrects for the step
  size if the answers from the two methods have a
  discrepancy at any stage of the calculation that
  will lead to a larger error than the allowed
  tolerance.

39

• To use this solver, we start by creating a
  function M-file that includes the system of
  equations under consideration.
• This function is then called from the command
  window with the ODE23 or ODE45 command.
• To compute the dynamics of the problem, we
  proceed into two steps.
• First, we generate the function M-file that
  contains the rate equations, and then proceed to
  solve these ODEs by calling the MATLAB ODE
  solver.

40
Example 6

• Specifically the function M-file representing the
  laser rate equations is given by
• function yplaser1(t,y)
• p30 pumping rate
• gamma10(-2) inverse natural lifetime
• B3 stimulated emission coefficient
• c30 inverse lifetime of photon in cavity
• yp(1,1)p-gammay(1,1)-By(1,1)y(2,1)
• yp(2,1)-cy(2,1)By(1,1)y(2,1)

41

• The script M-file to compute the laser dynamics
  and thus simulate the relaxation oscillations is
• tspan0 3
• yin1 1
• t,yode23('laser1',tspan,yin)
• subplot(3,1,1)
• plot(t,y(,1))
• xlabel('Normalized Time')
• ylabel('N')
• subplot(3,1,2)
• plot(t,y(,2))
• xlabel('Normalized Time')
• ylabel('n')
• subplot(3,1,3)
• plot(y(,1),y(,2))
• xlabel('N')
• ylabel('n')

42

• As can be observed in figure 6, the oscillations,
  as predicted, damp-out after a while and the
  dynamical variables reach a steady state. The
  phase diagram, shown in the bottom panel, is an
  alternate method to show how the population of
  the atomic higher excited state and the photon
  number density reach the steady state

43

• Figure 6

44

• A simplest form of the syntax for the MATLAB ODE
  solvers is
• t, y solver(fun, tspan, y0,
• where
• fun is a string containing name of the ODE m-file
  that describes the differential equation,
• tspan is the interval of integration, and y0 is
  the vector holding the initial value(s).

45
Example 7

• In the following example we will seek a numerical
  solution y
• at t 0, .25, .5, .75, 1 to the following
  initial value problem
• y' -2ty2,
• with the initial condition y(0) 1.
• We will use both the ode23 and the ode45 solvers.
  
• The exact solution to this problem is
• y(t) 1/(1 t2)

46

• The ODE m-file needed in these computations is
  named eq1.
• function dy eq1(t,y)
• The m-file for the ODE y' -2ty2.
• dy -2t.y(1).2

47

• then type MATLAB statement in Window Command
• format long
• tspan 0 .25 .5 .75 1 y0 1
• t1 y1 ode23('eq1', tspan, y0)
• t2 y2 ode45('eq1', tspan, y0)

48

• To compare obtained results let us create a
  three-column table holding the points of
  evaluation and the y-values obtained with the aid
  of the ode23 and the ode45 solvers.

49

• t1 y1 y2
• ans
• 0 1.00000000000000 1.00000000000000
• 0.25000000000000 0.94118221525751
  0.94117646765650
• 0.50000000000000 0.80002280597122
  0.79999999678380
• 0.75000000000000 0.64001788410487
  0.63999998775736
• 1.00000000000000 0.49999658522366
  0.50000000471194

50
Example 8

• Instead of writing the ODE m file for this
  system,
• y1'(t) y1(t) 4y2(t), y2'(t) -y1(t) y2(t),
• y1(0) 1 y2(0) 0.

51

• we will use MATLAB inline function
• dy inline('1 4-1 1y', 't', 'y')
• give this answer

52

• dy
• Inline function
• dy(t,y) 1 -4-1 1y

53

• The inline functions are created in the Command
  Window.
• Interval over which numerical solution is
• computed and the initial values are stored in the
  vectors tspan and y0, respectively
• tspan 0 1 y0 1 0
• Numerical solution to this system is obtained
  using the ode23 function
• t,y ode23(dy, tspan, y0)

54

• Give the answer
• t
• 0
• 0.00008000000000
• 0.00048000000000
• 0.00248000000000
• 0.01248000000000
• 0.04746845566255
• 0.10101102368152
• 0.16785081419689
• 0.24416687626714
• 0.32719104915831
• 0.41492011412819
• 0.50591928384684
• 0.59917400535493
• 0.69397526179837
• 0.78983224264696
• 0.88640761940849
• 0.98347075679869

55

• y
• 1.00000000000000 0
• 1.00008001600111 -0.00008000640060
• 1.00048057623966 -0.00048023052905
• 1.00249540908542 -0.00248616821334
• 1.01287361189987 -0.01263803003396
• 1.05333952580664 -0.04984951641514
• 1.12892127338309 -0.11249929331960
• 1.24996221189645 -0.20224167583097
• 1.43162623170992 -0.32413576345484
• 1.69437867683002 -0.48671811426647
• 2.06541997595762 -0.70251596876793
• 2.58095589060030 -0.98900563393263
• 3.28935422493395 -1.37004839736238
• 4.25539578621716 -1.87790924043463
• 5.56596129293983 -2.55602504924514
• 7.33763883947222 -3.46275783101774
• 9.72691149772029 -4.67645587101605

56

• Graphs of y1(t) (solid line) and y2(t) (dashed
  line) are shown below

57

• And the MATLAB statement are
• plot(t,y(,1),t,y(,2),'--'), legend('y1','y2'),
  xlabel('t'),
• ylabel('y(t)'),
• title('Numerical solutions y_1(t) and y_2(t)')

58

• The exact solution (y1(t), y2(t)) to this system
  is
• y1, y2
• y1
• 1/2exp(-t)1/2exp(3t)
• y2
• -1/4exp(3t)1/4exp(-t)
• Functions y1 and y2 were found using command
  dsolve which is available in the Symbolic Math
• Toolbox.

59
Example 9

• The following set of two nonlinear differential
  equations, called the van der Pol equations,

60

• Has been used as an illustrative example to solve
  ordinary differential equations using different
  solvers over last 18 years.
• Two m-files to solve these differential equations
  are given in MATLAB software, namely
• vdpode.m and vdp1000.m
• Both files are in the particular MATLAB
  directory.
• To perform numerical simulations, we run file
  vdpode.m by typing in the Command Window
• vdpode

61

• After running the results is presenting in figure
  below

62
Symbolic using MATLAB
63
Symbolic using MATLAB
64
Symbolic using MATLAB
65
Symbolic using MATLAB
66
Symbolic using MATLAB
67
Symbolic using MATLAB
x 2, 3
x 1, 3 y 1, 3/2
4/3
68
ODE solution using Symbolic MATLAB
dsolve('Dy3y','x')
C1exp(3x)
69
ODE solution using Symbolic MATLAB
dsolve('D2y-3Dy2y4exp(3x)','x')
2exp(3x)C1exp(x)C2exp(2x)
70
ODE solution using Symbolic MATLAB
dsolve('2xy(x2cos(y))Dy0','x')
C1x2ysin(y)0
71
Optimization using MATLAB

• A company must produce 5 types of products namely
  R1, R2, R3, R4, R5 using 3 ranges of raw
  materials M1, M2, M3. Consumption for these raw
  materials to make the products, the unitary
  profits and the available funds for each raw
  materials for the entirely period analyzed are
  given in the table below

72
Optimization using MATLAB
73
Optimization using MATLAB

• Knowing that the society have the capacity for
  any quantity of products that will be produced
  and all the products will be sell on the market,
  write a linear model for determining optimum
  quantity that will be produced for all these 5
  products, as the total profit for selling to be
  maximum and solve these mathematical model using
  MATLAB.

74
Optimization using MATLAB

• Solving We noted with xj , j 1, , 5 the
  quantity of products that will be produced Rj .
  The mathematical model of the problem contain
• - restrictions regarding the availability of
  limited raw materials

75
Optimization using MATLAB

• - restrictions of non-negativity imposed to
  variables of model given by financial contain

76
Optimization using MATLAB

• - objective function to be optimize

77
Optimization using MATLAB

• The mathematical model make, in Matlab consist
  by
• A0 1.2 1 1 0.5
• 0 1 0.5 1 1
• 1 1.3 0 2 1.5
• b1000 500 800
• x00 0 0 0 0
• f-50 -300 -30 -60 -40
• lbzeros(5,1)
• x,profitlinprog(f,A,b,,,lb)

78
Optimization using MATLAB

• after running we obtain
• x
• 150.0000
• 500.0000
• 0.0000
• 0.0000
• 0.0000
• profit 1.5750e005
• That, the optimum quantity that will be produced
  by society to obtain the maximum benefit of 157
  500 a.u., are
• for product 1 , 150 pieces
• for product 2, 500 pieces

79
Where to Get More Information

• Other training sessions
• List books, articles, electronic sources
• Consulting services, other sources

80

• Thanks for listening.
• Have a nice day.