外文数值解论文A comparison of numerical solution
偏微分方程数值解PPT课件

从(1)得到:
u(ti)u(ti1)hu(ti)O(h)
精选
14
从(2)得到:
u(ti)u(ti1)hu(ti)O(h)
从(1)-(2)得到:
u(ti)u(ti1)2 hu(ti1)O (h2)
从(1)+(2)得到:
u (ti)u (ti 1) 2 u h (2 ti) u (ti 1 ) O (h 2)
精选
15
对经典的初值问题
du
dt
f (t,u )
u ( 0 ) u 0
t (0,T)
满足Lipschitz条件
4
常微分方程的数值解
大气科学中
常微分方程和偏微分方程的关系
1. 大气行星边界层(近地面具有湍流运动特性的大 气薄层,1—1.5km), 埃克曼(V.W.Ekman)(瑞典) 螺线的导出;
2. 1963年,美国气象学家Lorenz在研究热对流的 不稳定问题时,使用高截断的谱方法,由 Boussinesq流体的闭合方程组得到了一个完全确 定的三阶常微分方程组,即著名的Lorenz系统。
2. Curtis F.Gerald and Patrick O., Applied Numerical Analysis, Person Education, Inc., 2004.
3. Eugenia Kalnay, Atmospheric Modeling, Data Assimilation and Predictability, the press Syndicate of the University of Cambridge,2003.
ìïïïïïïïïïíïïïïïïïïî
x1¢=
x
¢
A Numerical Method for Solution of Ordinary Differential Equations of Fractional Order

