Tip:
Highlight text to annotate it
X
In the previous lecture, we have been exploring how to use Markov property of response processes
associated with dynamical systems, which are driven by white noise excitations. We will
continue that topic.
So, what we discussed in the previous lecture was, we considered differential equations
governing the dynamical systems and cast them in the form of Ito stochastic differential
equations. And, for n-dimensional dynamical system, the equation has this form -- f is
the drift factor; G is the diffusion matrix; dB is increments of Brownian motion processes.
And, we showed that formal derivative of a Brownian motion process can be interpreted
as Gaussian white noise. And, the time illusion of the conditional probability density function
of the system states conditioned on the initial conditions is given by this partial differential
equation. And, this alpha j and the alpha i j are the incremental moments, which are
related to the drift and diffusion as shown here. In the previous lecture, I demonstrated
how the governing equation, the so-called equation Fokker-Planck equation can be derived
for different dynamical systems. And, this equation can be viewed as equation of motion
for time evolution of probability density function.
The application of Markov vector approach requires that the excitations have to be modeled
as white noises or as outputs of dynamical systems to white noise inputs. And, the equations
need to be represented in the state space form. So, these are some of the primary requisites
for applying the Markov vector approach.
We will consider few questions now in this lecture. How to solve the governing FPK equations?
And, I will present a few selected examples for which exact solutions are possible. And
then, we will also consider the question, how can we derive equations for evolution
of response moments? The equation gives the equation for time evolution of probability
density function. But, if you are happy with lesser level of description of response, it
is useful to ask what the equations that govern the response moments are. So, we will consider
these questions in the present lecture.
Now, one of the mathematical tools that I will be needing is the so-called Lagrange's
method for solving linear partial differential equations. This you must have studied in your
under graduate curriculum. So, we consider partial differential difference equations
of the form P into dou z by dou x plus Q into dou z by dou y is equal to R of x, y, z. Now,
to obtain an integral of the above equation, we consider the auxiliary equation dx by P,
dy by Q and dz by R. This is a sense of Lagrange's method. And from this set of equations, we
construct two independent solutions to be u of x, y, z equal to a and v of x, y, z equal
to b; where a and b are constants. Then, the solution to this partial differential equation
can be expressed as some phi (u, v) equal to 0. Alternatively, u as a function of v
is also a solution. So, we will be using this in our solution of FPK equations.
A quick example, a simple example say xz dou z by dou x plus yz dou z by dou y equal to
xy. Now, the auxiliary equation here will be dx by xz, dy by yz and dz by xy. Now, we
consider these two pairs of equations and the first two of these equations lead to the
condition that y by x is a constant. And similarly, the first and the third leads to the equation
z square minus xy is b. Therefore, the general solution of this partial differential equation
can be written as phi of y by x comma z square minus xy equal to 0. So, this can be verified.
I leave that as an exercise. But, you need to recall some of the details associated with
this method.
Now, we will start by considering linear dynamical systems under additive white noise. For which
we know the exact solution through a normal mode expansion and convolution integral approach.
We have already discussed that solution. But, in this part of the discussion, we will consider
how we can use the Markov property of the response vector and obtain the solution by
solving the governing Fokker-Planck equation. So, we consider n degree of freedom system;
X is N cross 1 and it is acted upon by m cross 1 vector of white noise excitations and gamma
is a kind of influence matrix, which is n by m. And, system starts from initial condition;
X of 0 is X naught and X dot of 0 is X naught dot. W of t is 0 mean white noise; its covariance
is given as 2D ij, which is m by m matrix, which is a direct delta function.
Now, first step in implementing the Markov vector approach for solving this problem would
be to recast the problem in the state space form. So, to achieve that, what we do is we
premultiply this equation by M inverse and I get X double dot plus M inverse CX dot plus
M inverse KX is equal to M inverse gamma W. Now, I introduce the state vector Y 1 and
Y 2 as X and X dot and recast this equation into a set of 2N first order equations; dY
I by dt is actually X dot; X dot is Y double I dt; and, d Y double I dt is obtained from
this equation and it is displayed here Now, these two equations can be combined into a
single form, which is dY t is minus PY dt plus QdB t with y of 0 be Y naught. The matrix
P here is a 2N by 2N matrix, which consists of M inverse K, M inverse C; and, Q is the
multiplier for the noise stumps and this is again is 2N cross m matrix.
So, the governing equation is displayed in the standard format as shown here and the
sizes of various quantities are displayed here. And, this equation dY t is minus PY
plus QdB t can be interpreted as an Ito stochastic differential equation of size 2N cross 1.
And, response vector Y consisting of Y I and Y double I will have Markov property and we
can derive the associated Fokker-Planck equation for the time evolution of the conditional
probability density function of Y conditioned on Y of 0. Before we do that, we would like
to diagonalize this matrix P that would facilitate certain numerical calculations. So, we consider
the eigenvalue problem P phi is lambda phi and we denote by capital phi the 2N cross
2N matrix of eigenvectors. And, capital lambda is the 2N cross 2N diagonal matrix of complete
set of eigenvalues of P. Now, P phi is phi lambda. Therefore, phi inverse P phi will
be capital lambda; and, capital lambda is a diagonal matrix. So, now, it could be this;
I introduce a transformation Y as phi Z.
And, use this relation phi inverse P phi is a diagonal matrix and uncouple the equation.
So, I substitute that; phi dZ t is minus P phi Z t dt plus QdB t. Now, I premultiply
by phi inverse; and, phi inverse phi is I. And, here I get phi inverse P phi Z plus phi
inverse QdB as an equation; and, phi inverse P phi is capital lambda, which is a diagonal
matrix. And, this phi inverse Q, I call it as G, which is 2N cross m matrix. This alpha
j is the drift coefficient, which will be minus lambda into z. And, alpha ij are GDG
transpose. These are the coefficients in the Fokker-Planck equations. So, the Fokker-Planck
equation itself is given in this form; where p is the conditional probability density function
of z evaluated at t conditioned on z at t equal to 0. The initial condition is in terms
of products of direct delta functions; and, boundary conditions, we assume that z i tends
to plus minus infinity; the function goes to 0.
Now, we need to solve this problem. So, this is a linear partial differential equation
with constant coefficients. So, the integral transform methods appear eminently suited
to handle this problem and indeed that is the way we look at it. We will use a characteristic
function, which is a fourier transform of p and try to solve this problem.
So, to facilitate that, we define the conditional characteristic function M theta 1, theta 2,
theta 2N. I denote this as M of theta tilde colon t. Now, this is actually equal to the
expectation of exponential i j to 1 to 2N theta 1 z j. This is the definition of the
characteristic function. So, this I write compactly in terms of a vector integral; z
tilde is a vector. And now, if I differentiate this with respect to theta, I get z k p k
-- z k p exponential this. See I would need this, because I need in this equation the
fourier transform z j into p, is what I would be needing. So, that I have to handle. Therefore,
you can see that the fourier transform z k p will be the derivative of the characteristic
function with respect to theta k.
So, what that means is M is given in terms of this. Therefore, the probability density
function is inverse fourier transform defined with respect to characteristic function. And,
since dou M by dou theta k is given by this expression, I get the fourier transform z
k into p, will be the inverse transform of this function, which is 1 by 2 phi i, etcetera
here. So, I am now ready with using characteristic function to solve the Fokker-Planck equation.
There I would need dou by dou z k of this and that is given in terms of characteristic
function as shown here.
So, now, I will also need a second derivative, which is shown here and derivative with respect
to time, which is shown here. And, I substitute all these into the governing Fokker-Planck
equation.
And, I get the equation in terms of characteristic function as shown here. Now, this equation
is again a partial differential equation, where independent variables are now thetas
and time. Thetas are the parameters appearing in the definition of characteristic function.
So, this equation I will now solve using Lagrange's method. So, I will write the auxiliary equation
in this form as shown here and I will solve them pairwise; and, construct the solution
and substitute into the general solution here, the final expression here and we will be able
to tackle this. Some of these details you can understand by studying the accompanying
powerpoint slides. So, I will briefly run through this. So, I will write this equation
as dt by 1 is d theta i by lambda i theta i. And, from this, I get theta i t is theta
i 0 exponential lambda i t for i running for 1 to 2N. This I recast in the matrix form
as shown here, where capital omega is a diagonal matrix with entries exponential minus lambda
i t.
Now, I consider the last term in that auxiliary equation and try to solve for M. So, I consider
two terms: d theta i by lambda i theta i and d theta l by lambda l theta l. And, I will
manipulate. Suppose if this equation is a and this equation is b. I will multiply a
by lambda i theta i theta l; and, equation b by lambda l theta i theta l. And, that leads
to left-hand side to be this and right -and side to be this. So, this can be recast as
d of theta l theta i divided by lambda i plus lambda l; and, this expression on the right-hand
side. Now, I do further manipulations. I multiply both sides by lambda i l by 2 and sum over
i and l to get the equation dM by M is equal to this.
And, this I can solve for M and I get M as in terms of M naught exponent of this. Now,
I will introduce capital XI to denote lambda i l by lambda i plus lambda l. And, from this,
I get the characteristic function in a compact form as M naught exponential minus half theta
tilde transpose XI theta tilde. Now, we have to determine M naught and for that, I use
initial conditions. And, if you manipulate this, we will be able to find M naught.
And finally, if you run through all these calculations, you can show that the characteristic
function here corresponds to that of a multidimensional Gaussian random vector. This is to be expected,
because we have a linear system, which is driven by Gaussian excitation. So, we know
by an alternate means, which we have already studied that response is a Gaussian random
process. So, we expect beforehand based on that knowledge that a characteristic function
corresponds to the double Gaussian random process.
Now, once this is z coordinate, if I take the inverse transform, I will get the Gaussian
probability distribution function, density function. And, from this characteristic equation,
I will also be able to find mean vector and covariance matrix. And, once I find the mean
and covariance matrix, I can use the transformation, which I introduced. Our original variables
were Y, but the solution has been obtained in Z; the transformation can be made directly
at the level of mean vector and covariance and I will get the requisite multidimensional
normal Gaussian density function for the response process. This is valid for all time; it is
not just for steady state. So, it is a complete solution. So, this type of complete solution
for Fokker-Planck equation is rarely obtainable; linear system under additive Gaussian noise
is one instance. It is not always possible to do this.
So, we will make some remarks. For linear systems, the exact solution can also be obtained
using convolution integral approach discussed earlier in this course. The Markov vector
approach does not offer any special advantage here, because we already know how to solve
this problem by an independent means. Now, I consider one of the limitation that we could
attribute to Markov vector method, is that the excitations need to be Gaussian white
noises. But, this is not a very stringent requirement. In the sense, if you do not have
a white excitation, you can always make a model for the given random excitation by passing
white noise through a filter. For example, if you consider now x double dot plus 2 eta
omega x dot plus omega square x equal to f; f is not a white noise process, but it is
output of a linear system to white noise excitation. Suppose there is another single degree freedom
system with natural frequency lambda and damping xi; receives white noise; and, output of this
is f of t. So, f of t can always be modeled as output of linear systems to white noise
excitations; under various general conditions it is possible.
Now, what happens is I define now an extended vector that consists of x, x dot, f and f
dot. So, at the expense of handling a larger dimensional state space I will be still be
able to use Markov vector approach. So, here obviously, the vector consisting of x and
x dot will not have Markov property. But, if you consider the extended vector x, x dot,
f, f dot, that vector will have Markov property. And therefore, I can cast this pair of equations
into a single Ito equation, which is of the form dX of t is equal to minus PX dt plus
QdB t. And, I have already discussed how to solve this problem; that means the formulation
that I just now presented can also be used if excitation needs not white, but it can
be modeled as output of a linear system, which is driven by white noise excitations. And,
since the solution is valid for all times, if there is any non-stationarity in terms
of a deterministic time envelope multiplying these excitations, the method that we have
discussed should be possible in principle to use that for this kind of excitations also.
Now, what are the other types of problems that we can solve using Fokker-Planck equation?
By solve I mean exactly. Now, I will consider a sequence of problems and just indicate how
the problems could be tackled. So, we return to the problem of a first order nonlinear
systems -- x dot plus beta of x is driven by white noise, w of t. And, this is an Ito's
representation and this is the Fokker-Planck equation. In the previous lecture, I indicated
that the transient solution to this can be obtained by using method of variable suppression.
And, by using the eigenfunction expansion, we would be able to solve this problem. And,
there are problems of this kind, which have been tackled. I will briefly mention that
later. But, if we restrict our attention to only
steady state; suppose the system admits steady state. It is not necessary that every system
should admit a steady state; it depends on nature of beta of x. Suppose the system admits
steady state, then as time becomes large, dou p by dou t will be equal to 0. p is actually
the probability density function of x at the time t -- single time t. So, if a process
is stationary, that probability density function would become independent of time. Therefore,
dou p by dou t would be 0. In that case, I get a simplified equation or a reduced equation,
which essentially characterizes only the stationary solution. So, if you are not interested in
transient, if you interested only in steady state, I can directly tackle this equation.
Since in this equation, the two independent variables are x and t, and by bringing in
the motion of steady state, I have eliminated one of the independent variables, namely,
t, time. I am left with ordinary differential equation of the form shown here
Now, this can be solved under the boundary conditions that x becomes plus infinity and
minus infinity; this function goes to 0. And, I can rewrite this equation as D of dp by
dx is equal to minus beta of x into p, because there is D by dx that we can put that out
and consider this equation. Now, this equation can easily be solved and I get p of x to be
C into exponential minus 1 by D 0 to x beta of s d s. Now, we have to select this C, such
that the normalization condition is satisfied. We will just give you an example; if beta
of x is equal to a of x plus bx cube say, then the probability density function will
have this form. If of course, b equal to 0, we get the Gaussian
density function, which is expected because if it is x dot plus alpha x equal to w of
t, then the solution is known to be Gaussian. You can find out its variance; you can do
power spectral analysis or you can use Duhamel integral analysis and find out the steady
state solutions. So, that solutions would essentially match with this solution with
b equal to 0. This is as it should be . Because system is linear and it is driven by Gaussian
white noise, the response is going to be Gaussian. And, in fact, it will exactly match with the
solution that you will get by using convolution integral approach.
We will consider now single degree freedom systems with nonlinear damping and nonlinear
stiffness. The equation of motion that we are considering is x double dot into x dot
some function of H -- I will explain what this H is -- plus g of x is equal to w of
t; t is greater than or equal to 0 and initial conditions are specified. w of t as 0 mean
and its covariance -- it is a delta correlated process. Therefore, expectation of w t into
w of t plus tau is 2D delta of tau. H is actually the energy in the system, which is kinetic
energy -- x dot square by 2 plus the strain energy -- 0 to x g of u du. It is the total
energy. Assume that the nonlinear damping is a function
of the energy. It is still nonlinear, but it obeys this particular functional form.
If all these assumptions of this model are respected, then it is possible to show that
the steady state response of this system can be exactly determined by using Fokker-Planck
equation approach. Let us see how we can do that. I introduce a vector X 1, X 2 state
vector, which is displacement and velocity. I recast the equation into the Ito differential
form; dX 1 of t is X 2 dt and dX 2 of t is minus X 2 -- this is f of H minus g of X 1
dt plus dB of t. Therefore, H is X 2 square by 2 plus 0 to x 1 g of u du. So, this is
the statement of the governing equation of motion.
Now, what is the governing Fokker-Planck equation? It is dou p by dou t equal to minus x 2 dou
p by dou x 1 minus dou by dou x 2... This is alpha 1; this is alpha 2; alpha 1 -- you
can see is actually X 2 and alpha 2 will be minus X 2 f of H and minus g of X 1. These
are the drift coefficients. And, we get those first two terms here and the diffusion is
capital D, because white noise is multiplied by a constant unity. Therefore, the diffusion
term will simply be D. So, p is a transition probability density function as displayed
here. And, initial conditions are in terms of product of two direct delta functions and
the boundary conditions we assume that at x 1 and x 2 plus minus infinity, the probability
density function goes to 0. Obtaining transient solution for this problem
is not straight forward. Currently, no exact solution exists for the complete solution
of this problem, exact solution But, an approximate solution can always be generated, for example,
by using method of weighted residuals or finite element method. You can directly tackle this
partial differential equation and approaches, such as that have been discussed in the existing
literature. But, in this discussion, we will restrict our attention to steady state. That
would mean dou p by dou t is 0. So, I am interested in now, p would be become the joint density
function between displacement and velocity evaluated at the same time. And, in this steady
state, that would be independent of time. Therefore, dou p by dou t would be equal to
0. So, I am left with only this right-hand side, which is equal to 0.
This equation is still a partial differential equation, because there are two independent
variables here. In the original problem, there are three independent variables: t x 1 and
x 2. By assuming steady state, I have reduced one of the dimensions, namely, time. But,
still I am left with two independent variables: x 1 and x 2. So, we need to solve this problem.
Now, we expand this and slightly rearrange these terms and bring it in the form as shown
here -- minus x 2 dou p by dou x 1; this is dou x 2; second term -- g of x 1 is dou p
by dou x 2 plus D into dou square p by dou x 2 square equal to 0. Now, we are looking
at stationary solutions. So, x of t as t tends to infinity becomes a stationary random process.
For a stationary random process, the process and its derivative are uncorrelated at the
same time. So, if process is Gaussian, it automatically implies that x and x dot are
independent also. But, here process is not Gaussian, but still x and x dot are uncorrelated.
But, we could explore... If x and x dot, we can find out a solution, where a kind of variable
separation becomes possible. So, what we will demand is we will demand that the first term
and this term; sum of these two is suppose 0; and, the remaining two terms are separately
0. I am looking for a special solution. Now, we will consider the first part of this
equation If these two are independently 0, and I find a p that would automatically satisfy
this equation also; the same p, which satisfies these two equations would satisfy this also.
Now, we will consider first of this equation -- minus x 2 dou p by dou x 1 plus g of x
1 dou p by dou x 2 equal to 0. Now, using Lagrange's equation, Lagrange's approach,
we can write this as dx 1 by minus x 2; dx 2 by g of x 1; dp by 0 -- this is auxiliary
equation for this problem.
If I consider now the first two of these equations, I get dx 1 by minus x 2 is equal to dx 2 by
g x 1. And, from this, I get the condition x 2 square by 2 plus 0 to x 1 g of u du is
a constant a. Now, dx 2 by g x 1 is dp by 0. Therefore, p itself should be a constant.
So, based on these, I can write a form of a solution, which is p -- some arbitrary function
capital phi of this function -- x 2 square by 2 plus 0 to x 1 g of u du. The argument
of this function in fact is H. So, I can write this as phi of H. We will now consider the
next equation and write it as dou by dou x 2 of minus x 2 F of H p plus D dou p by dou
x 2 square equal to 0.
Now, I will demand that the term inside the parenthesis is 0 and I will get... Now, for
p, I write phi of H and go through this equation. And, I will now get an equation, which is
d phi by dH is equal to minus 1 by D F of H phi of H. This phi of H now can be solved,
because this is now independent variable; here is H. So, we can solve this. And therefore,
p x 1 comma x 2 will be expressed in this particular form; where this H, which appears
as a limit, is this function. So, this is a deliberate effort to consider only that
class of dynamical systems for which the Fokker-Planck equation can be solved. So, if you are given
a dynamical system and if you set up a Fokker-Planck equation, that equation depending on the nature
of nonlinearity etcetera may not be in a form that permits exact solution. So, I am not
discussing generic procedures to solve Fokker-Planck equation, but I am doing a kind of an inverse
investigation on what could be the nature of nonlinear dynamical systems, which permit
exact solutions in steady state to the governing FPK equations. So, that is the theme of our
discussion.
Now, we can start by considering F of H to be 2 eta omega. That means H to the power
of 0 some constant into H to the power of 0. So, this is a class of system, where damping
is linear, but our spring force or the elastic forces are non-linear. So, this could model,
for example, structure with geometric nonlinearity. So, we will show that these types of problems
are amenable for exact solutions. So, we can write now... Since we have got the exact solution,
I will now write in the exact solution, F of psi to be 2 eta omega and I get this as
the solution. And, in terms of x and x dot, this will be the solution.
Now, you can see here that this form of the density function is such that it can be expressed
as product of a function of x alone and a function of x dot alone, because this is exponential
of a term involving x dot plus an exponential of term involving only x. So, it can be viewed
as product of exponentials, of two exponentials: the first one is only function of x dot; the
second one is function of x. Or, in other words, the form of this expression suggests
that the joint density function, x and x dot are stochastically independent.
Now, to make these discussions slightly more specific, we can consider the duffing system,
where the system has cubic nonlinearity. So, this is the joint density function between
displacement and velocity of a duffing's oscillator driven by white noise. The solution is in
steady state. This is an exact solution. A duffing system under deterministic excitations
admits no exact solutions; whereas, under white noise excitations and if you focus on
steady state solutions, it admits an exact solution. So, this is a very interesting feature
associated with using Markov method for solving this class of problems. This constant phi
naught is still floating around and that we have to select, so that the probability density
function is suitably normalized; the volume under the probability density function is
unity. So, we will make quick remarks. This is an
exact solution. Now, if alpha equal to 0, if this time goes off, then this corresponds
to a 2-dimensional Gaussian density function, which is in fact exact solution for a linear
single degree freedom system under white noise excitation. So, this we have derived by using
convolution integral approach earlier. So, we can verify that the solution that we are
getting is exactly the same as what we got by an alternate approach. Next, we can verify
that the joint density function can be written as product of p x 1 and p x 2; that would
mean that x and x dot are independent. We know that x of t and x dot of t are uncorrelated,
but in this particular case, they are also independent and they are non-Gaussian.
Now, we can verify that the probability density function of displacement has this form. It
has in the exponent, a term containing x to the power of 4. That would mean this is non-Gaussian;
whereas, velocity, that is, probability density function of velocity at a given time instant
t is Gaussian. We should carefully interpret this result; we should not conclude that velocity
is a Gaussian random process. If velocity were to be a Gaussian random process, it automatically
implies displacement is also a Gaussian random process, because velocity and displacement
are related through linear transformation. But, we know that displacement is a non-Gaussian
random process. So, we can only say that the velocity is a Gaussian random variable while
displacement is a non-Gaussian random variable. It should be noted that velocity is not a
Gaussian random process. If it were so, the displacement would also be Gaussian, which
is not the case. Therefore, here in this case, velocity is a non-Gaussian random process
with a first order PDF, that is, Gaussian. So, it is an interesting feature about the
response.
Now, one more remark that we could make is, in the expression for H, the first term is
kinetic energy and the second term is potential energy. So, this represents total energy.
And, you look at the joint probability density function between displacement and velocity;
we get in the exponent, an exponential distribution with respect to energy. That means H is exponentially
distributed. So, these types of results are encountered in Maxwell-Boltzmann equations
for gas dynamics and so on and so forth. So, it has certain resemblance to those physical
laws.
This is a plot of joint density function between displacement and velocity for a nonlinear
duffing oscillator. It appears Gaussian-like; it is uni-modal.
And, if you look at the marginal density function of the displacement, this red line is a Gaussian
density function and the blue line is the result for PDF of displacement. So, clearly,
we can see that it is non-Gaussian although it has a uni-modal character.
Through this analysis, valid albeit only for the steady state; we are still able to determine
the joint density function between displacement and velocity at the same time t. Now, if you
recall, much of the discussion that we made concerning level crossing statistics, first
passage time and maximum, the basic quantity of interest was the join density function
between process and its derivative at the same time. So, it is important to note that
the Fokker-Planck equation actually provides that. So, for a duffing oscillator under white
noise, I have exactly got the join density function between displacement and velocity
at the same time. Therefore, now I can solve for number of times the level alpha is crossed
in 0 to t; I can find out its average. And, if levels are high, I can use a Poisson model
for that and derive the probability distribution function for the first passage time. And,
based on the distribution for PDF of first passage time, I can also get a model for X
t. So, that would mean the joint density function between process and its derivative at the
same time is the key descriptor if you are interested in reliability related measures
of the response.
So, for this particular example, we can actually find out the average rate of the crossings
of level alpha with positive slope. And, we get for the duffing system in terms of the
non-Gaussian density function and an integral on a Gaussian density function. So, this can
easily be evaluated and we can go ahead and do the other calculations.
Now, we will consider another example, where F of H now I take it as 1 minus x square and
x dot square; and, g of x to be linear. So, this is a nonlinear system with linear stiffness
characteristic, but with nonlinear energy dissipation characteristics. Again, it is
driven by white noise, delta correlated and all that. And, H is a total energy now, which
is x dot square by 2 plus x square by 2. Since the stiffness terms are linear, we get this
expression. So, we can go ahead and find out the joint density function between displacement
and velocity in the steady state. And, that has this form.
Now, I can substitute for F of H, which is actually minus 1 minus 2H, because this is
x square minus x dot square minus of that; it is 1 minus 2H. So, that is what I have
written here. And, if you now substitute into this known exact solution, I get this final
answer, which is the joint density function between the displacement and velocity in steady
state. Now, you carefully look at this equation; you will see that x and x dot are not independent.
You can find out marginal densities of x and x dot and multiply them and you will be able
to show that the product would not lead to this.
Now, x and x dot are non-Gaussian, but they are dependent although they are uncorrelated.
x and x dot need to be uncorrelated, because the process in which steady state... For a
stationary random process, a process and its derivative should be uncorrelated at the same
time instant. So, this is an example, where process is non-Gaussian, that is, uncorrelated;
x and x dot are uncorrelated, but they are jointly non-Gaussian and dependent; mutually
dependent. Now, we need to take a look at the system, in the absence of noise, how does
the system behave. I will discuss this now and show that the joint probability density
function has a model line in the neighborhood of the limit cycle. And, we will see that
the PDF of displacement and velocity, which are the marginal density functions, are not
uni-modal.
So, to see that, we need to look at this system bit more closely. If I were to plot the phase
plane, that is, x of t and x dot of t, I will plot and I will eliminate time. And, I consider
free vibration characteristics of these systems here. You look at the second term here. For
small oscillations when x and x dot are small, the term inside this parenthesis would be
positive, because x and x square plus x dot square will be still less than 1. And, therefore,
the net sign of this is negative. Therefore, small motions tend to grow. That means origin
is unstable. But, if x and x dot are large enough so that the term inside this parenthesis
is positive, the system becomes positively damped and all large oscillations tend to
decay. So, small oscillations grow and large oscillations decay. And, what happens is small
oscillations grow will like this and large oscillations will decay like this; and in
between, there will be a closed trajectory. This is known as limit cycle. That means for
this system, the free vibration consist of oscillatory motions.
In fact, you can show that for this system, x of t is cos t plus phi; where phi is an
arbitrary phase angle. You substitute into this; you should be able to verify that this
is an exact solution, because you will immediately see here 1 minus x square minus x dot square
is 0 and x double dot plus x is satisfied by this equation. Therefore, this equation;
this is the solution for this equation . So, this is actually a circle in the phase plane
plot. So, this is the status of the system behavior when there is no noise. Now, upon
this if you impose a noise, you could expect that most of the motion will take place around
the steady state, will be around the limit cycle, the stable limit cycle. So, if you
plot the probability density function of x, you may expect this kind of probability density
function. There will be two modes: x... And, this corresponds to the amplitude of the limit
cycle system. In the absence of noise, the motion is essentially circular. So, you add
a small noise, the motion tends to remain around that circle.
So, if you plot the joint probability density function of displacement and velocity in presence
of noise, you see here a kind of a behavior, where there is a depression at x at the origin;
and, this circular feature is the reminiscent of the limit cycle in absence of noise. Since
for this system origin in absence of noise is unstable, there will be minima in the probability
density function at the origin. Unlike when we get the Gaussian models, where origin is
stable and a imposed noise, you get a mode in the probability density function, a maxima;
whereas for systems where the origin is unstable, we expect a minima in the probability density
function and maxima elsewhere, where there are stable fixed points are limit cycles.
So, that feature is being displayed here.
And, this is a probability density function of the displacement process. You can see here
that it is essentially non-Gaussian. And, there seems to be several modes here, but
if you carefully plot, it may have appearance something like this with two modes. So, as
noise becomes small, indeed there will be this kind of behavior in the probability density
function. So, this is one of the features of effect of noise on limit cycle systems.
And, if you vary some of the parameters of the nonlinear systems and see how the quality
of the response changes, in the contest of deterministic analysis, you will see that
you will study fixed points and their stability. And, as the parameter of the problem is varied,
the number of fixed points and their nature of their stability would change. And, that
we call as bifurcation. These bifurcations manifest in the stochastic case in terms of
maxima and minima of their associated probability density functions. And, that part of the interpretations
can be made from the exact solutions that we have obtained for this system.
Now, that was a single degree freedom system under white noise excitation, a class of nonlinear
models, where nonlinearity could be in damping or in stiffness. How about multi-degree freedom
systems? What type of problems are amendable for exact solutions? In fact, the family of
problems, which are amendable for exact solutions by their time means the steady state solutions
derivable from FPK equations are exact. That is a wide clause; I am just speaking few examples
from it. So, here I consider a multi-degree freedom system, where non-linearity in stiffness
and that is in terms of derivative of potential energy, dou U by dou X j. So, these equations,
the coefficient of accelerations and velocities, the matrices are diagonal. Nevertheless, all
those equations are coupled through this nonlinear term.
This W j of t is a vector of white noise processes for j
running from 1 to n. And, we assume that they are all independent. W j t is independent
of W k t. And also, we impose another condition that the parameter m j, which is the initial
parameter and D jj, which is noise parameter; and, beta j, which is the damping parameter.
For all j, this ratio is constant, which is independent of j. Now, under only
these exceptional situations, we get an exact solution.
Now, we can rewrite this equation in terms of state space form. I will first use first
n terms for displacements and the next n terms for velocities and we get this equation. And,
we can derive the drift and diffusion terms in terms of the coefficients of the governing
Ito differential equation. And, I will leave it as an exercise for you to show that the
join density function between n displacements and n velocity processes is given by this
expression in the steady state. This is a fairly comprehensive solution, but valid only
for a limited class of systems under certain idealizations. So, this type of exact solutions
will serve as benchmarks to test approximate methods. And also, it gives us insights into
the influence of noise on system behavior. Because these solutions are exact, the interpretations
we make have substantial value.
So far, we have been discussing about the evolution of probability density function.
Now, how about moments? Can we directly derive the equations for moments. So, that we can
do. Conceptually, it is not very difficult. So, we start by considering this, where there
is a non-linear drift and diffusion terms and I write the governing FPK equation like
this. Now, suppose if I am interested in this expectation, X 1 to the power of k 1 X 2 to
the power of k 2, etcetera, I can write in terms of the probability density function
and the initial probability... the transitional probability density function and the initial
density function as shown here. If I differentiate this with respect to time, I get m dot and
this is dou by dou t of p of this transitional probability density function. This is given
by the right-hand side of the Fokker-Planck equation. So, I can substitute there.
And, in principle, I can get this equation, because this dou by dou t of p is given in
terms of the two terms in the FPK equation. And, this can be integrated by parts and we
can get the moment equations. Same principle is possible. So, conceptually, there is no
difficulty to just derive the time evolution of moments.
A more general approach would be to consider an non-linear function h of X comma t. We
will assume that h is well-behaved; that is, in the sense, its derivatives -- second derivative
with respect to X i and X j and first derivative with respect to time exist. And, we consider
now this increment delta h, which is h of delta plus h of X plus delta X, t plus delta
t minus h of X comma t. Based on Taylor's expansion, I derive delta h in terms of h
and its gradients as shown here.
And now, if I consider the conditional expectation of delta h conditioned on X, this is a kind
of terms that I will need in the FPK equation. We can show that this is nothing but expectation
of this right-hand side . Now, we will digress briefly. If you consider two random variables
x and y and if you find conditional expectation of Y conditioned on X equal to x, it is given
by this. This actually is a random variable. And, if you take its expectation, it will
be actually expected value of Y. Therefore, this is a conditional expectation. Now, I
will take conditional expectation of this.
If I consider the expectation of this, my idea is to get expected values of this, the
governing equation for this. So, if we use this idea, we can show that the moment equations
are given by the equation shown here. This uses the properties of the Brownian motion
process and their increments. And, we are also using the drift and diffusion coefficients
in the Fokker-Planck equation and we can show that these are the exact equations for time
evolution of the moments.
Now, a simple example would throw light on this suppose if I consider X dot plus beta
X is equal to w of t. dX t is minus beta X t plus dB t. So, f will be minus beta x. This
is the drift term and this is the equation for the evolution of expectation of h of X
comma t. Now, suppose h is X to the power of k; m k, I denoted as expected value of
X k. So, m k dot is given by this equation. So, you can substitute these terms here and
you will see that m k dot is this. Now, k equal to 1. m 1 dot is minus beta m 1; m 2
dot is minus beta m 2 plus 2D; m 3 dot -- I have to put here; X k equal to 3 and we will
be able to solve this and I get a series of equations, which I can solve. Now, what happens
if we look at steady state? So, the time derivative, all these moments vanish and I am left with
an algebraic set of equations, which can be solved.
So, I immediately get the mean response, which is 0; the variance, which is this; skewness,
third moment, which is 0; and, fourth moment, which is this. And, we show that this is actually
three times sigma x to the power of 4, which is property of a Gaussian random variable.
So, in this particular case, is very straightforward to write these moment equations
. And, if you are interested only in steady state, you are only going to solve algebraic
equation. So, there is no convolution integral, there is no selection of partial differential
equation; you can directly get the steady state solution in a straightforward manner.
Now, one thing that we should notice here is that the moment equations here are closed.
In the sense, you want to find the mean, you need not have to know any moment, which is
of higher order. And, another thing that we should notice is that we know that the response
is Gaussian and it is displaying features; the response moments are displaying the features
of Gaussian random variable with 0 mean, for which mean is 0, skewness is 0 and variance,
the fourth order moment is three times the 4th power of standard deviation. So, this
is skewness, is 3, what I am getting. So, this matches with known features of Gaussian
response.
Now, we will consider the familiar single degree freedom system under white noise. So,
we can go through this formulation, write the Ito's equation and write these moment
equations. You can write equation for.... So, here I am introducing the notation -- this
notation -- X 1 to the power of m X 2 to the power of n; expected value I am writing as
m mn.
So, m 10 is expected value of X 1; m 01 is expected value of X dot. This is nothing but
expected value of X. So, m 20 is expected value of X square. m 11 is XX dot, so on and
so forth. So, you can write these equations. This is variance of mean square value of the
velocity. So, on the right-hand side, you see that if you are looking at mean, m 10
is coupled to m 01, but equation for m 01 has only m 10 and m 01. So, these two equations
can be solved together and I can find out the mean. At the level of characterizing mean,
I need not have to know anything more than the mean of the response. Similarly, second
order moments, m 20, m 11, m 02 -- all these three can be solved using the properties at
the second order level. So, these can be solved. So, we say that the moment equations are closed.
Now, steady state -- the time derivative of the moments vanish; mean of displacement is
0; mean of velocity is 0. The correlation coefficient between displacement and velocity
is 0, because the response is in steady state. Now, if we analyze this, we get that the variance
of the displacement is D by 4 eta omega cube; and, the variance of velocity is D by 4 eta
omega, and the process and time derivative are uncorrelated. These are the exact solution
that we obtained by two different methods: one by using convolution integral; other by
doing spectrum analysis. So, now, this is the third approach, which leads to the same
answer. So, these results agree with the exact solutions obtained earlier using convolution
integral and spectrum analysis approaches. So, it is quite satisfying.
So, in the next lecture, what I will do is I will consider these moment equations for
slightly more general class of problems; and then, we will also consider questions on first
passage times -- how to use Markov process theory to characterize first passage times.
And next, we will also consider questions on how to characterize enveloped processes
using Markov property. So, we will take these topics in the next lecture and we will conclude
this lecture at this juncture.