 If you recall in our first lecture we said data mining, data assimilation and prediction are parts of a continuum. Data mining relates to development of models almost all the models that come to be originally arose from considerations that can be related in one way or the other to the concept of data mining as we understand it today. Once the notion of a model was well established the notion of being able to use the model to create prediction came about. Once we are able to predict the interest became improving the quality of prediction became fundamental became a problem of fundamental importance. How do you improve the quality of prediction? We need to use the data. So data assimilation is a process by which you can fit the model to the data so that the model can be calibrated to make reasonably good prediction. So I have created the model I have observed the data I have created the assimilated model what is next to create prediction. So in this course we are going to see the various aspects of the quality of prediction who controls the quality what is predictability what are the why is the interest in predictability studies and that is the theme of this last set of modules. So first we would like to be able to split our discussion of predictability into deterministic predictability where deterministic refers to prediction generated by a deterministic model. Stochastic predictability relates to predictability of analysis of predictions arising from stochastic models. So as a starting point we are going to talk about what is called error dynamics or variational equation which we have already alluded to when we talked about 40 war methods. So we will quickly start developing some of the basic tools that we would need in the analysis of deterministic predictability that xk plus 1 is equal to m of xk alpha be a deterministic model x0 is the initial condition alpha are the parameters. Now let us consider the model running from one state called x0 bar another state called x0 so that is one prediction there is another prediction the trajectory is also called orbits. So let xk bar let the sequence xk bar the sequence xk be so what is the sequence x1 bar x2 bar xi bar here is x1 x2 xi. So let us consider two orbits or two trajectories generated by the model starting from two different initial conditions let epsilon k is equal to xk minus xk bar therefore the initial difference is epsilon not the initial difference is epsilon not then this difference is epsilon 1 this difference is epsilon i this difference is epsilon k. So what is epsilon k epsilon k is the difference in the forecast of the model at time k induced by the difference epsilon not in the difference in the initial condition that is used to generate the two different forecasts. So epsilon k is called the so if epsilon not is the error in the initial condition epsilon k is the induced error in the state of the system at time k. Why is this important? In many of the 4D war based methods as well as forward sensitive methods within the context of deterministic dynamic data simulation scheme our goal has been to be able to estimate the initial condition. The estimate of the initial condition based on finite samples will always have errors therefore if you used one set of observation I will get one estimate for the initial condition if you use another set of observations I will use another set of initial conditions. So no matter which method you use if you are trying to estimate the initial condition from a finite set of observation the estimate will not be equal to the true value that I would need to be able to make better forecast. In other words estimates have always errors embedded in them. So you can think of epsilon not as the error in the initial condition estimate. So if I want to give some life to this trajectory let X not bar be the true unknown initial state. So this could be the true trajectory this could be the trajectory this is the predicted trajectory this could be the predicted trajectory. The predicted trajectory is different that true trajectory because the initial state that is used to generate the prediction was different from the true state the initial state estimate. So you can think of X not as the estimate of the initial condition the epsilon not is the error in the initial condition. So the error in the initial condition however small it is it is going to be reflected in the forecast I am interested in the error dynamics. So the evolution of epsilon induced by the model is called the error dynamics that is the title of the slide. So I would like to be able to now derive an equation for the evolution of epsilon based on the evolution of the true state as of the predicted state. So what is that one would expect if the error in the initial condition is small if the prediction at time k the error is closer to the true state that means the error in the prediction at time k is also small in other words small initial errors leads to small errors in the forecast at future time then you would say the model forecast is more reliable. If the error in the initial state explodes in time then the quality of prediction given by the model deteriorates in time. So what does this relate to it relates to the sensitivity of the model to the variations are errors in the estimate of the initial condition. So if the initial errors are magnified by the model that means the model is very sensitive if the initial errors do not grow but grow goes down to 0 that means the models are more stable. So stable models can be used to create better prediction but not all models are stable. So the quality of prediction depends on essentially to the sensitivity of the model to errors in initial condition this is the fundamental theme of predictability analysis this is the fundamental theme of predictability analysis. So that is the pathway we are going to be taking in this set of slides that describes analysis of deterministic predictability. So let xk plus 1 actual model using is the actual model starting from the state x naught now to a first order approximation I can say the unknown true state is equal to the known prediction plus epsilon k. So it is again a first order theory. So xk is equal to this xk can be replaced by the true state plus epsilon k alpha. So here what is that I am assuming alpha is the parameter you know if I change the alpha the solution changes. So I am going to pretend for a time being I know alpha precisely and I have used alpha the precise value of alpha if there is going to be a forecast error it is all due to only errors in the initial condition and not the errors in the parameters. So I want to separate I do not want to assume errors in too many things again I want to emphasize the forecast errors can come from errors due to three sources one errors in the initial condition errors in the parameter are errors in the model I want to be able to enjoy analyze each one of them separately. So in this case what is that we are interested in assume the model is perfect deterministic assume I know all the parameters in the model absolutely precisely. So if there is any forecast error that is directly attributable to only errors in the initial condition and the errors in the initial condition where do they come from they come from estimates while data simulation schemes gives you reasonably good estimates optimal estimates in some sense because the estimates are derived out of finite number of samples the estimate from the finite set of samples may only be closer to the true value still that could be a non-zero error our aim is to be able to see how the model treats this more often inevitable error the error coming from estimation of the unknown initial condition based on finite samples that is the kick that is the real key. So if I now express this map in the form of a Taylor series it can be seen this is the base value this is the perturbation so from here if I from here if I so please also remember xk plus 1 bar is equal to m of xk bar plus alpha therefore this equation this equation becomes this because this is equal to xk plus 1 bar xk plus 1 minus xk plus 1 bar is epsilon k the right hand side is becoming the Jacobian of m at time k at time k means what at evaluated at xk times epsilon k so this becomes the dynamics for the evolution of the error where dk is essentially d of xk of m I did not want to complicate the notation so d of xk of m is simply dk so what does this tell you epsilon 1 is equal to d 0 m epsilon 0 epsilon 2 is equal to d 1 m epsilon 1 which is equal to d 1 m times d 0 m times epsilon not what is this this is the Jacobian evaluated the first state this is the evaluation at the next state so it is the product of Jacobians this equation 3 in mathematics is called the variational equation in meteorological it is also called the tangent linear system so tangent so in the context of in the context of 4d war we consider tangent linear system essentially with respect to understanding the propagation of perturbation here there we induced we thought of inducing a perturbation initially to be able to correct the forecast errors here there is no correction there is no data there is no we have already done the data simulation so I have a model x naught be the estimated initial state I am going to run the model from this estimate I am going to compare the performance of an assimilated model with respect to the unknown true state and that is the analysis we are trying to do so this is the post data simulation analysis of model forecast to be able to see how the initial errors are treated by the model so I am naturally interested in 3 3 efforts to the dynamics of evolution of the model errors so if I iterate 3 I get this expression what is this expression if I picked 2 instances in time yeah I am sorry this is 0 this is time yes and this is time k I would like to be able to express how the error at time yes is related to the error at time k and that is the relation so yes could be 0 or anything else so this is general expression though this tells you how the errors at 2 instances in times separated in time are are are are related by the dynamics d d of k colon s of m is simply the product of the Jacobian evaluated along the trajectory from time s to time k this matrix has a special name it is called a propagator matrix what does it do it tries to it tries to relate the errors between 2 different instances in time the propagator matrix is simply a product of the Jacobians now I am going to introduce a numerical measure please understand yep the error epsilon k so I have now understood how the errors propagate to a first order accuracy through the model now but the errors are vectors we would we understand certain scalar measure better than vectors so I am now going to introduce a scalar measure that tries to that tries to quantify the properties of errors and that is done by a quantity called Rayleigh coefficient the Rayleigh coefficient is simply the ratio of the square of the energy norm of the error at time k plus 1 to that at time 0 so in the in the in the in the previous slide we related the error at time k plus 1 to error at time s error time k plus 1 epsilon k plus 1 is a vector error time s is a vector s in general could be 0 or anything any number less than greater than 0 so s is greater than 0 less than or equal to k so that is essentially the range for s so epsilon k plus 1 is given by the product of Jacobian from k to s of m times epsilon s energy norm I am trying to evaluate a numerator with the energy norm based on an s p d a the bottom line s p d b so a and b let a and b be two given symmetric positive matrices evaluate the numerator using the energy norm based on the matrix a evaluate the denominator using the energy norm based on the matrix b in principle a and b could be equal but the treatment is simple is as simple when you do not have to consider a and b for example what could be the thing in meteorology I would I would like to be able to say hey initially I am interested in pressure but at time at a later time I may be interested in certain vorticity initially I may be interested in temperature difference at a later time I may be interested in rain so what do a and b bring to the bring to back on the problem I can I would like to be able to consider the sensitivity of the rain at a future time and it depends on a temperature distribution of the previous time I could be I may be I may be interested in some quantity at time at a later time and its sensitivity based on another quantity therefore at this stage when I am trying to talk about quantities sensitivities I do not have to make them to be the same therefore Rayleigh coefficient what is it it is simply the ratio of quadratic forms what is the quadratic form so the numerator so rk plus 1 epsilon is essentially equal to epsilon k plus 1 transpose a epsilon k plus 1 divided by epsilon s transpose b epsilon s that is exactly that is exactly the relation that that that is involved in here now please understand epsilon k plus 1 and epsilon s are related through the model the relation through the model is given by the product of the Jacobian along the lines so how does the initial so what is that you can think of it now suppose I have an grid I have an initial condition I have a model I am I am going to introduce a thermal bubble in a small scale in the at the initial time in involving certain small number of grid points so there is a perturbation of temperature distribution initially how does the model reacts to the thermal bubble at a time in the future how does this thermal bubble introduces other changes in the model so that could be one way of thinking about thinking about this ratio so A and B in general are two symmetric positive matrices this ratio is called is the celebrated Rayleigh coefficient in matrix theory and that comes to our that is a very useful measure for us to be able to consider predictability analysis so I would like to be able to understand the behavior of this ratio now let us see what this ratio represents this ratio in general represents a measure of the sensitivity you can you can see the energy in some quantity at a later time related to the energy of some other quantity at the initial time so this ratio captures the spirit of the sensitivity that is involved in the predictability analysis suppose this ratio remains less than 1 this ratio is a scalar right if this ratio remains small for all time beyond as what does it mean the error does not grow if the error does not grow what does it mean the predictions are pretty done accurate on the other hand if this ratio grows in time what does it mean the errors the initial errors magnify therefore analysis of the properties of this Rayleigh coefficient for appropriate choices of energy norms of the numerator the denominator essentially provides a very good clue to the quality of forecast generated by the model that is the that is the that is the secret of using Rayleigh coefficient so in order to further analyze the Rayleigh coefficient let us concentrate on the propagator matrix what is the propagator matrix it is a product of the Jacobian from time s to time k let us call that matrix A for the sake of simplicity in notation I am also going to assume the model Jacobian is non-singular if the model Jacobian non-singular means what the model is well formed non-singularity the Jacobian at every point along the trajectory goes to attest to the well formed nature of the model itself the model sometimes the model could be screwed up how do you measure what is the measure of the screw up in the model look at the Jacobian the model actually gives you a solution you try to evaluate the value of Jacobian along the trajectory if the value of the Jacobian remains full rank the model is well formed if the model if the rank of the Jacobian along the model varies then the model is not well formed model so I am assuming the model is well formed in the measure in the sense that Jacobian along the trajectory are non-singular okay now A is the square matrix please please please realize that A is a square matrix even though A is a square matrix in general A is not symmetric A depends on the starting time s and the starting time k so if I vary any one of the timings A changes but for the analysis I am keeping s fixed k fixed so A is fixed so A transpose A and A transpose are the two Grammians coming out of this model both of them are spd why both of them are spd if each of the components of a component matrices in this product are non-singular the product is non-singular if the if a matrix is non-singular its Grammian is full rank the Grammian is not only full rank is also symmetric and positive definite these are all immediate conclusions we have alluded to several times earlier especially within the context of development of svd singular value decomposition so let V1 to Vn be the Eigen values this is the Eigen vectors of the matrix A transpose A so A transpose AV is equal to V lambda VV transpose V transpose V is I now I am going to define a new vector UI which is equal to 1 over lambda I am using this similar formalism as we did in the svd the only differences earlier we defined svd for the context in the context of matrix H which is a triangular now I am talking svd within the context of the matrix A which is square H is a single matrix with forward operator A is a complex matrix is a product of matrices along the trajectory if I assume each of the matrices each of the Jacobian matrices are non-singular A is non-singular so I am just trying to reinforce the idea where A comes from so if VIs and lambda Is are the Eigen vector Eigen value pair for A I can now define a new make vector UI UI also is an n vector which is 1 over square root of lambda I times AVI from this derivation so I am now going to compute A transpose UI if you substitute the expression for UI and simplify it it turns out I get 9 so what does 9 say A A transpose UI is equal to lambda I UI so what does it mean if VIs are the Eigen vectors of A transpose A UIs are the Eigen vectors of A A transpose both the matrices share the same set of Eigen values. In this context the lambda Is are the Eigen values of A transpose A as well as A A transpose both of them are matrices of the same size so square root of lambda I is called the singular value that is the definition. So singular value of a matrix is the square root of the Eigen value of the Gramian that is the definition and the square root of the Eigen values of the Gramian are all going to be positive if the Gramian is symmetric positive definite a Gramian is symmetric positive definite if the component matrices are a full rank so full rank condition essentially tells you the problems are well formed and that is the condition we have been looking at all through full rank matrices all good things in life and that is exactly the kind of theory we have developed. So square root of lambda I are called the singular values VIs are called the right or forward singular vectors UIs are called left or backward singular vectors these are the general nomenclature that is well understood within the applied mathematics linear algebraic context. Now from 8 from 8 I would like to talk about some of the properties of equilibria from 8 I can rewrite my equation for the definition of UI this way there are answers if I collected all the relations the matrix form it becomes a matrix version so 11 is the matrix version of 10 therefore A transpose A is given by this decomposition because UU transpose U transpose U they are all identity U is orthogonal V is orthogonal therefore I get this expression so the columns of V are orthonormal system so what is that we are now going to do we have been doing analysis in the general coordinate system given by the standard basis which is EI what is the standard basis EI is equal to 0 1 0 and this is the F location I for 1 to n that is the standard basis if I have any other basis I can do the analysis in that basis why would you change the basis if the analysis in the new basis can be made simple nothing is lost in fact lot may be gained by changing the coordinate system from the standard to a given coordinate system now we are going to gain by changing the analysis from the given simple coordinate system to an orthogonal system defined by V the columns of V what are the columns of V they are the eigen vectors of the matrix A transpose A given a space that we should always have an n-n basis vectors if the basis vectors are orthogonal it becomes orthogonal basis if the basis vectors are not orthogonal it is it is a basis but doing arithmetic in non orthogonal basis is little bit more involved so we are simply changing from one orthogonal basis to another orthogonal basis and that is going to throw a lot of light on the behavior of the behavior of the Rayleigh coefficient so what is that what is the transformation we are going to be talking about now please understand epsilon s is the initial error at time s V is the new basis so I am going to transfer epsilon s to alpha so what is alpha alpha is the new vector alpha and epsilon n are referred to the same point in the space epsilon s is the coordinate of the same point in the original basis alpha is the coordinate of the point in the new basis whose columns whose whose whose basis vectors are columns of V and these two values of coordinates of the same point the two coordinate system are related by the expression in 13 therefore now I am going to re-express my value of the Rayleigh coefficient of time k plus 1 with respect to the Rayleigh coefficient of time s this is the definition of the Rayleigh coefficient in the standard basis which is 14 in the transformed domain my Rayleigh coefficient takes this form so what is that I am going to do please remember my a is equal to d colon s of m so a transpose is equal to d k s s transpose of m also remember my epsilon s is equal to is equal to V transpose alpha so if you substitute all this 14 becomes 15 14 becomes 15 what is the real kicker here the real kicker the the the the the numerator in the denominator get simplified what is that a transpose a V is equal to lambda V because V the columns of V are the what eigen vectors of that therefore if I multiply this by V transpose a transpose a V that is equal to lambda and that is exactly what is occurring in the numerator in the numerator therefore this complex matrix in the numerator becomes simplified as a diagonal matrix in the numerator we also know V transpose V is equal to V V transpose that again comes from the orthogonal eigen vector that is equal to I the denominator gets simplified by this therefore I have decoupled both the numerator and the denominator from the dependence and so in in in here what are the what are the various values of a I am assuming I am assuming a is identity I am I am sorry I am I am considering the energy norm I do not I do not think I want to go back right now sorry. So I hope this is clear now the the expression for the Rayleigh coefficient is given by 15 now what is the so denominator 15 what is that is the inner product of inner product of alpha with itself the the numerator is simply a quadratic form of alpha with a diagonal matrices so what is the so now what is that so alpha is a vector so alpha is any vector that is going to represent my my my epsilon s so I am going to now consider a normalized vector so the the the norm of alpha defined by this I am going to set to be 1 what does this mean if I have a vector here I can consider the normalized part of it. So that is essentially normalization the direction is more important than the magnitude itself that is the idea here so I am going to consider the ratio in 15 when I can find my attention to alpha unit vectors that is that is the whole that is the whole point of the game so in this case the numerator becomes 1 therefore rk plus 1 alpha is essentially is essentially the numerator which is alpha transpose lambda alpha which is again given by lambda I square root of square of alpha I square from here the following inequality becomes becomes follows directly this is the value the actual value since the sum since the sum of alpha I square is 1 since the sum of alpha I square is 1 so you can really see lambda I alpha I square I is equal to 1 to n is less than or equal to is less than or equal to maximum over I of lambda I times alpha I square I is equal to 1 to n and that is equal to 1 so that gives rise to this this is also right minimum over I of lambda I square times summation alpha I square I is equal to 1 to n that is equal to 1 therefore lambda 1 is the maximum eigenvalue lambda n is the minimum eigenvalue it readily follows rk plus 1 alpha is less than or equal to lambda 1 and lambda n where lambda 1 and lambda n are the maximum and minimum value of the eigenvalue so what is that we have done we have narrowed the values that the the ratio of the energy norms of the errors can take by by moving simply from the standard coordinate to the new coordinate formed by the eigenvectors of the Gramian A transpose A assuming A is full rank so this is a very very very very nice very nice treatment of the analysis of the analysis of the rather coefficient. Now what is that we are interested in we are interested in determining what is called the average rate of growth of errors from time s to k plus 1 then k equals to infinity that is of interest let us look at this why are we interested in the average growth of errors if I commit an initial error if the error grows how much error I can afford to correct my forecast such that forecast still meaning is still meaningful is the question for example in weather forecasting we have made tremendous improvement over the past 50 years but we still do not know how to make monthly forecast our models are some of the best known models in the last 100 years we have the best computers we have all kinds of observation systems we still cannot make predictions over the so the predictions that we make in today with all the things we know with all the technologies available computer sensors satellite models and everything else we can probably believe 3 to maximum 4 day forecast so what is 4 day 4 day is the predictability limit so what does it mean today I am going to make 1 day forecast 2 day forecast 3 day forecast 5 day forecast the error in 1 day forecast is there but smaller the error in 2 day forecast slightly larger error in the 5 day forecast is more why do you do not want to do the 10 day forecast is not that I do not know how to make a 10 day forecast I can make the 10 day forecast too but the error by the time 10 days into the future comes into being is so much that it dwarfs the signal so there is a signal on which the forecast error superimposed if the magnitude of the superimposing error is comparable to the signal then signal plus noise overwhelms the information so it is a classical study in engineering literature called signal to noise ratio so in any estimation problem in any detection problem we will always like to be able to analyze the signal to noise ratio and we would like to be able to maximize the signal to noise ratio so in some sense the Rayleigh coefficient tries to capture the spirit of the computation of signal to noise ratio so one measure one question is to how far does it take to be able to double the initial error so what is the norm of the initial error what is the norm of the error at a future time how long does it take for the initial error the norm of the initial error to double that is called error doubling time if the error doubles in 10 days after 10 days the error is even more than double so can I make a forecast that the error is doubled on the initial one no so what is the amount of error you may be able to tolerate how long how far the error grows how fast the error grows okay the error rate of growth varies from time to time to time so if the error rate of growth varies from time to time to time what is one measure average rate of growth of errors so to understand the quality of prediction we need to understand the predictability limit what is the predictability limit predictability limit is the limit beyond which errors overwhelm the signal if the errors overwhelm the signal the prediction is useless now here comes some of the basic things some processes are predictable perfectly the lunar solar eclipse IBM price I can predict for tomorrow but I don't think I can predict the price of an IBM stock day after tomorrow is much more volatile there are lots of factors that depend on how do you predict the price of a barrel of crude oil the barrel of crude oil depends not only on the cost of production but also the cost of transportation from place A to place B either you transport the unpurified base liquid or purified form but then the cost of crude oil is controlled by so many events in the world if there is a war in the Middle East or if there is a problem in the Middle East we believe that the supply may be may be may be affected the price immediately goes up so there are a number of factors that affect the price of a barrel of crude so the barrel of crude predictability of it is extremely difficult problem so likewise the weather prediction we have gained enough knowledge about the weather the behavior of it so that we can make good short term prediction long term prediction we try to do but we are not able to do long term prediction which are reliable for example I cannot predict seasonal weather precisely I do not let alone climate why it is not because we are not intelligent because the model we use exhibit extreme sensitivity to the initial condition so the initial condition you use to generate the forecast have even smaller errors if the model is very sensitive to small errors the errors blow up in time so that is the fundamental theme behind predictability analysis so the interest now goes to analyzing what is the average rate of growth of errors that is important so what is so given given yes as a starting point fixed this is k I would like to get to go to infinity so from a fixed starting point if I want to be able to make asymptotically large values of time if I want to be able to make forecast for such large values of time over the time what is the average rate of growth of errors as seen through the model to that end I am going to define a new matrix so what is this matrix you remember dksm which we have also called a I am going between a and this whenever we need to so this is simply a transpose a that is a gramium s is fixed it depends on m I am going to let k go to infinity so I am interested in the long term value of the gramium a transpose a I am also interested in the in the 2k-1 2 times k-1 the root of this matrix please understand please understand given a matrix a I can consider a square given matrix a I can consider a to the power of 1 half given a matrix a I can consider a to 1 over alpha for some alpha so this is a matrix I am I am I am I am considering 1 over 2 times k-1 s what is k-1 s k-1 s is the is the time difference between the starting point and the final time k goes to infinity so k-1 s goes to infinity so this this matrix plays a fundamental role the Eigen values of this limiting matrix are called Lyapunov vectors and the corresponding Eigen values of this limiting matrix are called Lyapunov numbers we are going to talk about these 2 but before we go there I would like you to be able to appreciate the definition of this matrix so what is this it is a product of the I am sorry it is a product of a transpose a a transpose a is the gramium a transpose a depends on s and and and k I am keeping s and m the model fixed I am letting k goes to infinity I am interested in the in the in the 2 times k-s the root of this product matrix and that is the matrix that limit matrix that limit matrix has an Eigen value Eigen vectors the Eigen values are called Lyapunov vectors Eigen values are called Lyapunov numbers Lyapunov is the famous Russian mathematician who who who who propounded the theory of stability in in in in the early decades of 20th century and and so in his honor it is called Lyapunov vectors Lyapunov numbers the natural logarithm of the Lyapunov numbers are called Lyapunov indices so let us look at this now so this is called Lyapunov indices Lyapunov indices so what is the Lyapunov index it is the let us look at this now let us look at this ratio what is this ratio it is the norm of the error at time so if I have time s if I have time k I have epsilon s I have epsilon k plus 1 the norm of the ratio of the so this is equal to d k times s of m of epsilon s so the norm of the error epsilon k plus 1 divided by epsilon s so what does it tell you by by running epsilon s through the model I am able to change the initial error so the the norm of this ratio is a number and this number depends on the time interval k minus s I am going to normalize this by 1 over k minus s s is fixed k goes to infinity so this ratio is the you can you can think of it as the average so this average value lambda is called the Lyapunov index and this Lyapunov index essentially tells you the average rate of growth of errors that makes sense I want to think about it now think about it now so epsilon k plus 1 by epsilon s what does it mean that is the rate of growth of error during the time from s to k I am trying to divide it by 1 over k minus s so that is the average rate of growth of errors so what is the meaning of this Lyapunov index I am going to introduce it by illustrating it on a scalar dynamics in fact the entire theory of dynamical chaos in deterministic system is centered around this notion of Lyapunov indices average rate of growth of errors sensitivity to initial conditions so what is the idea if a system exhibits is a extreme sensitivity to initial condition that system cannot be used to generate prediction for long intervals of time if the system exhibits less sensitivity to initial condition such systems can be used to create long term forecast it turns out most of the models of interest in climate studies most of the models of interest in oceanography atmospheric sciences exhibit do exhibit extreme sensitivity to initial condition that is the reason why we still are not being able to make 10 day forecast a monthly forecast seasonal forecast we still have not captured all aspects of the model sensitivity so the current models exhibit extreme sensitivity that is the reason why we are able to do what we are able to do so to understand this let us work as simple example let x dot is equal to f of x be a dynamical system x if x naught bar is a given initial condition that is the unknown that is the true value of the state so this is x naught bar so from here I can compute the solution x of t I have x naught this is x naught bar I can compute x of t so bar quantities are true unbarred quantities are coming from estimated values I have an initial difference which is epsilon naught I am doing the entire theory but some based in a simple one dimensional models so the dynamics of the error growth the variational equation for this case becomes y dot t is equal to Jacobian of f times y of t this y plays the role so instead of epsilon sorry instead of epsilon I am going to consider this as y naught so let y naught be the initial difference let y of t be the difference at time t epsilon be used in the context of a discrete time y of t I am going to use it in the case of continuous time so y of t is related to the rate of growth of y of t is given by this equation this equation is a linear dynamical system the linear dynamics varies along the trajectory again so this is the Jacobian if you discretize this in the interval 0 to t there are n sub intervals of time tau so n tau b t so this is 0 this is t I am going to divide it into intervals where the sub intervals are time of length tau and t is equal to n tap so now I am now going to assume that this equation which is given by 21 is such that during a given interval of time my d t of f remains constant that means the interval of time tau is so small that my Jacobian does not change too much within that small interval of time therefore I am now going to assume the d t of f in the in a small interval going from k tau to k plus 1 tau it remains constant and that constant value is lk so what is lk lk is the constant value of the Jacobian of the system during the time of evolution from k to k plus 1 or k tau to k plus 1 tau if tau is small this is a reasonably very good assumption to make therefore y dot therefore y dot d t of f y t by combining 22 so 22 when substituted in 21 becomes an equation like this sorry becomes an equation of this time that is a linear equation this equation depends on k you can see lk depends on this interval I can solve this equation very readily that is yk plus 1 is equal to e to the power of lk tau yk so if a time k this is time k plus 1 if yk is the initial condition I would like to be able to compute the solution at yk plus 1 the matrix lk remains the same so the the solutions at two intervals of at two end points are related by this equation 24 if I iterate this from y0 to yn I get this relation you can readily see by iterating 24 and using the definition of lk being constant in that in that domain yn is equal to e to the power of summation j is equal to 0 to n minus 1 lj times tau why not now from the definition of the Lyapunov index the Lyapunov index of time t I am assuming my my my my my forecast horizon is capital T so when my forecast horizon is is is capital T from the definition of 20 now it should be tau so this should be sorry as tau goes to 0 not not T so T is fixed tau goes to 0 as tau goes to 0 n goes to infinity n tau represent the total time capital T so n tau is always capital T please understand n tau is capital n tau is capital T so as tau goes to 0 n goes to infinity but the product is fixed T the product is fixed T now I am taking the logarithm of yn divided by y0 what is yn absolute value of yn is the error at time capital N y0 is the error at time 0 I am taking the ratio it is a logarithm of the ratio of the magnitude of the error at time n to time 0 so why is that the error is now I am interested in the in the in the errors that affect the prediction now if I if I if I substitute 25 the expression for yn in here and simplify what does it become it essentially becomes the average of lj that is l bar T that is a beautiful expression that comes from a simple so lambda T which is the function of T T is finite it becomes l bar T what is l bar T is the average value of the Jacobian along the trajectory that is beautiful that is beautiful so lambda T refers to the mean growth of errors during the interval 0 to T now what do I want to get the Lyapunov index let T go to infinity that is the definition of Lyapunov index so what is the Lyapunov index is the is the limit of capital is the limit of lambda T with the capital T going to infinity and that gives raise to that gives raise to the definition of average rate of growth of errors in unbounded intervals so what does that tell you if lambda is defined that way you can now see that yT at any time is equal to e to the power of lambda T times y0 clearly the error grows when lambda is greater than 0 so what does it mean one can analyze the quality of prediction by analyzing lambda if lambda is greater than 0 what does it mean like average rate of growth is greater than 0 the average rate of growth is greater than 0 error grows so in systems where lambda is greater than 0 there is an inherent predictability limit in systems where lambda is less than 0 there is no predictability limit I can predict for the whole future so that is the import of this analysis so I want you to go back and understand the theory reasonably well so we presented a general theory we are trying to illustrate the general theory based on a very simple scalar continuous time dynamics we have introduced the notion of average growth of errors so it is the average rate of growth of errors over asymptotically long intervals of time that helps to indicate the quality of prediction generated by the model so lambda depends on the model so lambda is going to be able to give us the guideline how to believe and how long I can believe the forecast generated by the model so how do we define the predictability limit the predictability limit Tp is given by Tp is equal to 1 over lambda so if lambda is positive my prediction can hold water only up to the forecast horizon Tp which is equal to 1 over lambda so what is 1 over lambda if y of t is equal to e to the power of lambda t y not if t is equal to Tp y not if Tp is equal to 1 over lambda y of t is equal to e times lambda y not that means the ratio of the magnitude of the initially ratio of the magnitude of the error time t to the initial error is the order of e this limit for atmospheric models used to be 1 or 2 days about 3 4 decades ago now has gone to the order of 5 to 7 days it is this improvement from 1 to 2 days to 5 to 7 days is the achievement by the meteorological community by bringing in the science of data simulation and analysis of analysis of predictions created the model so you can see the role of data simulation the role of errors in the initial condition the role of trying to measure the ratio of the errors at two different times and an understanding the rate at which the errors grow these are some of the beautiful of mechanisms using which one can create estimations of how good a forecast is that is the ultimate key so creating forecast is one thing trying to attest the goodness of the forecast is something else we have learned how to create forecast by so many different methods of data simulation the analysis of data simulation is not complete until we understand the predictability limit how long the predictability how long the prediction can hold water that is the idea here that is for time greater than tp tp for time greater than tp the error in the initial condition will overwhelm the signal corresponding to the base forecast clearly lambda is a function of x not so and the given value the parameter lambda varies as the function of the initial condition also the parameters now I am going to talk about how to compute lambda for real systems so we have seen the importance of a lack of index now I am going to talk about a renormalization strategy by which one can implement an algorithm on a any model to be able to compute the lack of index once a once a quantity is of great interest we need to be able to compute it so let us let us look at an algorithm now that epsilon not be the initial value of a perturbation based on the initial condition x not bar so initial condition is fixed given model is fixed alpha is fixed so please understand my model is xk plus 1 is equal to m of xk alpha so m is fixed alpha is fixed I have initial condition that is also fixed if I am interested in the ratio epsilon s to epsilon not sorry epsilon not epsilon s to epsilon not I can essentially express this a product like this so what is that epsilon s divided by epsilon s minus 1 times epsilon s minus 1 divided by epsilon s minus 2 likewise epsilon 1 divided by epsilon not all the cross terms will leaving behind this is equal to epsilon s divided by epsilon now therefore this ratio is equal to the product of the ratios between two successive times so what is that we are interested in we are interested in zero yes I am interested in trying to multiple of the product I am sorry I am interested in the ratios two successive times the product thereof the product thereof gives you the product from epsilon s to epsilon 0 so you can really see this product is given by this ratio so then from our definition what is lambda the Lyapunov index Lyapunov index is is the average of the log of the average of the log of the so I have to take the log of both sides if the log of the both sides log of the product with the sum of the logs so logarithm of ab is equal to log of a plus log of b so if I took the logarithm on both sides of 30 I get an expression so the ratio of the log of both sides as epsilon not goes to 0 so I would like my epsilon not to be as small as possible I want my n to go to infinity what is n the number of discrete intervals of time and that is the limit of the average of the logarithm of the amplification along the trajectory that is the basic idea so I would like to be able to compute this if I can compute 31 using a clever algorithm I am done and to be able to express 31 the key is in expressing the ratio of the error ratio the norms of errors the time yes to time 0 as the product of ratios that consecutive times so this is the product of ratios in the consecutive times the product of the ratios consecutive times the numerator denominator cancels leaving behind epsilon s and epsilon not so that is a very clever way of writing the ratio of the norms at time yes to time 0 so this is this is one of the ways of computing computing there so how do we how do we do that so here is an algorithm now so let epsilon epsilon not let epsilon I am sorry that is right let epsilon not be the initial error then run the model once if interval time compute the error epsilon 1 compute the ratio the norm of epsilon 1 by norm of epsilon 0 then run the model from epsilon 1 to epsilon 2 then multiply this by norm of epsilon 2 times epsilon 1 continue if I continued I get epsilon s divided by epsilon s minus 1 that is equal to norm of epsilon s divided by norm of epsilon 0 so by running the model once a better time by running the by computing the ratios by taking the product of those ratios taking the the logarithm of the ratio I can I can I can in a way estimate lambda so this provides an easy algorithm to be able to to be able to evaluate lambda so this is what is called renormalization strategy now I am going to give you well known information about certain facts about Lyapunov indices for forced okay I have not defined what a chaotic system now I let me let me let me talk about the definition of a chaotic system when do I say a dynamic system is chaotic if given a dynamical system if you compute the Lyapunov index if the Lyapunov index is positive in some part of the domain then it is called chaotic why it is chaotic in a chaotic system error grows on an average right if the error grows on an average rate then the predictability becomes very very very very difficult so what is the measure of the that predictability is Lyapunov index so Lyapunov index is being positive is one of the signatures of the chaotic system and also lambda greater than 0 is an indication of the fact that the system is extremely sensitive to initial condition because the initial condition errors go at a rate lambda so Lyapunov index sensitive to the initial condition being chaotic all these things are related concepts so given a forced chaotic system the first Lyapunov index is positive which is responsible for sensitivity the initial sensitivity responsible for the sensitivity the initial condition now I only talked about a Lyapunov index you may ask what is the what is the first Lyapunov index again that requires little bit of an explanation you remember we talked about the the the matrix let us go back in the definition of Lyapunov indices we talked about this matrix in 19 lambda m is a matrix a matrix have m different eigenvalues so I am now if I talk about the analysis of the first eigenvalue that leads to first Lyapunov index second eigenvalues leads to second Lyapunov index so if I system if if if the order of the system is is is 3 in general it should have 3 Lyapunov indices I only talked about the the the maximum value of the Lyapunov index or the Lyapunov index corresponding to the first eigenvalue so in n dimensional system has n eigenvalues n dimensional system can have n in Lyapunov indices so if there are n Lyapunov indices you can talk about the first second third and so on I am now going to talk about some of the well-known facts about Lyapunov indices for different types of systems. So for a dissipative system the sum of all the Lyapunov index must be negative that is a fact this can be proven it is proven on many good books on introduction to chaos theory for this for this class of system one of the intermediary Lyapunov indices is also 0 I am going to illustrate it further these are some of the well-known facts I am trying to summarize explaining each of this will take at least one or two lectures but here I am trying to collect many of the many of the results. The growth rate of the so what is the lambda 1 tells you so if I have a 3 dimensional system I have 3 Lyapunov indices lambda 1 tells you the growth rate of errors along one line lambda 1 plus lambda 2 refers to the growth rate of surface areas lambda 1 plus lambda 2 plus lambda 3 refers to growth rate of specific volumes that is where different Lyapunov indices come into play so likewise the sum of the first k Lyapunov indices tells you the rate of growth of k dimensional volumes so these are all very many simple facts. So now I am going to talk about some of the specific values for specific systems many of us are introduced the notion of Lorenz 1963 model it was Lorenz 1963 for the first time introduced or accidentally found the presence of chaotic system so for it is what there is this 1963 model of Lorenz is one of the most thoroughly analyzed understood systems of differential equation it consists of 3 ordinary differential equation coupled with nonlinearity it is also a dissipative system the because it is dissipative it is in general does not go to infinity it is the orbit is bounded the first Lyapunov index for this is 0.09 therefore e to the power of 0.9 is equal to 2.4596 so what does it mean any line segment any error along the line or any yeah any error along the line will magnify at the rate 2.4596 that is the average rate of growth of errors but the total Lyapunov indices for this is minus 11.9 that means the 3 dimensional volumes decrease at this rate decrease at this rate so what is the idea here the volumes decrease the attractor where the where the where the where the orbits of the system exist the you cannot predict where the orbits will be at what time because it is intrinsically chaotic because lambda 1 is positive. So that is the overall characteristics of what is called the Lorenz system and examples of values of the first Lyapunov index if the first or Lyapunov indices are positive then lambda prime p is equal to sum of all the Lyapunov indices which are positive two states there are infinitely close diverges the rate lambda e to the power of lambda bar p times t so this is again the average rate of growth this is an extension of y t being equal to e to the power of lambda d times y of t so in this case it is y t times e to the power of lambda p times t times y not where lambda p is the I am sorry lambda bar p sorry where the lambda bar p is equal to the sum of the first or some of the first or Lyapunov indices which are positive so the sum of the positive Lyapunov indices refers to the average rate of growth in this case the predictability limit is given by 1 over lambda to the 1 over lambda bar p now you can see if there are more than 1 Lyapunov index that is positive lambda bar p is larger if lambda p lambda bar p is larger 1 over lambda p bar is smaller therefore if a system exhibits larger number of positive Lyapunov indices in such system the predictability limit is much smaller if a system does not have any positive Lyapunov index then I can make predictions for a long periods of time predictions for the long periods of time so that is the ultimate in essence of the theory of Lyapunov index now I am going to give you these are exercise problems but I am going to quote you some of the values of the Lyapunov indices are some of the very well known mathematical models for the logistic model with the parameterization for the Lyapunov indexes 0.6931 that is positive so logistic models are chaotic because they are extremely sensitive to initial condition there is a model called Hinon model the Hinon model has the first Lyapunov index is a 2 I am sorry for the Hinon model I have the lambda 1 is 0.42 lambda 2 is minus 0.162 for the Lorenz's model we talked about first one is 0.9 second one 0 last one is this so these are some of the examples of simple systems there is a difference between the properties of the logistic model is a discrete type model that tries to capture some of the principles of population dynamics Hinon model is a mathematically oriented model it may not correspond to any particular physical system Lorenz's model how did he obtain the three differential equation he obtained it from starting with some vorticity equation in this case not vorticity equation I am sorry it is the heat transfer problem and from that he applied the spectral methods and he from that derived the three sets of equation whose Lyapunov indices are positive 0 and negative. So Lorenz's model theoretically for the first one to be able to exhibit the notion of Lyapunov index and the role of Lyapunov index in predictability studies here are some other properties of Lyapunov indices I am going to now create a comparison a comparison steady state behavior attractor Lyapunov index dimension of the attractor these are summaries are very well known result a steady steady a steady state for a dynamic system can be an equilibrium point what is an equilibrium point there is an attractor to which all the solutions come and settle down those steady state are called equilibrium point equilibrium point attractor is a point the Lyapunov indices are all negative therefore the dimension of the attractor 0 for a periodic orbit it is a cycle there is one Lyapunov index is 0 rest of them are negative the dimension of the attractor is 1 I know I have not talked about the dimension attractor much more elaborately but I believe I wanted to provide you a quick summary of the notion of predictability and some of the models for which the predictability can be answered. So what does it mean if a system exhibits an equilibrium predictability is complete if a system exhibits periodic orbit the predictability is complete I do not have to worry about predictability I can predict it for all for all time if the system exhibits a steady state with respect to periodic orbits the attractor is said to form a torus in which case two eigenvalues are 0 and the rest of the eigenvalues are negative the dimension of the attractor is 2 for a chaotic system the attractor is called fractal I have not introduced even the notion of fractal object fractal object is a complicated object in the case of a fractal object in the case of a chaotic attractor such as the Lorenzo's attractor the at least one Lyapunov index is positive the rest of them could be 0 or negative and the sum of all the overall Lyapunov indices totally together could be less than 0 that indicates the system is totally dissipate dissipate you and the dimension of the attractor is non non integer that means if they have what is called fractal dimension. So the notion of fractal attractor fractal dimension chaotic behaviour they are all interrelated with each other with that we conclude our discussion of deterministic predictability I would like you to pursue these of this exercise what is the exercise consider the simple logistic model this is a very good exercise consider alpha this is not alpha a the same parameter a as in here consider the range of parameter from 1 to 4 x0 is now is equal to 0 to 1 so I would like you to pick an x0 pick an alpha compute lambda. So plot this is a very good computing exercise Lyapunov index versus x0 for a given value of alpha sorry for a given value of for a given value of alpha so alpha is equal to 1 draw this curve alpha is equal to 2 draw this curve alpha is equal to 3 4 draw this curve you will see you will verify when alpha is 1 and 2 and 3 the system does not exhibit any sensitivity initial condition but when alpha is equal to 4 the system tends to exhibit extreme sensitivity initial condition and that is the model we talked about with this I conclude our discussion on predictability of deterministic system as a quick summary predictability analysis is of great interest after we assimilated data into the model once you assimilate data into the model I am going I have the capability to make prediction the question is how far for how long the prediction will hold water the answer to the question of how long the prediction will hold water depends on the intrinsic property of the model itself it relates to the sensitivity of the model solution with respect to the initial condition this sensitivity is essentially related to the forward sensitivity we have already defined so forward sensitivity analysis and analysis of Lyapunov index are intimately associated with each other even though I am not exploring that discussion in this in this in this in this set of lectures one way to be able to summarize the initial condition sensitivity is through a parameter called Lyapunov index Lyapunov index could be positive 0 or negative for dissipative system the sum of all the Lyapunov indexes must be less than 0 at least one Lyapunov index is 0 if one of the Lyapunov indices the greatest Lyapunov index is positive such systems are supposed to exhibit extreme sensitivity it turns out many many of the model that are currently being used in geophysical sciences to predict different types of geophysical phenomena have exhibited extreme sensitivity to initial condition that is the reason why we are not able to make long term predictions despite 50 plus years of progress in model building as well as data collection and data analysis and with all the advances we are able to increase the predictable limit to about 56 no more than 7 days these days so what does it mean the prediction problem continues it continues to dominate it continues to be a problem of great challenge and one of the goals of this study models data data simulation all relates to improving the predictability limit so what is the ultimate goal of doing all these things I would like to be able to make long range prediction very reliable long range prediction until such time we achieve the ability to make long range prediction our job is not done and there is no telling how long it may take though it all depends on it all depends on very many different aspects of very many different aspects of models data data simulation sensitivity to the model initial conditions so on and so forth so I have provided a simply a rudimentary ideas I simply scratch the surface on predictability is a very deep discipline I would encourage you hope this will encourage you to be able to look at some of the interesting books related to predictability analysis as well as chaos theory but as a starting point this exercise analyzing logistic model by choosing various alphas and choosing various x not and trying to compute will be a very good opening game in your understanding the theory of predictability thank you.