Transcript Lecture30.pptx

Recap
 Cubic Spline Interpolation
 Multidimensional Interpolation
 Curve Fitting
 Linear Regression
 Polynomial Regression
 The Polyval Function
 The Interactive Fitting Tools
 Basic Curve Fitting
 Curve Fitting ToolBox
 Numerical Integration
Numerical Integration
Example
 Here’s another example, using a
function handle and an anonymous
function, instead of defining the
function inside single quote
 First define an anonymous function
for a third-order polynomial
 fun_handle = @(x)-x.^3+20*x.^2-5
 Now plot the function, to see how it
behaves. The easiest approach is to
use fplot , since it also accepts a
function handle:
 fplot(fun_handle,[-5,25])
Example Continued….
 The integral of this
function is area under the
curve
25
 −5 −𝑥 3 + 20𝑥 2 − 5
 Finally, to evaluate the
integral we’ll use the
quad function, with the
function handle as input:
 quad(fun_handle,0,25)
 ans = 6.3854e+003
Solving Differential Equation
Numerically
 MATLAB includes a number of functions that solve ordinary
differential equations of the form

𝑑𝑦
𝑑𝑡
= 𝑓 𝑡, 𝑦
numerically
 In order to solve higher-order differential equations they must
be reformulated into a system of first-order expressions
 Not every differential equation can be solved by the same
technique, so MATLAB includes a wide variety of differential
equation solvers
 However, all of these solvers have the same format. This makes
it easy to try different techniques by just changing the function
name
Continued….
 Each solver requires the following three inputs as a minimum:
 A function handle to a function that describes the fi rst-order
differential equation or system of differential equations in terms
of t and y
 The time span of interest
 An initial condition for each equation in the system
 The solvers all return an array of t- and y -values:
 [t,y] = odesolver(function_handle,[initial_time, final_time],
[initial_cond_array])
 If the resulting arrays [t,y] are not specified, the functions
create a plot of the results
Function Handle Input
 A function handle can refer to either a standard MATLAB function, stored as an Mfile, or an anonymous MATLAB function
 For differentiation equation

𝑑𝑦
𝑑𝑡
= 𝑓 𝑡, 𝑦
a function handle is equivalent to 𝑑𝑦
𝑑𝑡
 Here’s an example of an anonymous function for a single simple differential
equation:
 dydt = @(t,y) 2*t corresponds to
𝑑𝑦
𝑑𝑡
= 2𝑡
 Although this particular function doesn’t use a value of y in the result (2 t ), y still
needs to be part of the input
 If you want to specify a system of equations, it is probably easier to define a
function M-file
 The output of the function must be a column vector of first-derivative values, as in




function dy=twofuns(t,y)
dy(1) = y(2);
dy(2) = -y(1);
dy=[dy(1); dy(2)];
Continued….
 This function represents the system
𝑑𝑦

𝑑𝑡
𝑑𝑥

𝑑𝑡
=𝑥
= −𝑦
which could also be expressed in a more compact notation as
 𝑦1 ′ = 𝑦2
 𝑦2 ′ = −𝑦1
where the prime indicates the derivative with respect to time,
and the functions with respect to time are 𝑦1 , 𝑦2 , and so on
 In this notation, the second derivative is equal to 𝑦 ′′ and the
third derivative is 𝑦 ′′′ :
 𝑦′
=
𝑑𝑦
,
𝑑𝑡
𝑦 ′′
=
𝑑2𝑦
,
𝑑𝑡 2
𝑦 ′′′
=
𝑑3𝑦
𝑑𝑡 3
Solving the Problem
 Both the time span of interest and the initial conditions for each
equation are entered as vectors into the solver equations, along
with the function handle
 To demonstrate, let’s solve the equation
𝑑𝑦

𝑑𝑡
= 2𝑡
 We created an anonymous function for this ordinary differential
equation in the previous section and called it dydt
 Evaluate y from -1 to 1 and specify the initial condition as
y(-1)=1
 If you don’t know how your equation or system of equations
behaves, your first try should be ode45 :
[t,y] = ode45(dydt,[-1,1],1)
 This command returns an array of t -values and a
corresponding array of y -values
Continued….
 User can either plot these
or allow the solver
function to plot them if
user don’t specify the
output array:
 ode45(dydt,[-1,1],1)
 The results are consistent
with the analytical
solution, which is
 𝑦 = 𝑡2
 The first derivative of this
function is 2 t and that y
= 1 when t = -1
Continued….
 When the input function or system of functions is stored
in an M-file, the syntax is slightly different
 The handle for an existing M-file is defined as
@m_file_name
 To solve the system of equations described in twofun we
use the command
 ode45(@twofun,[-1,1],[1,1])
 The time span of interest is from -1 to 1, and the initial
conditions are both 1
Solving Higher-Order Differential
Equations
 A higher-order differential equation can be expressed as a series of equations by making some
simple substitutions
 Consider the following equation:

