# MATLAB ch08a

Document Sample

```					 MATLAB 입문
CHAPTER 8
Numerical calculus and differential equations

ACSL, POSTECH                            1
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations                                            Acsl, Postech

Review of integration and differentiation

     Engineering applications
f(x)

   Acceleration and velocity:
A
b
v(b)   a (t )dt  v(0)
0

   Velocity and distance:
b
x(b)   v(t )  x(a)
a

x
   Capacitor voltage and current:
a                    b
1 b
i(t )dt  Q(a)
C a
v(b) 
                 


b                                                               Work expended:
A   f ( x)dx                                                                      d
a
W   f ( x)dx
0
d
W   kxdx
0

2
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations                                    Acsl, Postech

Integrals

   Properties of integrals
b                              b                    b
 [cf ( x)  dg ( x)]dx  c 
a                               a
f ( x)dx  d  g ( x)dx
a

b               c                  b
        f ( x)dx   f ( x)dx   f ( x)dx
a               a                c
This week's content
   Definite integrals                                                                               handles definite
x n1 x b b n1 a n1
  a
b
x dx 
n
     
n  1 x a n  1 n  1
n  1             integrals only
b    1           x b
a       x
dx  ln x
x a
 ln b  ln a
2
           sin xdx   cos x
x  2
 [cos2  cos ]  2      always numeric
x 

   Indefinite integrals

 cos xdx  sin x                                                          see chapter 9
3
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations   Acsl, Postech

Improper integrals and singularities

   Improper integrals

1
 x 1dx  ln | x 1 |
h    1
lim             dx  lim ln x  1 h
h    1-
0    x 1     h 1-
0

 lim(ln | h  1 |  ln | 1 |)
h   1-

 lim ln(1  h)  
h   1-

   singularity

4
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations                                                 Acsl, Postech

Derivatives
     Integration                                           differentiation
dg
g ( x)   f ( x)dx                                   f ( x) 
dx

dh          dg           df
h( x )  f ( x ) g ( x )                                  f ( x)      g ( x)               :    product rule
dx          dx           dx
df           dg
g ( x)  f ( x)
dh          dx           dx
h( x)  f ( x) / g ( x)                                               2
:    quotient rule
dx             g

dz dz dy
z  f ( y ) y  g (x)                                                                        :   chain rule
dx dy dx

example

dx n                                                             d ( x 2 sin x)
1)       nx n 1                                            2)                     x 2 cos x  2 x sin x
dx                                                                     dx
d ln x 1

dx      x
d sin x
 cos x                                                  d (sin 2 x)      dy
3)                 2y     2 sin x cos x
dx                                                                 dx           dx
d cos x
  sin x
dx
5
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations                                                 Acsl, Postech

Numerical integration
Rectangular integration                     trapezoidal integration

y                                           y
0
sin dx


   0
sin dx   cos x   cos 0  cos   2 (exact solution)
0
y=f(x)                                      y=f(x)

···                                          ···
(use of the trapz function)
x                                           x
a         ···         b                     a          ···       b

Numerical integration functions
Command                              Description
quad (‘function’,a,b,tol)            Uses an adaptive Simpson’s rule to compute the integral of the function ‘function’ with
a as the lower integration limit and b as the upper limit. The parameter tol is optional.
Tol indicates the specified error tolerance.
quadl (‘function’, a,b,tol)          Uses Lobatto quadrature to compute the integral of the function ‘function’. The rest of
the syntax is identical to quad.
trapz (x,y)                          Uses trapezoidal integration to compute the integral of y with respect to x, where the
array y contains the function values at the points contained in the array x.
6
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations                         Acsl, Postech

Rectangular, Trapezoidal, Simpson Rules
f(x2)
N 1                    f(x1)
I  h f i
    Rectangular rule                                                        f(x)
 Simplest, intuitive
 wi = 1
i 0

 Error=O(h)     Mostly not enough!
    Trapezoidal rule                                                               x0 =a x1 x2         xN=b
 Two point formula, h=b-a
 Linear interpolating polynomial
f(x)       f0              f1
I  h( f 0  f1 ) / 2  O(h 3 )
    Simpson’s Rule
 Three point formula, h=(b-a)/2
x0 =a                x1=b
 Lucky cancellation due to right-left symmetry
1        4        1
I  h( f 0  f1  f 2 )  O(h5 )

3        3        3
Significantly superior to Trapezoidal rule
    Problem for large intervals -> Use the extended formula

7
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations           Acsl, Postech

Matlab's Numerical Integration Functions
   trapz(x,y)                             When you have a vector of data
    trapezoidal integration method
    not as accurate as Simpson's method
    very simple to use

   quad('fun',a,b,tol)                                     When you have a function
    integral of ’fun’ from a to b
    tol is the desired error tolerance and is optional

    integral of ’fun’ from a to b
    tol is the desired error tolerance and is optional                                   8
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations                    Acsl, Postech

Comparing the 3 integration methods

test case:                      0
sin dx   cos x   cos 0  cos   2
0
(exact solution)
Trapezoid Method:
x=linspace(0,pi, 10);                                          Simpson's Method:
(

trapz(x,y);
Lobatto's Method

9
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations                           Acsl, Postech

Integration near a Singularity
Trapezoid:
x=[0:0.01:1];
y=sqrt(x);
A1=trapz(x,y);
y x
dy / dx  0.5 / x
Simpson:                                                                  (The slope has a singularity at x=0)

Lobatto:

10
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations   Acsl, Postech

1) Create a function file:
function c2 = cossq(x)
% cosine squared function.
c2 = cos(x.^2);

Note that we must use array exponentiation.

2) The quad function is called as follows:

3) The result is 0.6119.
11
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations   Acsl, Postech

Problem 1 and 5 -- Integrate

   Hint: position is the integral of velocity*dt
   Hint: Work is the integral of force*dx

12
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations   Acsl, Postech

Problem 11

   First find V(h=4) which is the Volume of the full cup
   Hint: Flow = dV/dt, so the integral of flow (dV/dt) from 0 to t = V(t)
   Set up the integral using trapezoid (for part a) and Simpson (for b)
   Hint: see next slide for how to make a vector of integrals with changing b values
   So you can calculate vector V(t), i.e. volume over time for a given flow rate
    V = quad(....., 0, t) where t is a vector of times
   and then plot and determine graphically when V crosses full cup                  13
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations                          Acsl, Postech

To obtain a vector of integration results:

   For example,
   sin(x) is the integral of the cosine(z) from z=0 to z=x
   In matlab we can prove this by calculating the quad integral in a loop:
for k=1:101
x(k)= (k-1)* pi/100;                             % x vector will go from 0 to pi
sine(k)=quad('cos',0,x(k));                      % this calculates the integral from 0 to x
end
plot(x, sine)                                        % this shows the first half of a sine wave
% calculated by integrating the cosine

14
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations   Acsl, Postech

Problem 14

   Hint: The quad function can take a vector of values for the endpoint
   Set up a function file for current
   Set up a vector for the time range from 0 to .3 seconds
   Then find v = quad('current',0,t)     will give a vector result
   and plot v
15
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations                              Acsl, Postech

Numerical differentiation

dy       y
 lim
True slope                         dx Δx 0 x
y3
y=f(x)
y2

y 2  y1 y 2  y1
B

A
C
mA                          : backward difference
y1
x 2  x1   x
Δx            Δx                            y3  y 2 y3  y2
mB                          : forward difference
x1            x2            x3
x3  x 2   x
mA  mB 1  y 2  y1 y 3  y 2  y 3  y1
mC                                         : central difference
2    2  x          x       2x

16
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations                        Acsl, Postech

The diff function

   The diff Function                                           example    x=[5, 7, 12, -20];
d= diff(x)
d= diff (x)                                                            d=2      5   -32

   Backward difference & central
difference method
example

17
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations                              Acsl, Postech

Try that: Compare Backward vs Central

  Construct function with noise                                              Central Difference
x=[0:pi/50:pi];                                                            d2=(y(3:n)-y(1:n-2))./(x(3:n)-x(1:n-2));
n=length(x);                                                               subplot(2,1,2)
% y=sin(x) +/- .025 random error                                           plot(x(2:n-1),td(2:n-1),x(2:n-1),d2,'o');
y=sin(x) + .05*(rand(1,51)-0.5);                                           xlabel('x'); ylabel('Derivative');
td=cos(x); % true derivative                                               axis([0 pi -2 2])
title('Central Difference Estimate');
   Backward difference using diff:
d1= diff(y)./diff(x);
subplot(2,1,1)
plot(x(2:n),td(2:n),x(2:n),d1,'o');
xlabel('x'); ylabel('Derivative');
axis([0 pi -2 2])
title('Backward Difference Estimate')                                                                           18
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations   Acsl, Postech

Problem 7: Derivative problem

   Hint: velocity is the derivative of height

19
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations   Acsl, Postech

Problem 17

20
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations                                                                       Acsl, Postech

Polynomial derivatives -- trivia
f ( x)  a1 x n  a 2 x n 1  a3 x n  2    an  1 x  an                             Example p1= 5x +2
df
 na1x n1  (n  1)a 2 x n2    an  1                                                            p2=10x2+4x-3
dx
 b1x n1  b2 x n2    bn  1

       b = polyder (p)
p = [a1,a2,…,an]
b = [b1,b2,…,bn-1]                                                                    Result:       der2 = [10, 4, -3]
prod = [150, 80, -7]
       b = polyder (p1,p2)
num = [50, 40, 23]
       [num, den] = polyder(p2,p1)                                                                         den = [25, 20, 4]

Numerical differentiation functions
Command                          Description

d=diff(x)                        Returns a vector d containing the differences between adjacent elements in the vector x.

b=polyder(p)                     Returns a vector b containing the coefficients of the derivative of the polynomial represented by the vector p.

b=polyder(p1,p2)                 Returns a vector b containing the coefficients of the polynomial that is the derivative of the product of the polynomials
represented by p1 and p2.       cf. equivalent comment: b=polyder(conv(p1,p2))

[num, den]=                      Returns the vector num and den containing the coefficients of the numerator and denominator polynomials of the
polyder(p2,p1)                   derivative of the quotient p2/p1, where p1 and p2 are polynomials.

21
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations   Acsl, Postech

Analytical solutions to differential equations(1/6)

   Solution by Direct Integration
   Ordinary differential equation (ODE)

example
dy 2
t
dt
t dy       t        t3 t t3
0 dt dt  0 t dt  3 0  3
2

t3
y (t )  y (0) 
3
dy
y (t ) 

dt
d2y
(t )  2
y
dt

   Partial differential equation (PDE) -- not covered

22
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations                                                 Acsl, Postech

Analytical solutions to differential equations(2/6)

   Oscillatory forcing function
dy
 f (t )        Forcing function                    example      f (t )  sin wt
dt
dy
t         t              cos wt t 1  cos wt
0 dt
dt   sin wtdt  
0                 w 0

w
t                  1  cos wt
y(t )  y(t )  y(0) 
0                       w
1  cos wt
y (t )  y (0) 
w
   A second-order equation
d2y 3
t
dt 2
t d2y       dy           t        t4
0 dt 2 dt  dt  y(0)  0 t dt  4
3


dy t 4
  y (0)

dt 4
t t     
4
t dy                                     t5
0 dt dt  y(t )  y(0)     y(0)dt 

0 4

     20
 ty(0)


y (t )  t 5 / 20  ty (0)  y (0)

23
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations                                                     Acsl, Postech

Analytical solutions to differential equations(3/6)
     Substitution method for first-order                                            suppose
equations                                                                      f (t )  0   for    t 0

dy                                                                                  M at        t=0
       y  f (t ) ( y (t )  Ae st , dy / dt  sAe st    )   f (t )  0
dt                                                                                Such a function is the step function

dy                                                                              y (t )  Ae st  B
  y  sAest  Aest  0
dt                                                                              B  y(0)  A
s  1  0               Characteristic equation                                y (t )  Ae st  y (0)  A
s  1/                 Characteristic root                                   ( s  1  0 ,     A  y(0)  M )
y (0)  Ae 0  A                 to find A                                      y (t )  y (0)e t /  M (1  e t / )   Forced response

y (t )  y (0)e  t /          Solution, Free response

∵The solution y(t) decays with time if τ>0
τ is called the time constant.

Natural (by Internal Energy) v.s. Forced Response (by External Energy)
Transient (Dynamics-dependent) v.s. Steady-state Response (External Energy-dependent)

24
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations   Acsl, Postech

Problem 18 – solve on paper then plot
(next week we solve with Matlab)
   Use the method in previous slide

25
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations   Acsl, Postech

Problem 21 – solve on paper then plot
(next week we solve with Matlab)

  Re-arrange terms:
mv' + cv = f , OR
(m/c) v' + v = f/c

now it's in the same form as 2 slides back
26
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations   Acsl, Postech

Analytical solutions to differential equations(4/6)

        Nonlinear equations
example
y  5 y  y  0
y      
y  sin y  0

y y 0


        Substitution method for second order equations(1/3)

1)     c 2 y  0
y                  (   y (t )  Ae st   )

( s 2  c 2 ) Ae st  0
s  c
y (t )  A1e ct  A2 e  ct         general solution

suppose that       c  2, y(0)  6, y(0)  4

y(0)  6  A1  A2
y (0)  2 A1  2 A2  4

A1  4, A2  2             solution

27
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations      Acsl, Postech

Analytical solutions to differential equations(5/6)

        Substitution method for second order equations(2/3)
2)      2 y  0
y                      (   y (t )  Ae st      )
it             it
y (t )  A1e          A2e                      general solution

Euler’s identities            eit  cost  i sin t
y(t )  B1 sin t  B2 cos t
B1  y(0) / , B2  y(0)


y(0)
y(t )              sin t  y(0) cost

2
period     P                   frequency            f  1/ P

3)   m  cy  ky  f (t )
y                                   (   f (t )  0, y (t )  Ae st   )

(ms 2  cs  k ) Ae st  0
1. Real and distinct: s1 and s2.
y (t )  A1e s1t  A2 e s2t
2. Real and equal: s1.
y (t )  ( A1  tA2 )e s1t
3. Complex conjugates:  s  a  i
y (t )  B1e sin t  B2 e at cos t
at

28
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations   Acsl, Postech

Analytical solutions to differential equations(6/6)

        Substitution method for second order equations(3/3)
1.    real, distinct roots: m=1,c=8, k=15.
characteristic roots s=-3,-5.

y (t )  A1e 3t  A2 e 5t

2.    complex roots: m=1,c=10,k=601
characteristic roots   s  5  24i
y (t )  B1e 5t sin 24 t  B2 e 5t cos 24 t
P  2 / 24   / 12

3.    unstable case, complex roots: m=1,c=-4,k=20
characteristic roots   s  2  4i
y (t )  B1e 2t sin 4t  B2e 2t cos 4t
P  2 / 4   / 2
4.    unstable case, real roots: m=1,c=3,k=-10
characteristic roots s=2,-5.
y (t )  A1e 2t  A2e 5t

29
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations   Acsl, Postech

Problem 22

    Just find the roots to the characteristic equation and determine which
form the free response will take (either 1, 2, 3 or 4 in previous slide)

    This alone can be a helpful check on matlab's solutions—if you don't
get the right form, you can suspect there is an error somewhere in
MATLAB 입문 : Chapter 8. Numerical calculus and differential equations   Acsl, Postech

Next Week:

   we learn how to solve ODE's
with Matlab

31

```
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
 views: 10 posted: 10/4/2012 language: Latin pages: 31
How are you planning on using Docstoc?