Loading...

PPT – Numerical Methods PowerPoint presentation | free to download - id: 4906d6-MGY2Z

The Adobe Flash plugin is needed to view this content

Numerical Methods

- Numerical Methods are techniques by with

mathematical problems are

formulated so that they can be solved on

computers - Powerful analytical methods doesnt yield easily

to non-linear, complex geometrical and high

dimension problems - Numerical methods are extremely powerful

problem-solving tools. They are capable of

handling large systems of equation,

nonlinearities and complicated geometries

impossible to solve analytically. - MATLAB and its toolbox contains an extensive

library for solving many practical numerical

problems e.g. root-finding, interpolation,

numerical calculus, solving systems of linear and

non linear equation , Ordinary Differential

equation.

Numerical Methods

- Contents
- Solving linear equation of system
- Matrix Factorization/Decomposition Techniques
- Curve fitting and polynomials operations
- Interpolation
- Numerical Integration and Differentiation
- Ordinary Differential Equation
- Non-linear equations

Linear Algebra

- Solving a linear system
- 5x 3y - 2z 10
- 8y 4z 3x 20
- 2x 4y - 9z 9
- Cast in standard matrix form
- A x b

Linear Equations

- A 5 -3 2 -3 8 4 2 4 -9
- b 10 20 9
- Left Division Method
- x A \ b
- Check solution
- c Ax
- c should be equal to b.

Gaussian Elimination

- Pivoting
- C A b augmented matrix
- row reduced echelon form
- Cr rref(C)
- Cr
- 1.0000 0 0 3.4442
- 0 1.0000 0 3.1982
- 0 0 1.0000 1.1868

Eigenvalues and Eigenvectors

- The problem A v ? v
- With pencil-paper, we must solve for
- det (A - ?) 0
- then solve the linear equation
- MATLAB way
- V, D eig(A)
- V,D EIG(X) produces a diagonal matrix D of

eigenvalues and a - full matrix V whose columns are the corresponding

eigenvectors so that XV VD.

Eigenvectors Eigenvalues

Interpretation of Eigenvectors Eigenvalues (1)

Matrix Factorizations

- LU decomposition
- L, U lu(A)
- such that LU A, L is a lower triangular

matrix, - U is an upper triangular matrix.
- DSP blockset (Simulink)

Matrix Factorization

- QR factorization Application in Curve fitting,
- Q, R qr(A) mgtn over determined
- If nrank(A) then A?Q,R ? Anxn
- such that QR A Q is orthogonal matrix

,independent row/coloums and R is upper

triangular matrix. QTQ-1, QTQI - DSP blockset (Simulink)

Matrix factorization

- Singular Value Decomposition
- U, D, V svd(A)
- such that UDV A, where U and V are orthogonal

and D is diagonal matrix Eigen values. - DSP Blockset
- A 1 1 1 2 -1 2 -1 4 3 4 2 1 3 -3 4
- b 1 2 -1 4 8
- S svd(A)
- U,S,V svd(A)
- USV'

Sparse Matrix

- See
- help sparfun
- for detail
- The sparse function generates matrices in the

MATLAB sparse storage organization. - S sparse(A) converts a full matrix to sparse

form by squeezing out any zero elements. If S is

already sparse, sparse(S) returns S. - A1 0 1 0 2 0 3 3 0
- Ssparse(A)
- is 1 3 2 3 1,js 1 1 2 2 3,as 1 3 2 3 1
- Ssparse(is,js,as) Tabular..

Curve Fitting Polyfit Polyval-

- Polynomial curve fitting
- apolyfit(x,y,n)
- do a least-square fit with polynomial of degree

n. x and y are data vectors, and a is

coefficients of the polynomial, a an, an-1, ,

a1, a0 - Store data and perform the fit
- x 1 2 4 5 data
- y 1 2 2 3 data
- c polyfit(x,y,n) n1 degree
- Evaluate and plot the fit with
- xfit0 6
- yfit c(1)xfit c(2) yc(1)x c(2)
- plot(x,y,o ,xfit,yfit,-)
- grid

Curve Fitting- Non Linear Functions

- y polyval(a, x)
- compute the value of polynomial at value x
- y a(1) xn a(2) x(n-1) a(n) x

a(n1) - The Polynomial is evaluated

at x5, 7, and 9 with - p 3 2 1
- polyval(p,5, 7, 9)
- CURVE FITTING TOOLBOX

Curve Fitting- Least Square Fit

Normal Equations Consider the problem Ax?b, where

A is an (m,n) matrix with m?n , rank(A) n, b is

an (m,1) vector, and x is the (n,1) vector, to be

over determined. More equation than unknowns. In

such case exact solution is not obtained

Axb b ?c Least

square provides the solution such that ?

Normal Equation Least Square curve Fit

Fit x,y data to x 0.955 1.380 1.854 2.093

2.674 3.006 3.225 3.940 4.060 y5.722 4.812

4.727 4.850 5.011 5.253 5.617 6.2282

6.255 A1./x x Coefficients matrix of

overdetermined system c(AA)\(Ay) Solve