𝑑2 𝑦
𝑑𝑡 2
+
𝑑𝑦
𝑑𝑡
=𝑦+𝑡
 We can reformulate it into a system of equations by introducing a new variable, z
 Let
 𝑧=
𝑑𝑦
𝑑𝑡
 It’s then easy to see that

𝑑𝑧
𝑑𝑡
=
𝑑2 𝑦
𝑑𝑡 2
 Substituting into the original equation we get

𝑑𝑧
𝑑𝑡
+
𝑑𝑦
𝑑𝑡
=𝑦+𝑡
which is a first-order differential equation
 Effectively we’ve replaced

𝑑2 𝑦
𝑑𝑡 2
+
𝑑𝑦
𝑑𝑡
=𝑦+𝑡
with the following two equations, which have been rearranged to solve for the first derivative
of our two dependent variables, y and z

𝑑𝑦
𝑑𝑡
=𝑧
𝑑𝑧
𝑑𝑡
=𝑦+𝑡−
𝑑𝑦
𝑑𝑡
Continued….
 Now all we need to do is create an M-file function to use in one of




the ode solvers
The function should have two inputs, which are typically called t and
y
The variable t is the independent variable, and the variable y is an
array of dependent variables
In this example y (1) corresponds to the y used in the hand
formulation, and y (2) corresponds to z
The function containing the system of equations should look like
this:




function dydt = twoeq(t,y)
dydt(1) = y(2);
dydt(2) = y(1) + t - dydt(1);
dydt = dydt'
 Notice that the function output has been formulated as a column
vector, as required by the ode solvers
 Also recall that the function name is arbitrary. We could have called
it anything, but twoeq is descriptive
Continued….
 Once the system of equations is defined in a function M-
file it is available to use as input to an ode solver
 For example: if the range of time is defined as -1 to +1
and the initial conditions are defined as y = 0 and z=0,
then the command becomes
 ode45(@twoeq,[-1,1],[0,0])
which gives the results
 A problem where the starting values are known is called
an initial value problem
Boundary Value Problems
 bvp4c function is used to solve boundary value problems
 The bvp4c function requires three inputs:
 A function handle to the system of ode’s to be solved
 A function handle to a function that solves for the residual
values of the function
 A set of guesses for the initial conditions
Continued….
 The first function handle is exactly the same as we used for the ode solver set of
functions
 It should contain the equations for the derivatives of interest and the results must be
a column vector
 To solve the problem a guess is made for the initial value of all the derivatives, then
the program checks to see how it did by comparing the calculated boundary values
with the actual values
 For example, if:
 at t = -1, y = 0 and
 at t = 1, y = 3
 the program would solve the system of equations based upon an initial guess of d y
/d t , and would then check to see how close the result is at t = -1
 This is accomplished using a boundary condition function where the equations are
arranged so that if the correct boundary condition is calculated, the function values
are zero
 In the case of our example,
 function residual=bc(y_initial, y_final)
 residual(1) = y_initial(1) + 0;
 residual(2) = y_final (1) - 3;
 residual = residual';
Continued….
 If this function is executed for values of y_initial = 0 and
y_final = 3, the result will be a column of zeros
 Any other result means that the program has calculated the
wrong values for y_initial and y_final , and the guesses for
the initial conditions must be updated according to the
function’s algorithm, which is a finite difference strategy
Continued….
 The last input to the bvp4c function is a mesh of guesses for the problem







solution, which are used as the starting point in the solution
MATLAB provides a helper function, bvpinit , to help create this mesh,
which is stored as a structure array
It requires two inputs;
an array of values corresponding to the independent variable
initial guesses for each of the variables defined in the ode system of
equations
In our case there are two equations, so we’ll need a guess for y and d y /d t
The mesh need not be particularly fi ne, and the initial guesses need not be
very good
For example:
 initial_guess = bvpinit(-1:.5:1, [0, -1])
 specifies five t values from -1 to 1 (-1, -0.5, 0, 0.5, 1) and initial guesses of
y = 0 and dydt = -1 at all values of t
Continued….
 Once the function describing the system of ode’s, the function defining the residuals and the
initial guesses created with bvpinit have been created, the bvp4c function can be executed.
 bvp4c(@twoeq, @bc, initial_guess)
 which returns
 ans =
 x: [1x9 double]
 y: [2x9 double]
 yp: [2x9 double]
 solver: 'bvp4c'
 The result is a structure array, where x is the value of the independent variable and where an
array of y values corresponds to the solutions to the system of ode’s. In this case, y and dy /dt .
 To access the array of x values simply use the structure syntax, ans.x
 If we had chosen to assign a name such as solution to our result instead of defaulting to ans ,
the structure would be called solution , and the x values would be stored in solution.x
 The values of most interest are the y values, which can also be accessed using structure syntax,
such as solution.y
 To plot the results in a manner similar to that displayed by the odesolvers use the code
 plot(ans.x,ans.y, '-o')
