MATLAB CHAPTER 8 Numerical calculus and differential equations

1
MATLAB ??CHAPTER 8 Numerical calculus and
differential equations

2
Review of integration and differentiation

• Engineering applications
• Acceleration and velocity
• Velocity and distance
• Capacitor voltage and current
• Work expended

3
Integrals

• Properties of integrals
• Definite integrals
• Indefinite integrals

This week's content handles definite integrals
only the answers are always numeric
see chapter 9
4
Improper integrals and singularities

• Improper integrals
• singularity

5
Derivatives

• Integration
  differentiation
• example

product rule
quotient rule
chain rule
1)
2)
3)
6
Numerical integration

• Rectangular integration trapezoidal
  integration
• Numerical integration functions

y
y
(exact solution)
yf(x)
yf(x)


(use of the trapz function)
x
x


a
a
b
b
7
Rectangular, Trapezoidal, Simpson Rules

• Rectangular rule
• Simplest, intuitive
• wi 1
• ErrorO(h) Mostly not enough!
• Trapezoidal rule
• Two point formula, hb-a
• Linear interpolating polynomial
• Simpsons Rule
• Three point formula, h(b-a)/2
• Quadratic interpolating polynomial
• Lucky cancellation due to right-left symmetry
• Significantly superior to Trapezoidal rule
• Problem for large intervals -gt Use the extended
  formula

f0
f1
f(x)
x1b
x0 a
8
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
• adaptive Simpsons rule
• integral of fun from a to b
• tol is the desired error tolerance and is
  optional
• quad1('fun',a,b,tol) Alternative to quad
• adaptive Lobatto integration
• integral of fun from a to b
• tol is the desired error tolerance and is
  optional

9
Comparing the 3 integration methods

• test case
• Trapezoid Method
• xlinspace(0,pi, 10)
• ysin(x)
• trapz(x,y)

(exact solution)

• Simpson's Method
• quad('sin',0,pi)
• Lobatto's Method
• quadl('sin',0,pi)

(
The answer is A1A20.6667, A30.6665
10
Integration near a Singularity

• Trapezoid
• x00.011
• ysqrt(x)
• A1trapz(x,y)
• Simpson
• A2quad('sqrt',0,1)
• Lobatto
• A3quadl('sqrt',0,1)

(The slope has a singularity at x0)
11
Using quad( ) on homemade functions

• 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
• quad('cossq',0,sqrt(2pi))
• 3) The result is 0.6119.

12
Problem 1 and 5 -- Integrate

• Hint position is the integral of velocitydt
• Hint Work is the integral of forcedx

13
Problem 11

• First find V(h4) 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

14
To obtain a vector of integration results

• For example,
• sin(x) is the integral of the cosine(z) from z0
  to zx
• In matlab we can prove this by calculating the
  quad integral in a loop
• for k1101
• 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

15
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

16
Numerical differentiation
backward difference
forward difference
central difference
17
The diff function

• The diff Function
• d diff (x)
• Backward difference central difference method
• example

example x5, 7, 12, -20
d diff(x) d 2 5
-32
18
Try that Compare Backward vs Central

• Central Difference
• d2(y(3n)-y(1n-2))./(x(3n)-x(1n-2))
• subplot(2,1,2)
• plot(x(2n-1),td(2n-1),x(2n-1),d2,'o')
• xlabel('x') ylabel('Derivative')
• axis(0 pi -2 2)
• title('Central Difference Estimate')

• Construct function with noise
• x0pi/50pi
• nlength(x)
• ysin(x) /- .025 random error
• ysin(x) .05(rand(1,51)-0.5)
• tdcos(x) true derivative
• Backward difference using diff
• d1 diff(y)./diff(x)
• subplot(2,1,1)
• plot(x(2n),td(2n),x(2n),d1,'o')
• xlabel('x') ylabel('Derivative')
• axis(0 pi -2 2)
• title('Backward Difference Estimate')

19
Problem 7 Derivative problem

• Hint velocity is the derivative of height

20
Problem 17
21
Polynomial derivatives -- trivia
Example p1 5x 2
p210x24x-3 Result der2 10, 4, -3
prod 150, 80, -7
num 50, 40, 23 den
25, 20, 4

• b polyder (p)
• p a1,a2,,an
• b b1,b2,,bn-1
• b polyder (p1,p2)
• num, den polyder(p2,p1)

Numerical differentiation functions
22
Analytical solutions to differential
equations(1/6)

• Solution by Direct Integration
• Ordinary differential equation (ODE)
• example
• Partial differential equation (PDE) -- not
  covered

23
Analytical solutions to differential
equations(2/6)

• Oscillatory forcing function
• A second-order equation

example
Forcing function
24
Analytical solutions to differential
equations(3/6)

• Substitution method for first-order equations

Natural (by Internal Energy) v.s. Forced Response
(by External Energy) Transient (Dynamics-dependent
) v.s. Steady-state Response (External
Energy-dependent)
25
Problem 18 solve on paper then plot(next week
we solve with Matlab)

• Use the method in previous slide

26
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

27
Analytical solutions to differential
equations(4/6)

• Nonlinear equations
• Substitution method for second order
  equations(1/3)

example
(
)
1)
general solution
suppose that
solution
28
Analytical solutions to differential
equations(5/6)

• Substitution method for second order
  equations(2/3)

2)
(
)
general solution
Eulers identities
period
frequency
3)
(
)
1. Real and distinct s1 and s2.
2. Real and equal s1.
3. Complex conjugates
29
Analytical solutions to differential
equations(6/6)

• Substitution method for second order
  equations(3/3)

• real, distinct roots m1,c8, k15.

characteristic roots s-3,-5.
2. complex roots m1,c10,k601
characteristic roots
3. unstable case, complex roots m1,c-4,k20
characteristic roots
4. unstable case, real roots m1,c3,k-10
characteristic roots s2,-5.
30
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
  solutionsif you don't get the right form, you
  can suspect there is an error somewhere in your
  solution

31
Next Week

• we learn how to solve ODE's with Matlab