Normal equation xf linspace(min(x),

max(x)) yfc(1)./xf c(2)xf Af1./xf

xf yfAfc plot(x, y,o,xf, yf,-) grid

Curve Fit Multivariate Least Square

- Fitting Data to a plane
- The equation can be solve for cs as
- overdetermined
- mgtn

A c Z c Z/A

Curve Fit Multivariate Least Square

- clf
- x-414'
- y4 -3 5 -4 1 -3 4 -1 3'
- z18.74 -1.10 19.88 -5.71 6.20 -10.37 4.96 -5.30

1.54' - Ax y ones(size(x))
- cA\z
- xglinspace(-5,5,10)
- ygxg
- X,Ymeshgrid(xg,yg)
- Zc(1)Xc(2)Yc(3)
- plot3(x,y,z,'k.')
- hold
- surf(X,Y,Z)
- grid
- hold off

Polynomials

- Polynomials are widely used in engineering
- How does Matlab represent a polynomial
- How can we use this?

p1 3 -6 -2 7 -12 p p 1 3 -6

-2 7 -12

polyval(p,1) ans -9 polyval(p,-1) ans

-13 polyval(p,3.32) ans 537.4667 polyval(p,0)

ans -12

Polynomial Operations

- There are a number of things we might like to do

with polynomials - find roots roots(p)
- reconstruct from roots poly(p)
- multiplication conv(p1,p2)
- division q,rdeconv(pnum,pdenom)
- addition or subtraction youre on your own!
- derivative polyder(p)
- rational well leave this for an advanced class!
- curvefitting polyfit(x,y,ndeg)

Types of Interpolation

- Interpolation involves computing approximate

values between limiting or endpoint values - constant (zeroth order, nearest neighbor)
- linear
- cubic
- spline
- Equivalent to what is often called table lookup
- Matlab yiinterp1(x,y,xi,method ) function
- x,y are vectors defining the table of values
- xi is the independent variable
- yi is the interpolated result
- method is linear (default), nearest, cubic,

spline - this function can treat xi,yi as vectors

Interpolation

- This example compares different interpolation

methods for 6 points defining a sine function

echo on hold off clf xlinspace(0,2pi,6) ysin(x

) plot(x,y,'or') hold on pause xilinspace(0,2pi

,100) yiinterp1(x,y,xi,'nearest') plot(xi,yi) p

ause yi2interp1(x,y,xi) plot(xi,yi2,'g') pause y

i3interp1(x,y,xi,'cubic') plot(xi,yi3,'k') pause

yi4interp1(x,y,xi,'spline') plot(xi,yi4,'c') pa

use plot(xi,sin(xi),'r') legend('data','nearest',

'linear','cubic','spline','sine') xlabel('X'),

ylabel('Y')

2D Interpolation interp2()

- A bivariate function z-sin(x2y2) is interpolated

on the square -1?x ?1, -1?y ?1 using linear and

cubic methods - x,ymeshgrid(-1 0.251)
- zsin(x.2 y.2)
- xi, yimeshgrid(-10.051)
- ziinterp2(x,y,z,xi,yi,linear)
- surf(xi, yi, zi),title(Bilinear-interpolant to

sin(x2y2)) - ziinterp2(x,y,z,xi,yi,cubic)
- surf(xi, yi, zi),title(Bicubic-interpolant to

sin(x2y2))

Numerical Integration

Numerical Integration

Numerical Differentiation

- .Backward Difference
- .Forward Difference
- .Average of mA and mB
- .diff() function Central Diff.
- x0pi/50pi
- nlength(x)
- tdcos(x)
- ysin(x)- 0.05(rand(1,51)-0.5)
- d1diff(y)./diff(x) Backward Diff.
- subplot(2,1,1), plot(x(2n),
- td(2n),x(2n),d1,o)
- xlabel(x), ylabel(Derivative)
- axis(0 pi 2 2)
- 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(De

rivative), - axis(0 pi 2 2),title(Central Diff.),grid

Numerical methods for Root finding

- Numerical methods are recipes that allow us to

solve generic problems, e.g. Root finding - Bisection, false position, secant method,

Newtons Method - To use these methods effectively requires we
- Have a basis for comparing between alternatives
- Know when each will and wont work and why
- Study of numerical methods involves
- Development of algorithms
- Study of their properties, e.g. convergence,

stability, etc.

Recap - Newtons method

- The iteration formula for Newtons method is
- Can think of Newtons method as a fixed point

iteration scheme with

- The derivative of g(x) is

- Hence for Newtons Method to converge

Roots of Polynomials and function

Secant methodRoots

General form of iteration equation

Secant method .Roots

Point The secant method can be interpreted as

Newtons method with an approximated gradient.

Roots .

Differential Equations

- Ordinary differential equations (ODEs) arise in

almost every field - ODEs describe a function y in terms of its

derivatives - The goal is to solve for y

Numerical Solution to ODEs

- In general, only simple (linear) ODEs can be

solved analytically - Most interesting ODEs are nonlinear, must solve

numerically - The idea is to approximate the derivatives by

subtraction