(6)
In the general case the Riemann-Liouville fractional derivatives do not commute
β α+β β α α (a Dt (a Dt y ))(t) = (a Dt )(t) = (a Dt (a Dt y ))(t) .
(7)
Extending our considerations, the integer operator commutes with the fractional operator
β α−β α (a Dt (a It y ))(t) = (a Dt )(t) ,
(9)
but they do not commute in the opposite way. We turn our attention to the composition rule in two following ways. First of all we found in literature [10,11] the fact, that authors neglect the general property of fractional derivatives given by eqn. (7). Regarding to solution of fractional differential equations we will apply above properties in the next sections.
α y )(t) = (a It
Econometric and Statistical Computing Using Ox

Econometric and Statistical Computing Using OxFRANCISCO CRIBARI–NETO1and SPYROS G.ZARKOS21Departamento de Estat´ıstica,CCEN,Universidade Federal de Pernambuco,Recife/PE,50740–540,Brazil E-mail:cribari@npd.ufpe.br2National Bank of Greece,86Eolou str.,Athens10232,GreeceE-mail:s.zarkos@primeminister.grAbstract.This paper reviews the matrix programming language Ox from the viewpoint of an econometri-cian/statistician.We focus on scientific programming using Ox and discuss examples of possible interest to econometricians and statisticians,such as random number generation,maximum likelihood estimation,and Monte Carlo simulation.Ox is a remarkable matrix programming language which is well suited to research and teaching in econometrics and statistics.Key words:C programming language,graphics,matrix programming language,maximum likelihood estima-tion,Monte Carlo simulation,OxOne of the cultural barriers that separates computer scientists from regular scientists and engineers is a differing point of view on whether a30%or50%loss of speed is worth worrying about.In many real-time state-of-the art scientific applications,such a loss is catastrophic.The practical scientist is trying to solve tomorrow’s problem with yesterday’s computer;the computer scientist,we think, often has it the other way.Press et.al.(1992,p.25) 1.IntroductionApplied statisticians,econometricians and economists often need to write programs that implement estimation and testing procedures.With computers powerful and affordable as they are nowadays,they tend to do that in programming environments rather than in low level programming languages.The former(e.g.,GAUSS,MATLAB,R,S-PLUS)make programming accessible to the vast majority of researchers,and,in many cases,can be combined with the latter(e.g.,C,Fortran)to achieve additional gains in speed.The existence of pre-packaged routines in statistical software that is otherwise best suited to perform data analysis(such as in S-PLUS)does not make the need for“statistical comput-ing”any less urgent.Indeed,many newly developed techniques are not rapidly implemented into statistical software.If one wishes to use such techniques,he/she would have to program them.Additionally,several techniques are very computer-intensive,and require efficient pro-gramming environments/languages(e.g.,bootstrap within a Monte Carlo simulation,double bootstrap,etc.).It would be nearly impossible to perform such computer-intensive tasks with traditional statistical software.Finally,programming forces one to think harder about the problem at hand,the estimation and testing methods that he/she will choose to use.Of course,the most convincing argument may be the following quote from the late John Tukey:“In a world in which the price of calculation continues to decrease rapidly,but the price of theorem proving continues to hold steady or increase,elementary economics indicates that we ought to spend a larger fraction of our time on calculation.”1The focus of our paper is on the use of Ox for‘econometric computing’.That is,we discuss features of the Ox language that may be of interest to statisticians and econometricians,and exemplify their use through examples.Readers interested in reviews of Ox,including the language structure,its syntax,and its advantages and disadvantages,are referred to Cribari–Neto(1997),Keng and Orzag(1997),Kusters and Steffen(1996)and Podivinsky(1999).1 2.A Brief Overview of OxOx is a matrix programming language with object-oriented support developed by Jur-gen Doornik,a Dutch graduate student(at the time)at Nuffield College,Oxford.The development of Ox started in April1994.Doornik’s primary goal was to develop a matrix programming language for the simulations he wished to perform for his doctoral dissertation. The veryfirst preliminary version of Ox dates back to November1994.In the summer of 1995,two other econometricians at Nuffield College started using Ox for their research:Neil Shephard and Richard Spady.From that point on,the development of Ox became a serious affair.The current Ox version is numbered3.00.Ox binaries are available for Windows and severalflavors of UNIX(including Linux)and can be downloaded from /Users/Doornik/,which is the main Ox web page.All versions are free for educational purposes and academic research,with the exception of the‘Professional Windows version’.This commercial version comes with a nice interface for graphics known as GiveWin(available for purchase from Timberlake Consultants, ).The free Ox versions can be launched from the command line in a console/terminal win-dow,which explains why they are also known as‘console versions’.Doornik also distributes freely a powerful text editor for Windows:OxEdit(see also the OxEdit web page,which is currently at ).It can be used as a front-end not only to Ox(the console version)but also to other programs and languages,such as C,C++,T E X,L a T E X,etc.The Ox syntax is very similar to that of C,C++and Java.In fact,its similarity to C (at least as far as syntax goes)is one of its major advantages.2One characteristic similarity with C/C++is in the indexing,which starts at zero,and not at one.This means that thefirst element of a matrix,say A,is accessed as A[0][0]instead of as A[1][1].A key difference between Ox and languages such as C,C++and Java is that matrix is a basic type in Ox. Also,when programming in Ox one needs to declare the variables that will be used in the program(as is the case in C/C++),but unlike in C/C++,one does not have to specify the type of the variables that are declared.Ox’s most impressive feature is that it comes with a comprehensive mathematical and statistical function library.A number of useful functions and methods are implemented into the language,which makes it very useful for scientific 1A detailed comparison involving GAUSS,Macsyma,Maple,Mathematica,MATLAB,MuPAD,O-Matrix,Ox, R-Lab,Scilab,and S-PLUS can be found at http://www.scientificweb.de/ncrunch/ncrunch.pdf(“Com-parison of mathematical programs for data analysis”by Stefan Steinhaus).Ox is the winner when it comes to speed.2Other important advantages of Ox are the fact that it is fast,free,can be easily linked to C,Fortran, etc.,and can read and write data in several different formats(ASCII,Gauss,Excel,Stata,Lotus,PcGive, etc.).2programming.Ox comes with a comprehensive set of helpfiles in HTML form.The documentation of the language can be also found in Doornik(2001).A good introduction to Ox is Doornik, Draisma and Ooms(1998).3.A Few Simple IllustrationsOurfirst example is a very simple one,and intends to show the similarity between the Ox and C syntaxes.We wish to develop a program that produces a small table converting temperatures in Fahrenheit to Celsius(from0F to300F in steps of20F).The source of this example is Kerninghan and Ritchie(1988).The C code can be written as follows./****************************************************************PROGRAM:celsius.c**USAGE:To generate a conversion table of temperatures(from*Fahrenheit to Celsius).Based on an example in the*Kernighan&Ritchie’s book.****************************************************************/#include<stdio.h>int main(void){int fahr;printf("\nConversion table(F to C)\n\n");printf("\t%3s%5s\n","F","C");/*Loop over temperatures*/for(fahr=0;fahr<=300;fahr+=20){printf("\t%3d%6.1f\n",fahr, 5.0*(fahr-32)/9.0);}printf("\n");return0;}The output produced by compiled C code using the gcc compiler(Stallman,1999)under the Linux operating system(MacKinnon,1999)is:[cribari@edgeworth c]$gcc-O2-o celsius celsius.c[cribari@edgeworth c]$./celsiusConversion table(F to C)F C0-17.820-6.7340 4.46015.68026.710037.812048.914060.016071.118082.220093.3220104.4240115.6260126.7280137.8300148.9The next step is to write the same program in Ox code.The Ox transcription of the celcius.c program follows:/****************************************************************PROGRAM:celsius.ox**USAGE:To generate a conversion table of temperatures(from*Fahrenheit to Celsius).Based on an example in the*Kernighan&Ritchie’s book.***************************************************************/#include<oxstd.h>main(){decl fahr;print("\nConversion table(F to C)\n\n");print("\t F C\n");//Loop over temperaturesfor(fahr=0;fahr<=300;fahr+=20){print("\t","%3d",fahr);print("","%6.1f", 5.0*(fahr-32)/9.0,"\n");}print("\n");}The Ox output is:[cribari@edgeworth programs]$oxl celsiusOx version 3.00(Linux)(C)J.A.Doornik,1994-2001Conversion table(F to C)F C40-17.820-6.740 4.46015.68026.710037.812048.914060.016071.118082.220093.3220104.4240115.6260126.7280137.8300148.9The two programs above show that the Ox and C syntaxes are indeed very similar.Note that Ox accepts C style comments(/*...*/),and also C++like comments to the end of the line(//).3We also note that,unlike C,Ox accepts nested comments.The similarity between the Ox and C syntaxes is a major advantage of Ox over other matrix languages.Kendrick and Amman(1999)provide an overview of programming languages in economics.In the introduction of their paper,they give the following advice to users who are starting to program:“Begin with one of the high-level or modeling languages.(...)Then work downward in the chain and learn either Fortran,C,C++,or Java.”If a user then starts with Ox and‘works downward’to C or C++the transition will be smoother than if he/she starts the chain with other high level languages.As a second illustration of the use of Ox in econometrics and statistics,we develop a simple program thatfirst simulates a large number of coin tosses,and then counts the frequency (percentage)of tails.The code which is an Ox translation,with a smaller total number of runs,of the C code given in Cribari–Neto(1999),thus illustrates Kolmogorov’s Law of Large Numbers.We begin by writing a loop-based version of the coin tossing experiment./*******************************************************************PROGRAM:coin_loop.ox**USE:Simulates a large number of coin tosses and prints*the percentage of tails.**PURPOSE:The program illustrates the first version of the*law of large numbers which dates back to James*Bernoulli.******************************************************************/#include<oxstd.h>/*maximum number of coin tosses*/3Ox also borrows from Java;the println function,for instance,comes from the Java programming language.5const decl COIN_MAX=1000000;main(){decl j,dExecTime,temp,result,tail,s;//Start the clock(to time the execution of the program).dExecTime=timer();//Choose the random number generator.ranseed("GM");//Main loop:for(j=10;j<=COIN_MAX;j*=10){tail=0;for(s=0;s<j;s++){temp=ranu(1,1);tail=temp>0.5?tail:tail+1;}result=100.0*tail/j;print("Percentage of tails from",j,"tosses:","%8.2f",result,"\n");}print("\nEXECUTION TIME:",timespan(dExecTime),"\n");}The instruction tail=temp>0.5?tail:tail+1;does exactly what it does in C: it sets the variable tail equal to itself if the stated condition is true(temp>0.5)and to tail+1otherwise.We now vectorize the above code for speed.The motivation is obvious:vectorization usually leads to efficiency gains,unless of course one runs into memory problems.It is note-worthy that one of the main differences between a matrix programming language and a low level language,such as C and C++,is that programs should exploit vector and matrix opera-tions when written and executed in a matrix-oriented language,such as Ox.The vectorized code for the example at hand is:/*******************************************************************PROGRAM:coin_vec.ox**USE:Simulates a large number of coin tosses and prints*the percentage of tails.**PURPOSE:The program illustrates the first version of the*law of large numbers which dates back to James*Bernoulli.******************************************************************/6#include<oxstd.h>/*maximum number of coin tosses*/const decl COIN_MAX=1000000;main(){decl j,dExecTime,temp,tail;//Start the clock(to time the execution of the program).dExecTime=timer();//Choose the random number generator.ranseed("GM");//Coin tossing:for(j=10;j<=COIN_MAX;j*=10){temp=ranu(1,j);tail=sumr(temp.<0.5)*(100.0/j);print("Percentage of tails from",j,"tosses:","%8.2f",double(tail),"\n");}print("\nEXECUTION TIME:",timespan(dExecTime),"\n");}The output of the loop-based program is:[cribari@edgeworth programs]$oxl coin_loopOx version 3.00(Linux)(C)J.A.Doornik,1994-2001Percentage of tails from10tosses:40.00Percentage of tails from100tosses:53.00Percentage of tails from1000tosses:49.10Percentage of tails from10000tosses:49.69Percentage of tails from100000tosses:49.83Percentage of tails from1000000tosses:49.99EXECUTION TIME: 2.41whereas the vectorized code generates the following output: [cribari@edgeworth programs]$oxl coin_vecOx version 3.00(Linux)(C)J.A.Doornik,1994-2001Percentage of tails from10tosses:40.00Percentage of tails from100tosses:53.00Percentage of tails from1000tosses:49.10Percentage of tails from10000tosses:49.69Percentage of tails from100000tosses:49.83Percentage of tails from1000000tosses:49.99EXECUTION TIME:0.237Note that the empirical frequency of tails approaches1/2,the population mean,as predicted by the Law of Large Numbers.As far as efficiency goes,we see that vectorization leads to a sizeable improvement.The loop-based program yields an execution time which is over10 times greater than that of its vectorized version,on a DELL Pentium III1GHz computer with512MB RAM running on Linux.4Some languages,like C,operate faster on rows than on columns.The same logic applies to Ox.To illustrate the claim,we modify the vectorized code so that the random draws are stored in a column vector(they were previously stored in a row vector).To that end,one only needs to change two lines of code:for(j=10;j<=COIN_MAX;j*=10){temp=ranu(j,1);//1st changetail=sumc(temp.<0.5)*(100.0/j);//2nd changeprint("Percentage of tails from",j,"tosses:","%8.2f",double(tail),"\n");}This new vectorized code now runs in0.35second.That is,we see a speed penalty of over 50%when we transpose the code so that we work with a large column vector instead of working with a large row vector.4.Econometric ApplicationsMaximum likelihood estimates oftentimes need to be computed using a nonlinear op-timization scheme.In order to illustrate how that can be done using Ox,we consider the maximum likelihood estimation of the number of degrees-of-freedom of a Student t distri-bution.Maximization is performed using a quasi-Newton method(known as the‘BFGS’method)with numerical gradient,i.e.,without specifying the score function.(Note that this estimator is substantially biased in small samples.)It is noteworthy that Ox has routines for other optimization methods as well,such as the Newton-Raphson and the BHHH methods. An advantage of the BFGS method is that it allows users to maximize likelihoods without having to specify a score function.See Press et al.(1992,Chapter10)for details on the BFGS and other nonlinear optimization methods.See also Mittelhammer,Judge and Miller(2000,§8.13),who on page199write that“[t]he BFGS algorithm is generally regarded as the best performing method.”The example below uses a random sample of size50,the true value of the parameter is3,and the initial value of the optimization scheme is2.(We have neglected a constant in the log-likelihood function.)/**************************************************************PROGRAM:t.ox**USAGE:Maximum likelihood estimation of the number of*degrees of freedom of a Student t distribution.*************************************************************/4The operating system was Mandrake Linux8.0running on kernel2.4.3.8#include<oxstd.h>#include<oxprob.h>#import<maximize>const decl N=50;static decl s_vx;fLogLik(const vP,const adFunc,const avScore,const amHess) {decl vone=ones(1,N);decl nu=vP[0];adFunc[0]=double(N*loggamma((nu+1)/2)-(N/2)*log(nu)-N*loggamma(nu/2)-((nu+1)/2)*(vone*log(1+(s_vx.^2)/nu)));if(isnan(adFunc[0])||isdotinf(adFunc[0]))return0;elsereturn1;//1indicates success}main(){decl vp,dfunc,dnu,ir;ranseed("GM");vp= 2.0;dnu= 3.0;s_vx=rant(N,1,3);ir=MaxBFGS(fLogLik,&vp,&dfunc,0,TRUE);print("\nCONVERGENCE:",MaxConvergenceMsg(ir));print("\nMaximized log-likelihood:","%7.3f",dfunc);print("\nTrue value of nu:","%6.3f",dnu);print("\nML estimate of nu:","%6.3f",double(vp));print("\nSample size:","%6d",N);print("\n");}Here is the Ox output:[cribari@edgeworth programs]$oxl tOx version 3.00(Linux)(C)J.A.Doornik,1994-2001CONVERGENCE:Strong convergenceMaximized log-likelihood:-72.813True value of nu: 3.0009ML estimate of nu: 1.566Sample size:50The maximum likelihood estimate ofν,whose true value is3,is ν=1.566.This example shows that nonlinear maximization of functions can be done with ease using Ox.Of course, one can estimate more complex models in a similar fashion.For example,the parameters of a nonlinear regression model can be estimated by setting up a log-likelihood function,and maximizing it with a MaxBFGS call.It is important to note,however,that Ox does not come with routines for performing constrained maximization.The inclusion of such functions in Ox would be a great addition to the language.A number of people have developed add-on packages for Ox.These handle dynamic panel data(DPD),ARFIMA models,conditionally heteroskedastic models,stochastic volatil-ity models,state space forms.There is,moreover,Ox code for quantile regressions,and in particular, 1(i.e.,least absolute deviations)regressions.The code corresponds to the al-gorithm described in Portnoy and Koenker(1997)and is available at Roger Koenker’s web page(/roger/research/rqn/rqn.html).We consider,next,the G@RCH2.0package recently developed by S´e bastien Laurent and Jean–Philippe Peters,which is dedicated to the estimation and forecasting of ARCH,GARCH models.The GARCH add-on package comes in two versions,namely:(i)the‘Full Version’which requires a registered copy of Ox Professional3.00,since it is launched from OxPack and makes use of the GiveWin interface,and(ii)the‘Light Version’which only requires the free (‘console’)version of Ox.It relies on Ox’s object-oriented programming capabilities,being a derived class of Ox’s Modelbase type of class.The package is available for download at http://www.egss.ulg.ac.be/garch.We borrow the example program(GarchEstim.ox)in order to illustrate the use of the GARCH code(as with everything else,in the context of the console,i.e.free,version of Ox).The GARCH object(which is created with the source code provided by this add-on package)allows for the estimation of a large number of uni-variate ARCH-type models(e.g.,ARCH,GARCH,IGARCH,FIGARCH,GJR,EGARCH, APARCH,FIEGARCH,FIAPARCH)under Gaussian,Student–t,skewed Student and gen-eralized error distributions.Forecasts(one-step-ahead density forecasts)of the conditional mean and variance are also available,as well as several misspecification tests and graphics commands.#include<oxstd.h>#import<packages/garch/garch>main(){decl garchobj;garchobj=new Garch();//***DATA***//garchobj.Load("Data/demsel.in7");();garchobj.Select(Y_VAR,{"DEM",0,0});10garchobj.SetSelSample(-1,1,-1,1);//***SPECIFICATIONS***//garchobj.CSTS(1,1);//cst in Mean(1or0),cst in Variance(1or0)garchobj.DISTRI(0);//0for Gauss,1for Student,2for GED,3for Skewed-Student garchobj.ARMA(0,0);//AR order(p),MA order(q).garchobj.ARFIMA(0);//1if Arfima wanted,0otherwisegarchobj.GARCH(1,1);//p order,q ordergarchobj.FIGARCH(0,0,1000);//Arg.1:1if Fractionnal Integration wanted.//Arg.2:0->BBM,1->Chung//Arg.3:if BBM,Truncation ordergarchobj.IGARCH(0);//1if IGARCH wanted,0otherwisegarchobj.EGARCH(0);//1if EGARCH wanted,0otherwisegarchobj.GJR(0);//1if GJR wanted,0otherwisegarchobj.APARCH(0);//1if APARCH wanted,0otherwise//***TESTS&FORECASTS***//garchobj.BOXPIERCE(<5;10;20>);//Lags for the Box-Pierce Q-statistics.garchobj.ARCHLAGS(<2;5;10>);//Lags for Engle’s LM ARCH test.garchobj.NYBLOM(1);//1to compute the Nyblom stability test,0otherwisegarchobj.PEARSON(<40;50;60>);//Cells(<40;50;60>)for the adjusted Pearson//Chi-square Goodness-of-fit test,0if not computed//G@RCH1.12garchobj.FORECAST(0,100);//Arg.1:1to launch the forecasting procedure,//0elsewhere//Arg.2:Number of one-step ahead forecasts//***OUTPUT***//garchobj.MLE(1);//0:both,1:MLE,2:QMLEgarchobj.COVAR(0);//if1,prints variance-covariance matrix of the parameters.garchobj.ITER(0);//Interval of iterations between printed intermediary results//(if no intermediary results wanted,enter’0’) garchobj.TESTSONLY(0);//if1,runs tests for the raw Y series,prior to//any estimation.garchobj.GRAPHS(0);//if1,prints graphics of the estimations//(only when using GiveWin).garchobj.FOREGRAPHS(0);//if1,prints graphics of the forecasts//(only when using GiveWin).//***PARAMETERS***//garchobj.BOUNDS(0);//1if bounded parameters wanted,0otherwisegarchobj.DoEstimation(<>);garchobj.STORE(0,0,0,0,0,"01",0);//Arg.1,2,3,4,5:if1->stored.(Res-SqRes-CondV-MeanFor-VarFor)//Arg.6:Suffix.The name of the saved series will be"Res_ARG6"//(or"MeanFor_ARG6",...).//Arg.7:if0,saves as an Excel spreadsheet(.xls).//If1,saves as a GiveWin dataset(.in7)delete garchobj;}11We have run the above code to obtain the MLE and QMLE results of an ARMA(0,0)model in the mean equation and GARCH(1,1)model in the variance equation,assuming Gaussian distributed errors.Some portmanteau tests,such as the Box–Pierce Q-statistic and the LM ARCH test,the Jarque–Bera normality test etc,were also calculated for the daily observations on the Dow Jones Industrial Average(Jan.1982-Dec.1999,a total of4,551observations). The output follows.Ox version 3.00(Linux)(C)J.A.Doornik,1994-2001Copyright for this package:urent and J.P.Peters,2000,2001.G@RCH package version 2.00,object created on14-08-2001----Database information----Sample:1-4313(4313observations)Frequency:1Variables:4Variable#obs#miss min mean max std.devDEM43130-6.3153-0.0022999 3.90740.75333PREC4313000.4259250.82935SUCC4313000.418550.81568OBSVAR43130 3.3897e-060.567539.853 1.3569 **********************SPECIFICATIONS*********************Mean Equation:ARMA(0,0)model.No regressor in the meanVariance Equation:GARCH(1,1)model.No regressor in the varianceThe distribution is a Gauss distribution.Strong convergence using numerical derivativesLog-likelihood=-4651.57Please wait:Computing the Std Errors...Maximum Likelihood EstimationCoefficient Std.Error t-value t-probCst(M)0.0031860.0100190.31800.7505Cst(V)0.0178730.003216 5.5580.0000GARCH(Beta1)0.8702150.01168674.460.0000ARCH(Alpha1)0.1028470.00964210.670.0000Estimated Parameters Vector:0.003186;0.017873;0.870215;0.102847No.Observations:4313No.Parameters:4*************TESTS**12***********Statistic t-Test P-ValueSkewness-0.20031 5.37237.7733e-08Excess Kurtosis 1.868425.061 1.3133e-138Jarque-Bera656.19656.19 3.2440e-143---------------Information Criterium(minimize)Akaike 2.158856Shibata 2.158855Schwarz 2.164763Hannan-Quinn 2.160942---------------BOX-PIERCE:ValueMean of standardized residuals-0.00065Mean of squared standardized residuals0.99808H0:No serial correlation==>Accept H0when prob.is High[Q<Chisq(lag)] Box-Pierce Q-statistics on residualsQ(5)=17.7914[0.00321948]Q(10)=26.4749[0.00315138]Q(20)=44.9781[0.00111103]Box-Pierce Q-statistics on squared residuals-->P-values adjusted by2degree(s)of freedomQ(5)=8.01956[0.0456093]Q(10)=12.4119[0.133749]Q(20)=34.563[0.0107229]--------------ARCH1-2test:F(2,4306)= 2.7378[0.0648]ARCH1-5test:F(5,4300)= 1.5635[0.1668]ARCH1-10test:F(10,4290)= 1.2342[0.2632]--------------Diagnostic test based on the news impact curve(EGARCH vs.GARCH)Test ProbSign Bias t-Test 1.175980.23960Negative Size Bias t-Test 1.828560.06747Positive Size Bias t-Test0.975420.32935Joint Test for the Three Effects 4.468820.21509---------------Joint Statistic of the Nyblom test of stability: 1.77507Individual Nyblom Statistics:Cst(M)0.43501Cst(V)0.22234GARCH(Beta1)0.10147ARCH(Alpha1)0.10050Rem:Asymptotic1%critical value for individual statistics=0.75.Asymptotic5%critical value for individual statistics=0.47.---------------Adjusted Pearson Chi-square Goodness-of-fit testLags Statistic P-Value(lag-1)P-Value(lag-k-1)4078.06890.0002040.0000405089.05190.0004090.00010060103.25320.0003250.00008913Rem.:k=#estimated parameters---------------Elapsed Time: 4.67seconds(or0.0778333minutes).The stochastic volatility package(SvPack),written by Neil Shephard,is essentially a dy-namic link library for Ox of C code that deals with the implementation of likelihood inference in volatility models.The fact that it is written in C guarantees optimal speed,whereas the linking to Ox definitely improves usability.It requires the Ox state space package(SsfPack), which provides for Kalmanfiltering,smoothing and simulation smoothing algorithms of Gaus-sian multivariate state space forms(see Koopman,Shephard and Doornik,1999;Ooms,1999, and also ),as well as ARMS(Adaptive Rejection Metropolis Sam-pling),an Ox front-end for C code for adaptive rejection sampling algorithms(i.e.,routines for efficient sampling from complicated univariate densities)developed and documented by Michael Pitt(based on C code by Wally Gilks).The Arfima package is a set of Ox functions that create a class(an ARFIMA object) for the estimation and testing of AR(F)IMA models(Beran,1994).The models can be esti-mated via exact maximum likelihood,modified profile likelihood and nonlinear least squares. ArfimaSim is an additional simulation class included in the Arfima package that provides the means for Monte Carlo experiments based on the Arfima class.The Dynamic Panel Data package,DPD,like the Arfima and G@RCH packages,is a nice example of object-oriented Ox programming.They are derived classes written in Ox.DPD, which is entirely written in Ox,implements dynamic panel data models,as well as some static ones,and can handle unbalanced panels.Monte Carlo experimentation is possible with the simulation class DPSSim,included in this Ox add-on package.5.GraphicsOx has a number of commands that help create publication-quality graphics.This is, however,one of the areas where more progress is expected.The graphics capabilities of the console version of Ox are not comparable to those of,say,GAUSS,MATLAB,R or S-PLUS.It is important to note,however,that the professional version of Ox comes with an impressive interface for graphics:GiveWin.It allows users,for example,to modify a graph with a few clicks of the mouse.With GiveWin,it is possible to edit all graphs on the screen, manipulate areas,add Greek letters,add labels,change fonts,etc.Therefore,users who intend to make extensive use of the plotting capabilities of the language to produce publication quality graphics should consider using the professional version of Ox.5An alternative strategy would be to use Ox for programming,save the results to afile, read the resultsfile into R,which is also free,and then produce publication quality plots from there.6It is also possible to use GnuDraw,an Ox package written by Charles Bos (http://www2.tinbergen.nl/~cbos/).GnuDraw allows users to create gnuplot(http:// )graphics from Ox,extending the possibilities offered by existing OxDraw 5The newest,just released,version3.00of Ox has improved graphics capabilities.For instance,it now has built-in functions for producing3D plots.6For details on R,see .14。
数值分析论文范文

数值分析论文范文Title: An Overview of Numerical AnalysisIntroduction:Numerical analysis is a field of mathematics that deals with the development and application of algorithms to solve mathematical problems. In this paper, we will provide an overview of numerical analysis, including its history, important concepts, and applications in various fields.History of Numerical Analysis:Important Concepts in Numerical Analysis:2. Error Analysis: Since numerical methods involve approximation, it is essential to quantify and analyze theerrors associated with these approximations. Error analysis provides insights into the accuracy and efficiency of numerical algorithms. Different error measures, such as absolute error, relative error, and truncation error, are used to evaluate the quality of the approximate solutions.3. Numerical Algorithms: Numerical analysis relies on the development and implementation of efficient algorithms to solve mathematical problems. Iterative methods, such as Newton's method for finding roots of equations and the Jacobi method for solving linear systems, are widely used in numerical analysis.Applications of Numerical Analysis:3. Finance: In finance, numerical analysis is used to model and solve problems related to option pricing, portfolio optimization, risk management, and financial derivatives pricing. The Black-Scholes-Merton model, for instance, relies heavily on numerical methods for pricing options.Conclusion:。
New Newton's Method with Third-order Convergence for Solving Nonlinear Equations

New Newton’s Method with Third-order Convergence for Solving Nonlinear EquationsOsama Yusuf AbabnehAbstract—For the last years,the variants of the Newton’s methodwith cubic convergence have become popular iterative methods tofind approximate solutions to the roots of non-linear equations.Thesemethods both enjoy cubic convergence at simple roots and do notrequire the evaluation of second order derivatives.In this paper,wepresent a new Newton’s method based on contra harmonic mean withcubically convergent.Numerical examples show that the new methodcan compete with the classical Newton’s method.Keywords—Third-order convergence,Non-linear equations,Root-finding,Iterative method.I.I NTRODUCTIONS OLVING non-linear equations is one of the most impor-tant problems in numerical analysis.In this paper,weconsider iterative methods tofind a simple root of a non-linear equation f(x)=0,where f:D⊂R→R for anopen interval D is a scalar function.The classical Newtonmethod for a single non-linear equation is written asx n+1=x n−f(x n)(x n).(1)This is an important and basic method[8],which converges quadratically.Recently,some modified Newton methods with cubic convergence have been developed in[1],[2],[3],[4],[5], [6]and[7].Here,we will obtain a new modification of New-tons method.Analysis of convergence shows the new method is cubically convergent.Its practical utility is demonstrated by numerical examples.Letαbe a simple zero of a sufficiently differentiable function f and consider the numerical solution of the equation f(x)=0.It is clear thatf(x)=f(x n)+xx nf (t)dt.(2)Suppose we interpolate f on the interval[x n,x]by the con-stant f (x n),then(x−x n)f (x n)provides an approximation for the indefinite integral in(2)and by taking x=αwe obtain 0≈f(x n)+(α−x n)f (x n).Thus,a new approximation x n+1toαis given byx n+1=x n−f(x n) f (x n).Dr.Osama Yusuf Ababneh is with the Department of Mathematics,Irbid National University,Irbid,Jordan e-mail:(ababnehukm@).On the other hand,if we approximate the indefinite integral in(2)by the trapezoidal rule and take x=α,we obtain 0≈f(x n)+12(α−x n)(f (x n)+f (α)),and therefore,a new approximation x n+1toαis given byx n+1=x n−2f(x n)f (x n)+f (x n+1).If the Newton’s method is used on the right-hand side of the above equation to overcome the implicity problem,thenx n+1=x n−2f(x n)f (x n)+f (z n+1),(3) wherez n+1=x n−f(x n)f (x n)is obtained which is,for n=0,1,2,...,the trapezoidal Newton’s method of Fernando et al.[1].Let us rewrite equation(3)asx n+1=x n−f(x n)(f (x n)+f (z n+1))/2,n=0,1, (4)So,this variant of Newton’s method can be viewed as obtained by using arithmetic mean of f (x n)and f (z n+1)instead of f (x n)in Newton’s method defined by(1).Therefore,we call it arithmetic mean Newton’s(AN)method.In[3],the harmonic mean instead of the arithmetic mean is used to get a new formulax n+1=x n−f(x n)(f (x n)+f (z n+1))2f (x n)f (z n+1),n=0,1, (5)which is called harmonic mean Newton’s(HN)method and used the midpoint to getx n+1=x n−f(x n)((x n+z n+1)/2),n=0,1, (6)which is called midpoint Newton’s(MN)method.II.N EW ITERATIVE METHOD AND CONVERGENCEANALYSISIf we use the contra harmonic mean instead of the arithmetic mean in(4),we get new Newton methodx n+1=x n−f(x n)(f (x n)+f (z n+1))f 2(x n)+f 2(z n+1),n=0,1, (7)wherez n +1=x n −f (x n )f (x n ),n =0,1, (8)we call contra harmonic Newton’s (CHN)method.Theorem 2.1:Let α∈D be a simple zero of a sufficientlydifferentiable function f :D ⊂R →R for an open interval D.If x 0is sufficiently close to α,then the methods defined by (7)converge cubically.Proof Let αbe a simple zero of f .Since f is sufficiently differentiable,by expanding f (x n )and f (x n )about αwe getf (x n )=f (α) e n +c 2e 2n +c 3e 3n+... ,(9)andf (x n )=f (α) 1+2c 2e n +3c 3e 2n +4c 4e 3n+... ,(10)where c k =(1/k !)f (k )(α)/f(α),k =2,3,...and e n =x n −α.Direct division gives usf (x n )f (x n )=e n −c 2e 2n +2(c 22−c 3)e 3n +O (e 4n ),and hence,for z n +1given in (8)we havez n +1=α+c 2e 2n +2(c 3−c 22)e 3n +O (e 4n ).(11)Again expanding f (z n +1)about αand using (11)we obtainf (z n +1)=f (α)+(z n +1−α)f (α)+(z n +1−α)2!f(α)+...=f (α)+[c 2e 2n +2(c 3−c 22)e 3n +O (e 4n )]f(α)+O (e 4n )=f (α)[1+2c 22e 2n +4(c 2c 3−c 32)e 3n +O (e 4n )].(12)By using (10)we obtainf 2(x n )=f 2(α)[1+4c 2e n + 4c 22+6c 3 e 2+(12c 2c 3+8c 4)e 3+...].From (12),we getf 2(z n +1)=f 2(α) 1+4c 22e 2n + 8c 2c 3−8c 32 e 3+... ,andf 2(x n )+f 2(z n +1)=2f 2(α)[1+2c 2e n + 4c 22+3c 3 e 2n + 4c 4+10c 2c 3−4c 32 e 3n +...].From (10)and (12)we getf (x n )+f (z n +1)=2f (α)[1+c 2e n +c 22+32c 3e 2n +2 c 2c 3−c 32+c 4 e 3n +O (e 4n )],and using (9)to getf (x n )(f (x n )+f (z n +1))=2f 2(α)[e n +2c 2e 2n+2c 22+52c 3e 3n +O (e 4n )].Hence,f (x n )(f (x n )+f (z n +1))f 2(x n )+f 2(z n +1)=e n −2c 22+12c 3 e 3n +O (e 4n ),x n +1=x n −f (x n )(f (x n )+f (z n +1))f 2(x n )+f 2(z n +1),x n +1=x n −e n −2c 22+12c 3 e 3n +O (e 4n ) ,or subtracting αfrom both sides of this equation we gete n +1=2c 22+12c 3 e 3n +O (e 4n ),which shows that contra harmonic Newton’s method is ofthird order.III.N UMERICAL RESULTS AND CONCLUSIONSIn this section,we present the results of some numerical tests to compare the efficiencies of the new method (CHN).We employed (CN)method,(AN)method of Fernando et al.[1],and (HN)and (MN)methods in [3].Numerical computations reported here have been carried out in a MTHEMATICA environment .The stopping criterion has been taken as |x n +1−x n |<ε,We used the fixed stopping criterion ε=10−14and the following test functions have been used.f 1(x )=x 3+4x 2−10,α=1.365230013414097,f 2(x )=x 2−e x −3x +2,α=0.2575302854398608,f 3(x )=Sinx 2−x 2+1,α=1.404491648215341,f 4(x )=Cosx −x,α=0.7390851332151607,f 5(x )=(x −1)3−1,α=2.In Table 1and Table 2,we give the number of iterations (N)and total number of function evaluations (TNFE)required to satisfy the stopping criterion.As far as the numerical results are considered,for most of the cases HN and MN methods requires the least number of function evaluations.All numerical results are in accordance with the theory and the basic advantage of the variants of Newton’s method based on means or integration methods that they do not require the computation of second-or higher-order derivatives although they are of third order.R EFERENCES[1]S.Weerakoon,T.G.I.Fernando,A variant of Newtons method withaccelerated third-order convergence,Appl.Math.Lett.13(2000)87-93.[2]M.Frontini,E.Sormani,Some variants of Newtons method with third-order convergence,put.140(2003)419-426.[3] A.Y .¨o zban,Some new variants of Newtons method,Appl.Math.Lett.17(2004)677-682.[4]M.Frontini,E.Sormani,Modified Newtons method with third-orderconvergence and multiple roots,put.Appl.Math.156(2003)345-354.TABLE II TERATION NUMBER(N)f(x)NCN AN HN MN CHNf1,x0=164445f2,x0=154444f2,x0=265545f2,x0=375555f3,x0=175455f3,x0=375445f4,x0=153444f4,x0=1.754444f4,x0=−0.364555f5,x0=375555TABLE IIT HE TOTAL NUMBER OF FUNCTION EVALUATIONS(TNFE)f(x)TNEFCN AN HN MN CHNf1,x0=11212121215f2,x0=11012121212f2,x0=21215151515f2,x0=31415151515f3,x0=11415121515f3,x0=31415121215f4,x0=1109121212f4,x0=1.71012121212f4,x0=−0.31212151515f5,x0=31415151515[5]Changbum Chun,A two-parameter third-order family of methods forsolving nonlinear equations,Applied Mathematics and Computation189 (2007)1822-1827.[6]Kou Jishenga,LiYitianb,Wang Xiuhuac,Third-order modification ofNewtons method,Journal of Computational and Applied Mathematics 205(2007)1 5.[7]Mamta,V.Kanwar,V.K.Kukreja,Sukhjit Singh,On some third-orderiterative methods for solving nonlinear equations,Applied Mathematics and Computation171(2005)272-280.[8] A.M.Ostrowski,Solution of Equations in Euclidean and Banach Space,third ed.,Academic Press,NewYork,1973.。
mapreduce数据分析-文档资料

ABSTRACT:There is currently considerable enthusiasm around the MapReduce (MR) paradigm for large-scale data analysis. Although the basic control flow of this framework has existed in parallel SQL database management systems (DBMS) for over 20 years, some have called MR a dramatically new computing model. In this paper, we describe and compare both paradigms. Furthermore, we evaluate both kinds of systems in terms of performance and development complexity. To this end, we define a benchmark consisting of a collection of tasks that we have run on an open source version of MR as well as on two parallel DBMSs. For each task, we measure each system’s performance for various degrees of parallelism on a cluster of 100 nodes. Our results reveal some interesting trade-offs. Although the process to load data into and tune the execution of parallel DBMSs took much longer than the MR system, the observed performance of these DBMSs was strikingly better. We speculate about the causes of the dramatic performance difference and consider implementation concepts that future systems should take from both kinds of architectures.
数学英文论文

070451 Controlling chaos based on an adaptive nonlinear compensatingmechanism*Corresponding author,Xu Shu ,email:123456789@Abstract The control problems of chaotic systems are investigated in the presence of parametric u ncertainty and persistent external distu rbances based on nonlinear control theory. B y designing a nonlinear compensating mechanism, the system deterministic nonlinearity, parametric uncertainty and disturbance effect can be compensated effectively. The renowned chaotic Lorenz system subject to parametric variations and external disturbances is studied as an illustrative example. From Lyapu nov stability theory, sufficient conditions for the choice of control parameters are derived to guarantee chaos control. Several groups of experiments are carried out, including parameter change experiments, set-point change experiments and disturbance experiments. Simulation results indicate that the chaotic motion can be regulated not only to stead y states but also to any desired periodic orbits with great immunity to parametric variations and external distu rbances.Keywords: chaotic system, nonlinear compensating mechanism, Lorenz chaotic systemPACC: 05451. IntroductionChaotic motion, as the peculiar behavior in deterministic systems, may be undesirable in many cases, so suppressing such a phenomenon has been intensively studied in recent years. Generally speaking chaos suppression and chaos synchronization[1-4 ]are two active research fields in chaos control and are both crucial in application of chaos. In the following letters we only deal with the problem of chaos suppression and will not discuss the chaos synchronization problem.Since the early 1990s, the small time-dependent parameter perturbation was introduced by Ott,Grebogi, and Y orke to eliminate chaos,[5]many effective control methods have been reported in various scientific literatures.[1-4,6-36,38-44,46] There are two lines in these methods. One is to introduce parameter perturbations to an accessible system parameter, [5-6,8-13] the other is to introduce an additive external force to the original uncontrolled chaotic system. [14-37,39-43,47] Along the first line, when system parameters are not accessible or can not be changed easily, or the environment perturbations are not avoided, these methods fail. Recently, using additive external force to achieve chaos suppression purpose is in the ascendant. Referring to the second line of the approaches, various techniques and methods have been proposed to achieve chaos elimination, to mention only a few:(ⅰ) linear state feedback controlIn Ref.[14] a conventional feedback controller was designed to drive the chaotic Duffing equation to one of its inherent multiperiodic orbits.Recently a linear feedback control law based upon the Lyapunov–Krasovskii (LK) method was developed for the suppression of chaotic oscillations.[15]A linear state feedback controller was designed to solve the chaos control problem of a class of new chaotic system in Ref.[16].(ⅱ) structure variation control [12-16]Since Y u X proposed structure variation method for controlling chaos of Lorenz system,[17]some improved sliding-mode control strategies were*Project supported by the National Natural Science Foundation of C hina (Grant No 50376029). †Corresponding au thor. E-mail:zibotll@introduced in chaos control. In Ref.[18] the author used a newly developed sliding mode controller with a time-varying manifold dynamic to compensate the external excitation in chaotic systems. In Ref.[19] the design schemes of integration fuzzy sliding-mode control were addressed, in which the reaching law was proposed by a set of linguistic rules. A radial basis function sliding mode controller was introduced in Ref.[20] for chaos control.(ⅲ) nonlinear geometric controlNonlinear geometric control theory was introduced for chaos control in Ref.[22], in which a Lorenz system model slightly different from the original Lorenz system was studied considering only the Prandtl number variation and process noise. In Ref.[23] the state space exact linearization method was also used to stabilize the equilibrium of the Lorenz system with a controllable Rayleigh number. (ⅳ)intelligence control[24-27 ]An intelligent control method based on RBF neural network was proposed for chaos control in Ref.[24]. Liu H, Liu D and Ren H P suggested in Ref.[25] to use Least-Square Support V ector Machines to drive the chaotic system to desirable points. A switching static output-feedback fuzzy-model-based controller was studied in Ref.[27], which was capable of handling chaos.Other methods are also attentively studied such as entrainment and migration control, impulsive control method, optimal control method, stochastic control method, robust control method, adaptive control method, backstepping design method and so on. A detailed survey of recent publications on control of chaos can be referenced in Refs.[28-34] and the references therein.Among most of the existing control strategies, it is considered essentially to know the model parameters for the derivation of a controller and the control goal is often to stabilize the embedded unstable period orbits of chaotic systems or to control the system to its equilibrium points. In case of controlling the system to its equilibrium point, one general approach is to linearize the system in the given equilibrium point, then design a controller with local stability, which limits the use of the control scheme. Based on Machine Learning methods, such as neural network method[24]or support vector machine method,[25]the control performance often depends largely on the training samples, and sometimes better generalization capability can not be guaranteed.Chaos, as the special phenomenon of deterministic nonlinear system, nonlinearity is the essence. So if a nonlinear real-time compensator can eliminate the effect of the system nonlinearities, chaotic motion is expected to be suppressed. Consequently the chaotic system can be controlled to a desired state. Under the guidance of nonlinear control theory, the objective of this paper is to design a control system to drive the chaotic systems not only to steady states but also to periodic trajectories. In the next section the controller architecture is introduced. In section 3, a Lorenz system considering parametric uncertainties and external disturbances is studied as an illustrative example. Two control schemes are designed for the studied chaotic system. By constructing appropriate L yapunov functions, after rigorous analysis from L yapunov stability theory sufficient conditions for the choice of control parameters are deduced for each scheme. Then in section 4 we present the numerical simulation results to illustrate the effectiveness of the design techniques. Finally some conclusions are provided to close the text.2. Controller architectureSystem differential equation is only an approximate description of the actual plant due to various uncertainties and disturbances. Without loss of generality let us consider a nonlinear continuous dynamic system, which appears strange attractors under certain parameter conditions. With the relative degree r n(n is the dimension of the system), it can be directly described or transformed to the following normal form:121(,,)((,,)1)(,,,)(,,)r r r z z z z za z v wb z v u u d z v u u vc z v θθθθθθθθ-=⎧⎪⎪⎪=⎪=+∆+⎨⎪ ++∆-+⎪⎪ =+∆+⎪=+∆⎩ (1) 1y z =where θ is the parameter vector, θ∆ denotes parameter uncertainty, and w stands for the external disturbance, such that w M ≤with Mbeingpositive.In Eq.(1)1(,,)T r z z z = can be called external state variable vector,1(,,)T r n v v v += called internal state variable vector. As we can see from Eq.(1)(,,,,)(,,)((,,)1)d z v w u a z v w b z v uθθθθθθ+∆=+∆+ ++∆- (2)includes system nonlinearities, uncertainties, external disturbances and so on.According to the chaotic system (1), the following assumptions are introduced in order to establish the results concerned to the controller design (see more details in Ref.[38]).Assumption 1 The relative degree r of the chaotic system is finite and known.Assumption 2 The output variable y and its time derivatives i y up to order 1r -are measurable. Assumption 3 The zero dynamics of the systemis asymptotically stable, i.e.,(0,,)v c v θθ=+∆ is asymptotically stable.Assumption 4 The sign of function(,,)b z v θθ+∆is known such that it is always positive or negative.Since maybe not all the state vector is measurable, also (,,)a z v θθ+∆and (,,)b z v θθ+∆are not known, a controller with integral action is introduced to compensate theinfluenceof (,,,,)d z v w u θθ+∆. Namely,01121ˆr r u h z h z h z d------ (3) where110121112100ˆr i i i r r r r i i ii r i i d k z k k k z kz k uξξξ-+=----++-==⎧=+⎪⎪⎨⎪=----⎪⎩∑∑∑ (4)ˆdis the estimation to (,,,,)d z v w u θθ+∆. The controller parameters include ,0,,1i h i r =- and ,0,,1i k i r =- . Here011[,,,]Tr H h h h -= is Hurwitz vector, such that alleigenvalues of the polynomial121210()rr r P s s h sh s h s h --=+++++ (5)have negative real parts. The suitable positive constants ,0,,1i h i r =- can be chosen according to the expected dynamic characteristic. In most cases they are determined according to different designed requirements.Define 1((,,))r k sign b z v θμ-=, here μstands for a suitable positive constant, and the other parameters ,0,,2i k i r =- can be selected arbitrarily. After011[,,,]Tr H h h h -= is decided, we can tune ,0,,1i k i r =- toachievesatisfyingstaticperformances.Remark 1 In this section, we consider a n-dimensional nonlinear continuous dynamic system with strange attractors. By proper coordinate transformation, it can be represented to a normal form. Then a control system with a nonlinear compensator can be designed easily. In particular, the control parameters can be divided into two parts, which correspond to the dynamic characteristic and the static performance respectively (The theoretic analysis and more details about the controller can be referenced to Ref.[38]).3. Illustrative example-the Lorenz systemThe Lorenz system captures many of the features of chaotic dynamics, and many control methods have been tested on it.[17,20,22-23,27,30,32-35,42] However most of the existing methods is model-based and has not considered the influence ofpersistent external disturbances.The uncontrolled original Lorenz system can be described by112121132231233()()()()x P P x P P x w x R R x x x x w xx x b b x w =-+∆++∆+⎧⎪=+∆--+⎨⎪=-+∆+⎩ (6) where P and R are related to the Prendtl number and Rayleigh number respectively, and b is a geometric factor. P ∆, R ∆and b ∆denote the parametric variations respectively. The state variables, 1x ,2x and 3x represent measures of fluid velocity and the spatial temperature distribution in the fluid layer under gravity , and ,1,2,3i w i =represent external disturbance. In Lorenz system the desired response state variable is 1x . It is desired that 1x is regulated to 1r x , where 1r x is a given constant. In this section we consider two control schemes for system (6).3.1 Control schemes for Lorenz chaotic system3.1.1 Control scheme 1The control is acting at the right-side of the firstequation (1x), thus the controlled Lorenz system without disturbance can be depicted as1122113231231x Px Px u xRx x x x x x x bx y x =-++⎧⎪=--⎨⎪=-⎩= (7) By simple computation we know system (7) has relative degree 1 (i.e., the lowest ordertime-derivative of the output y which is directly related to the control u is 1), and can be rewritten as1122113231231z Pz Pv u vRz z v v v z v bv y z =-++⎧⎪=--⎨⎪=-⎩= (8) According to section 2, the following control strategy is introduced:01ˆu h z d=-- (9) 0120010ˆ-d k z k k z k uξξξ⎧=+⎪⎨=--⎪⎩ (10) Theorem 1 Under Assumptions 1 toAssumptions 4 there exists a constant value *0μ>, such that if *μμ>, then the closed-loop system (8), (9) and (10) is asymptotically stable.Proof Define 12d Pz Pv =-+, Eq.(8) can be easily rewritten as1211323123z d u v Rz z v v vz v bv =+⎧⎪=--⎨⎪=-⎩ (11) Substituting Eq.(9) into Eq.(11) yields101211323123ˆz h z d dv R z z v v v z v bv ⎧=-+-⎪=--⎨⎪=-⎩ (12) Computing the time derivative of d and ˆdand considering Eq.(12) yields12011132ˆ()()dPz Pv P h z d d P Rz z v v =-+ =--+- +-- (13) 0120010000100ˆ-()()ˆ=()d k z k k z k u k d u k d k z k d d k dξξξ=+ =--++ =-- - = (14)Defining ˆdd d =- , we have 011320ˆ()()dd d P h P R z P z v P v P k d=- =+- --+ (15) Then, we can obtain the following closed-loop system101211323123011320()()z h z dvRz z v v v z v bv d Ph PR z Pz v Pv P k d⎧=-+⎪=--⎪⎨=-⎪⎪=+---+⎩ (16) To stabilize the closed-loop system (16), a L yapunovfunction is defined by21()2V ςς=(17)where, ςdenotes state vector ()123,,,Tz v v d, isthe Euclidean norm. i.e.,22221231()()2V z v v dς=+++ (18) We define the following compact domain, which is constituted by all the points internal to the superball with radius .(){}2222123123,,,2U z v v d zv v dM +++≤(19)By taking the time derivative of ()V ςand replacing the system expressions, we have11223322*********01213()()(1)V z z v v v v dd h z v bv k P d R z v P R P h z d P v d P z v d ς=+++ =----++ +++-- (20) For any ()123,,,z v v d U ∈, we have: 222201230120123()()(1)V h z v b v k P dR z v PR Ph z d P v d d ς≤----+ ++++ ++ (21)Namely,12300()(1)22020V z v v dPR Ph R h R P ς⎡⎤≤- ⎣⎦++ - 0 - - 1 - 2⨯00123(1)()2Tb PR Ph P k P z v v d ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥0 ⎢⎥2⎢⎥++⎢⎥- - - +⎢⎥⎣22⎦⎡⎤⨯ ⎣⎦(22) So if the above symmetrical parameter matrix in Eq.(22) is positive definite, then V is negative and definite, which implies that system (16) is asymptotically stable based on L yapunov stability theory.By defining the principal minor determinants of symmetrical matrix in Eq.(22) as ,1,2,3,4i D i =, from the well-known Sylvester theorem it is straightforward to get the following inequations:100D h => (23)22004RD h =-> (24)23004R b D bh =-> (25)240302001()(1)(2)821[2(1)]08P M D k P D b PR Ph PR D Pb Ph R PR Ph =+-+++--+++>(26)After 0h is determined by solving Inequalities (23) to (25), undoubtedly, the Inequalities (26) can serve effectively as the constraints for the choice of 0k , i.e.20200031(1)(2)821[2(1)]8P M b PR Ph PR D Pb Ph R PR Ph k P D ++++ ++++>- (27)Here,20200*31(1)(2)821[2(1)]8P M b PR Ph PR D Pb Ph R PR Ph P D μ++++ ++++=-.Then the proof of the theorem 1 is completed. 3.1.2 Control scheme 2Adding the control signal on the secondequation (2x ), the system under control can be derived as112211323123x P x P x x R x x x x u xx x bx =-+⎧⎪=--+⎨⎪=-⎩ (28) From Eq.(28), for a target constant 11()r x t x =,then 1()0xt = , by solving the above differential equation, we get 21r r x x =. Moreover whent →∞,3r x converges to 12r x b . Since 1x and 2x havethe same equilibrium, then the measured state can also be chosen as 2x .To determine u , consider the coordinate transform:122133z x v x v x=⎧⎪=⎨⎪=⎩ and reformulate Eq.(28) into the following normal form:1223121231231zRv v v z u vPz Pv v z v bv y z =--+⎧⎪=-⎨⎪=-⎩= (29) thus the controller can be derived, which has the same expression as scheme 1.Theorem 2 Under Assumptions 1, 2, 3 and 4, there exists a constant value *0μ>, such that if *μμ>, then the closed-loop system (9), (10) and (29) is asymptotically stable.Proof In order to get compact analysis, Eq.(29) can be rewritten as12123123z d u v P z P v vz v bv =+⎧⎪=-⎨⎪=-⎩ (30) where 2231d Rv v v z =--Substituting Eq.(9) into Eq.(30),we obtain:1012123123ˆz h z d dv P z P v v z v bv ⎧=-+-⎪=-⎨⎪=-⎩ (31) Giving the following definition:ˆdd d =- (32) then we can get22323112123212301()()()()dRv v v v v z R Pz Pv Pz Pv v v z v bv h z d =--- =--- ----+ (33) 012001000ˆ-()d k z k k z k u k d u k dξξ=+ =--++ = (34) 121232123010ˆ()()()(1)dd d R Pz Pv Pz Pv v v z v bv h z k d=- =--- --+-+ (35)Thus the closed-loop system can be represented as the following compact form:1012123123121232123010()()()(1)zh z d v Pz Pv v z v bv d R Pz Pv Pz Pv v v z v bv h z k d⎧=-+⎪⎪=-⎪=-⎨⎪=---⎪⎪ --+-+⎩(36) The following quadratic L yapunov function is chosen:21()2V ςς=(37)where, ςdenotes state vector ()123,,,Tz v v d , is the Euclidean norm. i.e.,22221231()()2V z v v dς=+++ (38) We can also define the following compact domain, which is constituted by all the points internalto the super ball with radius .(){}2222123123,,,2U z v v d zv v dM =+++≤ (39)Differentiating V with respect to t and using Eq.(36) yields112233222201230011212322321312()(1)(1)()V z z v v v v dd h z P v bv k dP R h z d P z v z v v P b v v d P v d P z v d z v d ς=+++ =----+ +++++ ++--- (40)Similarly, for any ()123,,,z v v d U ∈, we have: 2222012300112133231()(1)(1)(2V h z P v b v k dPR h z d P z v v P b d P v d d M z dς≤----+ +++++ ++++ + (41)i.e.,12300()(12)22V z v v dPR M h P h P Pς⎡⎤≤- ⎣⎦+++ - -2 - 0 ⨯ 001230(12)(1)2TP b PR M h P k z v v d ⎡⎤⎢⎥⎢⎥⎢⎥ - ⎢⎥⎢⎥⎢⎥ ⎢⎥22⎢⎥⎢⎥ +++ - - -+⎢⎥⎣22⎦⎡⎤⨯ ⎣⎦(42) For brevity, Let1001(12)[(222)82(23)]P PR M h b PR P h M P b α=++++++ ++(43) 2201[(231)(13)]8P M P b b PR h α=+-+++ (44)230201(2)[2(12)8(2)(4)]PM P b P P PR M h P b Ph P α=++ +++ ++- (45)Based on Sylvester theorem the following inequations are obtained:100D h => (46)22004PD h P =-> (47)3202PMD bD =-> (48)403123(1)0D k D ααα=+---> (49)where,1,2,3,4i D i =are the principal minordeterminants of the symmetrical matrix in Eq.(42).*0k μ>*12331D αααμ++=- (50)The theorem 2 is then proved.Remark 2 In this section we give two control schemes for controlling chaos in Lorenz system. For each scheme the control depends on the observed variable only, and two control parameters are neededto be tuned, viz. 0h and 0k . According to L yapunov stability theory, after 0h is fixed, the sufficient condition for the choice of parameter 0k is also obtained.4. Simulation resultsChoosing 10P =,28R =, and 8/3b =, the uncontrolled Lorenz system exhibits chaotic behavior, as plotted in Fig.1. In simulation let the initial values of the state of thesystembe 123(0)10,(0)10,(0)10x x x ===.x1x 2x1x 3Fig.1. C haotic trajectories of Lorenz system (a) projected on12x x -plane, (b) projected on 13x x -plane4.1 Simulation results of control the trajectory to steady stateIn this section only the simulation results of control scheme 2 are depicted. The simulation results of control scheme 1 will be given in Appendix. For the first five seconds the control input is not active, at5t s =, control signal is input and the systemtrajectory is steered to set point2121(,,)(8.5,8.5,27.1)T Tr r r x x x b =, as can be seen inFig.2(a). The time history of the L yapunov function is illustrated in Fig.2(b).t/sx 1,x 2,x 3t/sL y a p u n o v f u n c t i o n VFig.2. (a) State responses under control, (b) Time history of the Lyapunov functionA. Simulation results in the presence ofparameters ’ changeAt 9t s =, system parameters are abruptly changed to 15P =,35R =, and 12/3b =. Accordingly the new equilibrium is changedto 2121(,,)(8.5,8.5,18.1)T Tr r r x x x b =. Obviously, aftervery short transient duration, system state converges to the new point, as shown in Fig.3(a). Fig.4(a) represents the evolution in time of the L yapunov function.B. Simulation results in the presence of set pointchangeAt 9t s =, the target is abruptly changedto 2121(,,)(12,12,54)T Tr r r x x x b =, then the responsesof the system state are shown in Fig.3(b). In Fig.4(b) the time history of the L yapunov function is expressed.t/sx 1,x 2,x 3t/sx 1,x 2,x 3Fig.3. State responses (a) in the presence of parameter variations, (b) in the presence of set point changet/sL y a p u n o v f u n c t i o n Vt/sL y a p u n o v f u n c t i o n VFig.4. Time history of the Lyapunov fu nction (a) in the presence of parameter variations, (b) in the presence of set point changeC. Simulation results in the presence ofdisturbanceIn Eq.(5)external periodic disturbance3cos(5),1,2,3i w t i π==is considered. The time responses of the system states are given in Fig.5. After control the steady-state phase plane trajectory describes a limit cycle, as shown in Fig.6.t/sx 1,x 2,x 3Fig.5. State responses in the presence of periodic disturbancex1x 3Fig.6. The state space trajectory at [10,12]t ∈in the presence ofperiodic disturbanceD. Simulation results in the presence of randomnoiseUnder the influence of random noise,112121132231233xPx Px x Rx x x x u xx x bx εδεδεδ=-++⎧⎪=--++⎨⎪=-+⎩ (51) where ,1,2,3i i δ= are normally distributed withmean value 0 and variance 0.5, and 5ε=. The results of the numerical simulation are depicted in Fig.7,which show that the steady responses are hardly affected by the perturbations.t/sx 1,x 2,x 3t/se 1,e 2,e 3Fig.7. Time responses in the presence of random noise (a) state responses, (b) state tracking error responses4.2 Simulation results of control the trajectory to periodic orbitIf the reference signal is periodic, then the system output will also track this signal. Figs.8(a) to (d) show time responses of 1()x t and the tracking trajectories for 3-Period and 4-period respectively.t/sx 1x1x 2t/sx 1x1x 2Fig.8. State responses and the tracking periodic orbits (a)&( b)3-period, (c)&(d) 4-periodRemark 3 The two controllers designed above solved the chaos control problems of Lorenz chaoticsystem, and the controller design method can also beextended to solve the chaos suppression problems of the whole Lorenz system family, namely the unified chaotic system.[44-46] The detail design process and close-loop system analysis can reference to the author ’s another paper.[47] In Ref.[47] according to different positions the scalar control input added,three controllers are designed to reject the chaotic behaviors of the unified chaotic system. Taking the first state 1x as the system output, by transforming system equation into the normal form firstly, the relative degree r (3r ≤) of the controlled systems i s known. Then we can design the controller with the expression as Eq.(3) and Eq.(4). Three effective adaptive nonlinear compensating mechanisms are derived to compensate the chaotic system nonlinearities and external disturbances. According toL yapunov stability theory sufficient conditions for the choice of control parameters are deduced so that designers can tune the design parameters in an explicit way to obtain the required closed loop behavior. By numeric simulation, it has been shown that the designed three controllers can successfully regulate the chaotic motion of the whole family of the system to a given point or make the output state to track a given bounded signal with great robustness.5. ConclusionsIn this letter we introduce a promising tool to design control system for chaotic system subject to persistent disturbances, whose entire dynamics is assumed unknown and the state variables are not completely measurable. By integral action the nonlinearities, including system structure nonlinearity, various disturbances, are compensated successfully. It can handle, therefore, a large class of chaotic systems, which satisfy four assumptions. Taking chaotic Lorenz system as an example, it has been shown that the designed control scheme is robust in the sense that the unmeasured states, parameter uncertainties and external disturbance effects are all compensated and chaos suppression is achieved. Some advantages of this control strategy can be summarized as follows: (1) It is not limited to stabilizing the embeddedperiodic orbits and can be any desired set points and multiperiodic orbits even when the desired trajectories are not located on the embedded orbits of the chaotic system.(2) The existence of parameter uncertainty andexternal disturbance are allowed. The controller can be designed according to the nominal system.(3) The dynamic characteristics of the controlledsystems are approximately linear and the transient responses can be regulated by the designer through controllerparameters ,0,,1i h i r =- .(4) From L yapunov stability theory sufficientconditions for the choice of control parameters can be derived easily.(5) The error converging speed is very fast evenwhen the initial state is far from the target one without waiting for the actual state to reach the neighborhood of the target state.AppendixSimulation results of control scheme 1.t/sx 1,x 2,x 3t/sL y a p u n o v f u n c t i o n VFig.A1. (a) State responses u nder control, (b) Time history of the Lyapunov functiont/sx 1,x 2,x 3t/sx 1,x 2,x 3Fig.A2. State responses (a) in the presence of parameter variations, (b) in the presence of set point changet/sL y a p u n o v f u n c t i o n Vt/sL y a p u n o v f u n c t i o n VFig.A3. Time history of the L yapu nov fu nction (a) in the presence of parameter variations, (b) in the presence of set point changet/sx 1,x 2,x 3Fig.A4. State responses in the presence of periodic disturbanceresponsest/sx 1,x 2,x 3Fig.A5. State responses in the presence of rand om noiset/sx 1x1x 2Fig.A6. State response and the tracking periodic orbits (4-period)References[1] Lü J H, Zhou T S, Zhang S C 2002 C haos Solitons Fractals 14 529[2] Yoshihiko Nagai, Hua X D, Lai Y C 2002 C haos Solitons Fractals 14 643[3] Li R H, Xu W , Li S 2007 C hin.phys.16 1591 [4]Xiao Y Z, Xu W 2007 C hin.phys.16 1597[5] Ott E ,Greb ogi C and Yorke J A 1990 Phys.Rev .Lett. 64 1196 [6]Yoshihiko Nagai, Hua X D, Lai Y C 1996 Phys.Rev.E 54 1190 [7] K.Pyragas, 1992 Phys. Lett. A 170 421 [8] Lima,R and Pettini,M 1990 Phys.Rev.A 41 726[9] Zhou Y F, Tse C K, Qiu S S and Chen J N 2005 C hin.phys. 14 0061[10] G .Cicog na, L.Fronzoni 1993 Phys.Rew .E 30 709 [11] Rakasekar,S. 1993 Pramana-J.Phys.41 295 [12] Gong L H 2005 Acta Phys.Sin.54 3502 (in C hinese) [13] Chen L,Wang D S 2007 Acta Phys.Sin.56 0091 (in C hinese) [14] C hen G R and Dong X N 1993 IEEE Trans.on Circuits andSystem-Ⅰ:Fundamental Theory and Applications 40 9 [15] J.L. Kuang, P.A. Meehan, A.Y.T. Leung 2006 C haos SolitonsFractals 27 1408[16] Li R H, Xu W, Li S 2006 Acta Phys.Sin.55 0598 (in C hinese) [17] Yu X 1996 Int.J.of Systems Science 27 355[18] Hsun-Heng Tsai, C hyu n-C hau Fuh and Chiang-Nan Chang2002 C haos,Solitons Fractals 14 627[19] Her-Terng Yau and C hieh-Li C hen 2006 C hao ,SolitonsFractal 30 709[20] Guo H J, Liu J H, 2004 Acta Phys.Sin.53 4080 (in C hinese) [21] Yu D C, Wu A G , Yang C P 2005 Chin.phys.14 0914 [22] C hyu n-C hau Fuh and Pi-Cheng Tu ng 1995 Phys.Rev .Lett.752952[23] Chen L Q, Liu Y Z 1998 Applied Math.Mech. 19 63[24] Liu D, R en H P, Kong Z Q 2003 Acta Phys.Sin.52 0531 (inChinese)[25] Liu H, Liu D and Ren H P 2005 Acta Phys.Sin.54 4019 (inChinese)[26] C hang W , Park JB, Joo YH, C hen GR 2002 Inform Sci 151227[27] Gao X, Liu X W 2007 Acta Phys.Sin. 56 0084 (in C hinese) [28] Chen S H, Liu J, Lu J 2002 C hin.phys.10 233 [29] Lu J H, Zhang S. 2001 Phys. Lett. A 286 145[30] Liu J, Chen S H, Xie J. 2003 C haos Solitons Fractals 15 643 [31] Wang J, Wang J, Li H Y 2005 C haos Solitons Fractals 251057[32] Wu X Q, Lu JA, C hi K. Tse, Wang J J, Liu J 2007 ChaoSolitons Fractals 31 631[33] A.L.Fradkov , R .J.Evans, 2002 Preprints of 15th IF AC W orldCongress on Automatic Control 143[34] Zhang H G 2003 C ontrol theory of chaotic systems (Shenyang:Northeastern University) P38 (in C hinese)[35] Yu-Chu Tian, Moses O.Tadé, David Levy 2002Phys.Lett.A.296 87[36] Jose A R , Gilberto E P, Hector P, 2003 Phys. Lett. A 316 196 [37] Liao X X, Yu P 2006 Chaos Solitons Fractals 29 91[38] Tornambe A, V aligi P.A 1994 Measurement, and C ontrol 116293[39] Andrew Y.T.Leung, Liu Z R 2004 Int.J.Bifurc.C haos 14 2955 [40] Qu Z L, Hu,G .,Yang,G J, Qin,G R 1995 Phys.Rev .Lett.74 1736 [41] Y ang J Z, Qu Z L, Hu G 1996 Phys.Rev.E.53 4402[42] Shyi-Kae Yang, C hieh-Li Chen, Her-Terng Yau 2002 C haosSolitons Fractals 13 767。
美赛外文数据库

美赛外文数据库评价模型:层次分析模糊数学灰色系统投资组合的熵理论和信息价值literaturetomotivateyourmodelandtovalidateyourreult.egregioulyimplifietheproblem,itmayprovideinighttothebaicbeha viorandinpireimprovedmodel.Thiprogreionleadtoarepectableandatifyingfinalmodel,wit hitvaliditytemmingfromtheprecedingmodel.建立多个模型。
即使一开始的模型假设太过过分,但可以使我们对原理有所了解。
因此,如果有的模型错了的话,可以写上,但要注明错在哪里the\A连续B离散问题三大组织:TheIntituteforOperationReearchandtheManagementScience(INFORM S)(管理科学与运用协会)TheSocietyforIndutrialandAppliedMathematic(SIAM)(美国工程与运用数学协会)TheMathematicalAociationofAmerica(MAA)(美国数学协会)国家科技图书文献中心中科院武汉情报中心:外文数据库”1同济图书馆:美国生态学会ESA运筹学和管理学研究协会是出题的组织之一发布最新的运筹学和管理学方法和应用。
此外,INFORMS还不定期组织众多的专业会议。
Paword:Inform05具体11个杂志的研究领域:1.DeciionAnalyi---决策分析本期刊创刊于2004年,是同行评审性季刊。
推动决策分析的理论、应用和教学。
本刊主要聚焦于发展和研究运筹决策制定方法,吸收决策理论和决策分析的各个方面,为决策制定者们提供实用的指导。
rmationSytemReearch信息系统研究国际领先的理论、研究和学术思想发展的期刊,聚焦于组织、机构、社会的信息系统,为人类组织和管理提供更多的信息技术生产应用知识关于运筹研究和计算交集的论文。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
errors, conflicting data and site heterogeneities all contribute to the complexities in site characterization. Detailed site information needed to fully characterize a site are usually lacking as large number of samples are required to completely capture a site's descriptions and its variabilities. Substantial data requirement means high costs and usually the site assessor makes decisions with much less data. Apart from issues related to site variabilities, multiple sources of information have to be compiled such as historical, geological and hydrogeologic information which may be in “soft” (descriptive) or “hard” (numeric) form. There may be data from laboratory studies, field investigations as well as from expert opinions. All these information are required to determine soil characteristics and site parameters such as hydraulic conductivity, storage coefficient, and porosity which
Kejiang Zhang a, Gopal Achari a,⁎, Hua Li b
a b
Department of Civil Engineering, University of Calgary, Calgary, AB, Canada T2N 1N4 Department of Mathematics, Zhengzhou University, Zhengzhou, Henan, PR China
Journal of Contaminant Hydrology 110 (2009) 45–59
Contents lists available at ScienceDirect
Journal of Contaminant Hydrology
j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j c o n h yd
A comparison of numerical solutions of partial differential equations with probabilistic and possibilistic parameters for the quantification of uncertainty in subsurface solute transport
1. Introduction Modeling contaminant fate and transport is an integral part of exposure assessment, a necessary step in environmental risk assessment. The various physical phenomena in groundwater flow and contaminant transport are represented by partial differential equations (PDEs) in space and time. The geologic and hydrogeologic site characteristics are included as parameters in the governing PDEs. To determine these parameters is a challenging part of modeling. Scarcity of information, uncertainty in the data collected, mea. / Journal of Contaminant Hydrology 110 (2009) 45–59
are then integrated into a groundwater flow and contaminant transport model. For contaminant fate and transport, additionally, the biogeochemical information about the target contaminant such as precipitation–dissolution, redox reaction, and biodegradation are needed. Uncertain, variable and multi-sourced “soft” and “hard” information not only require proper representations of different uncertainties but also their integration (Porter et al., 2000; Dubois et al., 2000; Sentz and Ferson, 2002; Dubois et al., 2004), and incorporation of the integrated information into the groundwater flow and contaminant transport modeling. Three types of uncertainty (Möller and Beer, 2004), stochastic, informal and lexical, are present in site characterization. Stochastic uncertainty is best described by classical probability theory. Informal uncertainty results from an information deficit, such as when only a small number of observations are available. Lexical uncertainty occurs when instead of numbers, words such as “high”, “medium”, and “low” are used. Opinions of experts are usually linguistic and will have this type of uncertainty. Probability theory is suitable for representation of probabilistic information. Fuzzy set theory (Zadeh, 1965), and possibility theory (Zadeh, 1978; Dubois and Prade, 1988), evidence theory (Shafer, 1976), and random sets (Matheron, 1975; Molchanov, 2005; Nguyen, 2006) are employed to represent informal and lexical uncertain information or possibilistic information (Liu and Peng, 2005). Evidence theory is more suitable when a combination of different types of uncertain information exist (Ross, 2004). This may occur when one theory best represents one parameter whereas another theory may be more suitable for another parameter. Both types of parameters are however required in one governing equation. This leads to the interesting question of hybrid uncertainties and their propagation. Recent research by Guyonnet et al. (2003) and Baudrit et al. (2004, 2007) identified a method of hybrid uncertainty propagation through the conventional risk assessment process. Chen et al. (2003), Li et al. (2007) and Li et al. (2008) use the hybrid fuzzy-stochastic modeling approach to evaluate risk posed by contaminated groundwater and air pollution. However, different types of uncertainties associated with contaminant fate and transport were all captured using probability theory and fuzzy sets were only used to quantify uncertainties inherent in environmental quality guidelines and health evaluation criteria. Baraldi and Zio (2008) used a combined Monte Carlo and possibilistic approach to propagate uncertainty in event tree analysis. Conventionally, uncertainties associated with different parameters in a PDE are modeled using stochastic partial differential equations (SPDEs) (Dagan and Neuman, 1997; Glimm and Sharp, 1998). There are two limitations of stochastic modeling: (1) it needs sufficient data to justify the assumed statistical distributions. Given sufficient statistical data, the probability theory can better capture the random uncertainties; (2) different types of uncertainty that may be associated with the parameters cannot all be properly represented by probability distributions. Chen et al. (2003), Li et al. (2007), and Maxwell et al. (1998) used stochastic methods to simulate the uncertainties in groundwater flow and contaminant transport. Both used the Monte Carlo method for obtaining the simulation results.