or, if the results were named solution
 plot(solution.x, solution.y, '-o'),
which gives the results
Summary of Chapter
 Tables of data are useful for summarizing technical information
 However, if a value is needed that is not included in the table,
must approximate that value by using some sort of
interpolation technique
 MATLAB includes such a technique, called interp1 . This
function requires three inputs: a set of x -values, a
corresponding set of y -values, and a set of x -values for which
you would like to estimate y -values
 The function defaults to a linear interpolation technique, which
assumes that approximate these intermediate y -values as a
linear function of x that is,
y = f(x) = ax + b
Continued….
 A different linear function is found for each set of two
data points, ensuring that the line approximating the data
always passes through the tabulated points
 The interp1 function can also model the data by using
higher-order approximations, the most common of which
is the cubic spline
 The approximation technique is specified as a character
string in a fourth optional field of the interp1 function
 If it’s not specified, the function defaults to linear
interpolation. An example of the syntax is
new_y = interp1(tabulated_x, tabulated_y, new_x, 'spline')
Continued….
 In addition to the interp1 function, MATLAB includes a two-





dimensional interpolation function called interp2 , a threedimensional interpolation function called interp3 , and a
multidimensional interpolation function called interpn
Curve-fitting routines are similar to interpolation techniques
However, instead of connecting data points, they look for an
equation that models the data as accurately as possible
Once you have an equation, you can calculate the
corresponding values of y
The curve that is modeled does not necessarily pass through the
measured data points
MATLAB’s curve-fitting function is called polyfit and models
the data as a polynomial by means of a least-squares regression
technique
Continued….
 The function returns the coefficients of the polynomial equation of
the form
 𝑦 = 𝑎1 𝑥 𝑛 + 𝑎2 𝑥 𝑛−1 + 𝑎3 𝑥 𝑛−2 + ⋯ … … … … + 𝑎𝑛 𝑥 + 𝑎𝑛+1
 These coefficients can be used to create the appropriate expression in
MATLAB, or they can be used as the input to the polyval function to
calculate values of y at any value of x
 For example: the following statements find the coefficients of a
second-order polynomial to fi t the input x–y data and then calculate
new values of y , using the polynomial determined in the first
statement:
 coef = polyfit(x,y,2)
 y_first_order_fit = polyval(coef,x)
 These two lines of code could be shortened to one line by nesting
functions:
 y_first_order_fit = polyval(polyfit(x,y,1),x)
Continued….
 MATLAB also includes an interactive curve-fitting capability that allows
the user to model data not only with polynomials, but with more
complicated mathematical functions
 The basic curve-fitting tools can be accessed from the Tools menu in the
figure window
 More extensive tools are available in the curve-fitting toolbox, which is
accessed by typing
 cftool




in the command window.
Numerical techniques are used widely in engineering to approximate both
derivatives and integrals
Derivatives and integrals can also be found with the symbolic toolbox
The MATLAB diff function finds the difference between values in adjacent
elements of a vector
By using the diff function with vectors of x - and y -values, we can
approximate the derivative with the command
 slope = diff(y)./diff(x)
 The more closely spaced the x and y data are, the closer will be the
approximation of the derivative
Continued….
 The gradient function uses a forward difference approach to




approximate the derivative at the first point in an array
It uses a backward difference approach for the final value in the
array, and a central difference approach for the remainder of
the points
In general, the central difference approach gives a more
accurate approximation of the derivative than either of the
other two techniques
Integration of ordered pairs of data is accomplished using the
trapezoidal rule, with the trapz function
This approach can also be used with functions, by creating a set
of ordered pairs based on a set of x values and the
corresponding y values
Continued….
Integration of functions is accomplished more directly with one of two

quadrature functions: quad or quadl
 These functions require the user to input both a function and its limits of
integration
 The function can be represented as a character string, such as
 'x.^2-1‘
as an anonymous function, such as
 my_function = @(x) x.^2-1
or as an M-fi le function, such as
 function output = my_m_file(x)
 output = x.^2-1;
 Any of the three techniques for defining the function can be used as input,
along with the integration limits—for example,
 quad('x.^2-1',1,2)
 Both quad and quadl attempt to return an answer accurate to within 1 ×
10−6
 The quad and quadl functions differ only in the technique they use to
estimate the integral
 The quad function uses an adaptive Simpson quadrature technique, and the
quadl function uses an adaptive Lobatto quadrature technique
Continued….
 MATLAB includes a series of solver functions for first-order
ordinary differential equations and systems of equations
 All of the solver functions use the common format
 [t,y] = odesolver(function_handle,[initial_time, final_time],
[initial_cond_array])
 A good first try is usually the ode45 solver function, which
uses a Runge–Kutta technique
 Other solver functions have been formulated for stiff
differential equations and implicit formulations
 The ode solver functions require that the user know the initial
conditions for the problem
 If, instead, boundary conditions are known at other than the
starting conditions, the bvp4 function should be used