Matlabs ODE solvers

- Matlab has several ODE solvers
- ode23 and ode45 are standard RK solvers
- ode15s and ode23s are specialized for stiff

problems - several others, check help ode23 or book

Euler Method

Euler Method

- Simplest ODE scheme, but not very good
- 1st order, explicit, multi-step solver
- General multi-step solvers

(weighted mean of f evaluated at lots of ts)

Runge-Kutta Methods

- Multi-step solvers--each N is computed from N at

several times - can store previous Ns, so only one evaluation of

f/iteration - Runge-Kutta Methods multiple evaluations of

f/iteration

Ordinary Differentional Equation-ODE

ODE Example

Example differential equation

General ODE This example is model of an

electronic circuit with a nonlinear vacuum tube,

called van den Pol equation.

Example differential equation

- Construct an m-file function for the right-hand

side of

fvdp.m

function dxdt fvdp(t,x) dxdt fvdp(t,x)

van der Pol equation dxdt x(2,) ...

(1-x(1,).2).x(2,)-x(1,)

Solve the ODE in the time span 0lttlt10 t,x

ode45(fvdp,0,10,20) plot(t,x) plot(x(,1),

x(,2))

OPTIONAL MATERIAL Solving ODEs

- Matlab includes a number of functions to solve

Ordinary Differential Equations (ODEs),

including Initial Value Problems (IVPs),

Boundary Value Problems (BVPs) and Partial

Differential Equations (PDE) - Lets consider a simple IVP in the form of a

familiar ODE (an sdof vibration problem) - Matlabs ode23() and ode(45) functions use the

Runge-Kutta-Fehlberg method to solve ODEs

expressed as

This is advanced material you will cover in your

next math classes. Matlab will be there when you

need it!

Solving an ODE setup

- We can convert an Nth order ODE into N first

order ODEs using a simple algorithm - In more compact forms

or

Matlab ode45( ) Syntax

gtgt help ode45 ODE45 Solve non-stiff

differential equations, medium order method.T,Y

ODE45(ODEFUN,TSPAN,Y0) with TSPAN T0 TFINAL

integrates the system of differential equations

y' f(t,y) from time T0 to TFINAL with initial

conditions Y0. Function ODEFUN(T,Y) must return a

column vector corresponding to f(t,y). Each row

in the solution array Y corresponds to a time

returned in the column vector T. To obtain

solutions at specific times T0,T1,...,TFINAL (all

increasing or all decreasing), use TSPAN T0 T1

... TFINAL. (truncated)

- T,Y are the returned values and each row

defines a value of t where the solution is

computed along with the corresponding solutions,

yi , in successive columns. - But we need to provide a function to compute

f(t,y) whenever ode45( ) needs it - We also need to specify the start and end times

and initial conditions in the arguments to ode45(

)

rhs( ) function for ode45( )

- m-function will compute f(t,y) for ode45( )
- returns the RHS column vector

NOTEIf you need to pass parameter values to

compute the RHS (e.g, zeta or g(t)), these can be

added to the ode45( ) function call (see help

ode45)

Solving the Problem

See help ode45 for more options

gtgt tt,yyode45('rhs', 0 35,1 0') gtgt whos

Name Size Bytes Class tt

205x1 1640 double array yy

205x2 3280 double array Grand total

is 615 elements using 4920 bytes gtgt

plot(tt,yy(,1))

You supply this m-function

- Note the sizes of the returned variables
- You can plot either column of yy as needed
- How would you construct a phase plane plot (e.g.,

y versus y)?

A More Interesting RHS

- Note how g(t) is formed here

Result is familiar square pulse with ringing

oscillations

Simulink Example

Best thing to do is to go through an example

2nd order, constant coefficient, linear

differential equation

Response to a step command

Simulink Example

Get an equivalent block diagram for the system

use mouse to drag blocks into the model window

and to connect blocks with arrows

use integrators to get dy/dt and y

Simulink Example

add gain and summer blocks

Simulink Example

add the step input block

Simulink Example

add the output block

Simulink Example

Now, double click the blocks to open and set the

blocks parameters

set gain value

set initial condition

set variable name

set output format to array

Simulink Example

To set the simulation parameters.

select Simulation -gt Simulation Parameters

set Start and Stop time (in seconds)

set numerical integration type

Simulink Example

Time to run the simulation

click the run button to begin the simulation

when the simulation is complete, Ready appears

at the bottom

Simulink Example

Simulink will automatically save a variable named

tout to the workspace.

This variable contains the time values used in

the simulation, important for variable time

integration types

Simulink also will create the output variable(s)

you specified

Simulink Example

plot(tout,yoft)

graph of the step response

Simulink Example

Another approach to solving the 2nd order single

DOF problem, is to cast it as a 1st order 2 DOF

problem

In Matrix (or State Space) form.

Simulink Example

1st Order State-Space Models

Simulink Example

Multi Input Multi Output Systems

use Mux and Demux blocks to combine and extract

vector signals

specify number of signals

Non-Linear Function

Non-Linear Function

Non-Linear Function

Non-Linear Function

END