数值积分法.doc

第7章 数值积分法7.1 实验目的 了解求积公式及代数精度概念，理解并掌握求定积分的求积公式的算法构造和计算，学习用计算机求定积分的一些科学计算方法和简单的编程技术和能用程序实现这些算法。7.2 概念与结论1. 求积公式计算定积分的如下形式的近似公式：称为求积公式。2代数精度 若求积公式 对一切不高于m次的 多项都准确成立，而对于m+1次多项式等号不成立，则称此求积公式 的代数精度为m。 代数精度越高，求积公式越好。3.求积余项4.Newton-Cotes求积公式的代数精度 n点Newton-Cotes求积公式的代数精度至少可以达到n-1,且当n为奇数时，可以达到n。5.Richardson外推定理 设函数F1(h)逼近量F*的余项为：F*-F1(h)=a1h p1+a2h p2+a k p k+式中p kp k-1p2p10, F*和a i (i=1,2, )都是与h无关的常数,且k1时，a k 0，则由：定义的函数F2(h)也逼近F*，且有F*-F2(h)= b2h p2+b k p k+6. 关于复合梯形公式的展开定理设f(x)在a,b区间上无穷次可微，则有如下展开式：T(h)=I+a1h2+a2h4+a3h6+amh2m+式中T(h)是函数f(x)在a,b区间上的复化梯形值Tn,7.3 程序中Mathematica语句解释 1. 随机函数 Random 随机给出闭区间0,1内的一个实数 RandomReal, xmax 随机给出闭区间0,xmax内的一个实数 RandomReal, xmin, xmax 随机给出闭区间xmin,xmax内的一个实数 RandomInteger 随机给出整数0或1 RandomInteger, xmin, xmax 随机给出xmin到xmax之间的一个整数 RandomComplex 随机给出单位正方形内的一个复数2a1,a2,an 表示由元素a1,a2,an组成的一个表，元素可以是任何内容。如：1，3，4，5，1,x,2,3,x+y，1，3，1，2，3，3，2，4等3listk 表list中的第k个元素4listi,j 表list中第i个元素中的第j个元素,此时list中的第i个元素应该也是一个表。 7.4 方法、程序、实验 在实际问题中，往往会遇到被积函数f(x)的原函数无法用初等函数来表示，或函数只能用表格表示，或有的虽然能用初等函数表示，但过分复杂，所以这些情形都需要去建立定积分的近似计算公式来做积分计算。数值积分是进行定积分计算的一种方法，它可以解决不能用定积分基本公式计算的所有定积分问题。数值积分涉及很多计算公式，这里主要介绍Newton-Cotes求积公式、复合求积公式、Romberg求积方法和Monte-Carlo方法的构造过程和算法程序。1. n点 NewtonCotes求积公式 n点 NewtonCotes求积公式又称为等距节点求积公式，它是利用被积函数f(x)在积分区间a,b的n个等分节点上的函数值构造的插值函数j(x)代替f(x)做定积分计算所构造求积公式。这个求积公式是通常做定积分近似计算的梯形公式和抛物线公式的推广，主要在理论上用的多些。1) n点 NewtonCotes求积公式的构造过程 将积分区间a,b 分为n-1等分，其中n个节点 xi=a+(i-1h, i=1,2,n，h=(b-a)/(n-1)，然后用f(x)在这n个节点上建立插值于f(x)的n-1次代数多项式Pn-1(x)，引入变换x=a+th, 0tn-1则有 带入定积分，有：Ck(n)称为Cotes(柯特斯)系数, 则得到n点 NewtonCotes求积公式:n点 NewtonCotes求积公式的求积余项为当n=2时，2点的 NewtonCotes求积公式就是如下梯形公式： 梯形求积公式求积余项为当n=3时，3点的NewtonCotes求积公式就是如下抛物线（Simpson）公式： Simpson求积公式求积余项为如果想得到其他的 NewtonCotes求积公式只要在有关书中查出Cotes系数表就可以马上得到相应的NewtonCotes求积公式。2) n点 NewtonCotes求积公式算法 1 输入被积函数f(x)及积分上下限a,b2 选择Cotes系数构造求积公式3 用求积公式求定积分3) n点 NewtonCotes求积公式程序Cleara,b,x,n,s;a=Inputa=b=Inputb=fx_=Input被积函数f(x)=n= Input求积节点个数n=;c=1/2,1/2,1/6,4/6,1/6,1/8,3/8,3/8,1/8,7/90,16/45,2/15,16/45,7/90,19/288,25/96,25/144,25/144,25/96,19/288,41/840,9/35,9/280,34/105,9/280,9/35,41/840,751/17280,3577/17280,1323/17280,2989/17280,1323/17280,3577/17280,751/17280,989/28350,5888/28350,-928/28350,10496/28350,-4540/28350,10496/28350,-928/28350,5888/28350,989/28350;h=(b-a)/(n-1);x=Tablea+k*h,k,0,n-1;s=(b-a)*Sumcn-1,k*fxk,k,1,nPrint定积分=,Ns,8;说明 本程序用n（n=2,3,4,5,6，7，8，9）点 NewtonCotes求积公式求a,b上的定积分近似值。程序执行后，按要求通过键盘输入积分下限a、积分上限b、被积函数f(x)和求积节点个数n后，计算机则给出定积分的近似值。程序中变量说明a：存放积分下限b: 存放积分上限fx: 存放被积函数f(x)n: 存放求积节点个数c: 存放Cotes系数s: 存放定积分近似值h: 存放节点步长x：存放节点xi注 语句c=1/2,1/2,1/6,4/6,1/6,1/8,3/8,3/8,1/8,7/90.16/45,2/15,16/45,7/90, 19/288,25/96,25/144,25/144,25/96,19/288的第i个分量表是具有i个节点的Cotes系数。4) 例题与实验例1. 用n=3和n=4的Newton-Cotes求积公式求定积分的近似值。解：执行n点 NewtonCotes求积公式程后，在输入的窗口中按提示分别输入1、3、Exp-x/2、3每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出用n=3的Newton-Cotes求积公式计算出的定积分结果：定积分=0.76705953 再执行n点 NewtonCotes求积公式程后，在输入的窗口中按提示分别输入1、3、Exp-x/2、4每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出用n=4的Newton-Cotes求积公式计算出的定积分结果： 定积分=0.76691628因此用n=3和n=4的Newton-Cotes求积公式求本题定积分近似值分别为0.76705953和0.76691628注意到本题的精确值为0.766800999.,可见n=4的Newton-Cotes求积公式计算结果较好。2. 复化求积公式 复化求积公式是把积分区间分成若干小区间，在每个小区间上采用次数不高的插值公式，如梯形公式或抛物线公式，构造出相应的求积公式，然后再把它们加起来得到整个区间上的求积公式。复化求积公式克服了高次Newton-Cotes公式计算不稳定的问题，其运算简单且易于在计算机上实现。常用的复化求积公式是复化梯形公式和复化抛物线公式，下面分别讨论。1) 复化梯形公式的构造过程 把区间a,b n等分，取节点xk=a+kh,k=0,1,.n, h=(b-a)/n,对每个小区间xk,xk+1用梯形求积公式，再累加起来得：公式就是复化梯形公式。它的求积余项为： 由复化梯形公式的余项，可以看到，当n不断变大时， Tn无限接近定积分值。因此复化梯形公式可以使定积分的计算达到任意精度。为了计算简单，提高效率，常用|T2n Tn|e来得到满足精度要求的定积分值。2) 复化抛物线公式的构造过程在每个小区间xk,xk+1上，有 把区间a,b n等分，取节点xk=a+kh,k=0,1,.n, h=(b-a)/n,对每个小区间xk,xk+1 用抛物线求积公式，再累加起来得：公式就是复化抛物线公式。它的求积余项为： 由复化抛物线公式的余项，可以看到，当n不断变大时，Sn无限接近定积分值。因此复化抛物线公式也可以使定积分的计算达到任意精度，且收敛的速度比复化梯形公式更快。为了计算简单，提高效率，常用| S2n Sn|e来得到满足精度要求的定积分值。3) 复化梯形求积公式算法1.输入被积函数f(x),积分上下限a,b,和求积精度e2. n1, 计算Tn3. 计算T2n4 判断|T2n Tn|e是否成立，如果成立，输出定积分近似值，停止5 否则， Tn T2n , n2n,转34) 复化抛物线求积公式算法1.输入被积函数f(x),积分上下限a,b,和求积精度e2. n1, 计算Sn3.计算S2n4.判断|S2n Sn|eps,Printn=,2n,定积分值为,Nt2,10;Print误差=,er;h=h/2;t1=t2;n=n+1;t2=t1/2+h*Sumfa+k*h,k,1,2n,2;er=t2-t1/N;Printn=,2n,定积分值为,Nt2,10;Print误差=,er说明 本程序从梯形公式T1开始，用复合梯形求积公式求a,b上满足精度小于eps要求的定积分近似值。程序执行后，按要求通过键盘输入积分下限a、积分上限b、被积函数f(x)和精度要求eps 后，计算机则给出满足精度要求的定积分近似值及中间计算值和误差。程序中变量说明a：存放积分下限b: 存放积分上限fx: 存放被积函数f(x)eps: 存放求积精度eh: 存放节点步长x：为函数fx:提供变量t1: 存放复合梯形值Tnt2: 存放复合梯形值T2n和定积分近似值n: 存放复合梯形公式的节点加密次数er:存放误差注：为提高计算效率，计算T2n时使用了公式6) 复化抛物线求积公式程序(*复合抛物线公式*）Cleara,b,x,n,f,s;a=Inputa=b=Inputb=fx_=Input被积函数f(x)=n=Input分割区间数n=;h=(b-a)/n;a1=a+h/2;s=h/6*(fa+fb+ 2Sumfa+k*h,k,1,n-1+4Sumfa1+k*h,k,0,n-1);Printn=,n;Print定积分值为,Ns,10说明 本程序用复合抛物线求积公式求a,b上定积分近似值。程序执行后，按要求通过键盘输入积分下限a、积分上限b、被积函数f(x)和分割区间数n 后，计算机则给出定积分的近似值Sn。程序中变量说明a：存放积分下限b: 存放积分上限fx: 存放被积函数f(x)h: 存放节点步长x：为函数fx提供变量n: 存放分割区间数s: 存放复合抛物线值Sn和定积分近似值注：Mathematica数学软件中也有一个求数值积分的命令，命令形式为：NIntegratefx,x,a,b 它可以求出a,b上的定积分近似值. 7) 例题与实验例2. 用复合梯形公式求定积分的近似值,要求误差小于10-7。解：执行复化梯形求积公式程后，在输入的窗口中按提示分别输入：0、3、Exp-x2、10(-7)每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出的计算结果如下：n=2 定积分值为0.7313702518 误差=0.0474305n=4 定积分值为0.7429840978 误差=0.0116138n=8 定积分值为0.7458656148 误差=0.00288152n=16 定积分值为0.7465845968 误差=0.000718982n=32 定积分值为0.7467642547 误差=0.000179658n=64 定积分值为0.7468091636 误差=0.000044909n=128 定积分值为0.7468203905 误差=0.0000112269n=256 定积分值为0.7468231972 误差=2.806710-6n=512 定积分值为0.7468238989 误差=7.0167610-7n=1024 定积分值为0.7468240743 误差=1.7541910-7n=2048 定积分值为0.7468241182 误差=4.3854810-8从计算结果可以得出经过11次加密分割积分区间得到满足精度要求的定积分近似值0.7468241182。此外，通过计算过程可以明显看到复合梯形求积公式的收敛性。例3. 分别用n=1,2,4,8,16,32的复合抛物线求积公式求定积分的近似值。解：执行复化抛物线求积公式程后，在输入的窗口中按提示分别输入：0、 Pi、Expx*Cosx、1每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出的计算结果如下：n=1定积分值为-11.59283955再重复如上操作，将分割区间数分别输入2，4，8，16，32，得如下输出结果：n=2定积分值为-11.98494402n=4定积分值为-12.06420896n=8定积分值为-12.06995132n=16定积分值为-12.07032146n=32定积分值为-12.07034476。注意到本题的定积分值为-12.07034631639，从中可以看到复合抛物线求积公式收敛还是很快的。你能通过修改复合抛物线求积公式一次得出如上所有计算结果吗？3.Romberg求积公式 Romberg又称逐次分半加速收敛法，它是对复合梯形求积公式采用Richardson加速技术得到的一种数值求积方法。1) Romberg求积公式的构造过程 由复合梯形公式的展开定理，有如下关系式：T1(h)-I=a1h2+a2h4+a3h6+amh2m+这里利用Richardson外推定理对T1(h)进行加速，注意到这里p1=2,取q=1/2有于是得到收敛于I的加速公式T2(h),如果再T2(h) 进行加速，可以有继续加速下去，可以得到加速序列式中m代表对I进行的加速次数。 此外，注意到T1(h/2)=T2n ，且直接计算可以知道T2(h)就是复合抛物线公式Sn,于是有计算关系式类似地，还有如下计算格式：和 于是我们至少可以得到如下3个收敛于定积分I的数列：T型值数列：T1、T2、T4、T8 、T16 、T2n 、Simpson数列：S1、S2、S4、S8 、S16 、S2n 、Cotes数列：C1、C2、C4、C8 、C16 、C2n 、Romberg数列：R1、R2、R4、R8 、R16 、R2n 、利用Romberg数列求定积分I的方法称为Romberg求积方法。 如果给定求积精度e，可以用作为终止计算的条件，并取最后算出的R值作为积分值。2) Romberg求积算法1.输入被积函数f(x),积分上下限a,b,和求积精度e2.按顺序计算T1、T2、T4、T8 并做 t1 T1、t2 T2、t3 T4、t4 T8、 s1 (4t2-t1)/3,s2 (4t3-t2)/3,s3 (4t4-t3)/3;c1 (16s2-s1)/15,c2 (16s3-s2)/15,R2 (64c2-c1)/63;3.n1,R1 c2 , t1t4,s1s3 ,c1c24.判断|R2 R1|eps, r1=r2; Printr(,k,)=,Nr2,20, eps=,Ner,10; h=h/2;m=2m; t2=t1/2+h*Sumfa+i*h,i,1,m,2/N; s2=(4t2-t1)/3; c2=(16s2-s1)/15; r2=(64c2-c1)/63; t1=t2; s1=s2; er=r2-r1; k=k+1; c1=c2;Printr(,k,)=,Nr2,20, eps=,Ner,10;说明 本程序从用Romberg求积方法求a,b上满足精度小于eps要求的定积分近似值。程序执行后，按要求通过键盘输入积分下限a、积分上限b、被积函数f(x)和精度要求eps 后，计算机则给出满足精度要求的定积分近似值及中间计算值和误差。程序中变量说明a：存放积分下限b: 存放积分上限fx: 存放被积函数f(x)eps: 存放求积精度eh: 存放节点步长x：为函数fx:提供变量t,t1 ,t2: 存放复合梯形值Tns1 ,s2: 存放Simpson数列值Snc1 ,c2: 存放Cotes数列值CnR1 ,R2: 存放Romberg数列值Rnm: 存放复合梯形公式的节点加密次数er:存放误差注：为编程方便，程序中先算出了产生Romberg数列的值：T1、T2、T4、T8 、S1、S2、S4、C1、C2使，其中t=Table0,4提供一个4元素的表用于存放T1、T2、T4、T8。 4) 例题与实验例4.用Romberg求积方法求定积分的近似值,要求误差小于10-12。解：执行Romberg求积方法程序后，在输入的窗口中按提示分别输入：0、1、 4/(1+x2）、10(-12)每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出的计算结果如下：r(1)=3.141585783761874 eps=-8.31036401510-6r(2)=3.141592638396796 eps=6.85463492210-6r(3)=3.141592653590029 eps=1.51932333410-8r(4)=3.141592653589794 eps=-2.35811370410-13因此得到满足精度要求的定积分值为3.141592653589794 ，误差=-2.35811370410-134. Monte-Carlo求积方法 Monte-Carlo求积方法以概率论大数定理为依据式，借助随机数来求定积分的近似值的一种方法，该方法特别适用于求多重积分和积分区域复杂的重积分，不足之处是收敛较慢。1) Monte-Carlo求积方法的构造过程 假设重积分是有限的，式中WRn，G(x)是定义在W上的多元函数, P(x)是定义在W上的分布函数，则有只要式中x(k)是W上关于分布P(x)的n个任意独立选取的随机点列即可，这个结论可以有概率论中的大数定律： 以概率1成立。 利用以上结果求重积分的方法就称为Monte-Carlo求积方法。对具体的积分计算有如下计算公式：1） 求定积分的计算公式式中xk是a,b均匀分布的n个任意独立选取的随机数列。2） 求二重积分的计算公式式中(xk,yk)是 平面区域D中均匀分布的n个任意独立选取的随机点列，|D|表示区域D的面积。3） 求三重积分的计算公式式中(xk,yk,zk)是空间区域W中均匀分布的n个任意独立选取的随机点列，|W|表示区域W的体积。 Monte-Carlo求积方法的收敛阶为O(n1/2).2) Monte-Carlo求积方法算法1.输入被积函数f(x),积分区域边界和随机点的个数n2.在积分区域中产生n个随机点列3.计算被积函数f(x)在随机点列上的函数值并相加求平均4.用平均值与积分区域的测度之积作为积分近似值，这里测度就是长度、面积、体积之类的度量。3) Monte-Carlo求积方法程序 （*求定积分*） Cleara,b,x,n,f; a=Inputa= b=Inputb= fx_=Input被积函数f(x)= n=Input随机点个数n=;s=(b-a)*SumfRandomReal,a,b,n/nPrintn=,n, 积分,Ns,10;说明 本程序从用Monte-Carlo求积方法求a,b上的定积分近似值。程序执行后，按要求通过键盘输入积分下限a、积分上限b、被积函数f(x)和随机点个数n 后，计算机则给出定积分近似值。程序中变量说明a：存放积分下限b: 存放积分上限fx: 存放被积函数f(x)n: 存放随机点个数n（*求二重积分*） Cleara,b,x,y,n,f; a=Inputa= b=Inputb= c=Inputc= d=Inputd= fx_,y_=Input被积函数f(x,y)= n=Input随机点个数n=;s=(b-a)*(d-c)*SumfRandomReal,a,b, RandomReal,c,d,n/nPrintn=,n, 积分,Ns,10;说明 本程序从用Monte-Carlo求积方法求a,bc,d上的定积分近似值。程序执行后，按要求通过键盘输入积分变量x的下限a、积分上限b和积分变量y的下限c、积分上限d、被积函数f(x)和随机点个数n 后，计算机则给出定积分近似值。程序中变量说明a：存放积分变量x下限b: 存放积分变量x上限c：存放积分变量y下限d: 存放积分变量y上限fx: 存放被积函数f(x)n: 存放随机点个数n4) 例题与实验例5.用Monte-Carlo求积方法求定积分的近似值,随机点个数取为200。解：执行Monte-Carlo求积方法求定积分程序后，在输入的窗口中按提示分别输入：0、1、Sinx/x、200每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出的计算结果如下：n=200 积分=0.9463437539本题积分准确值为0.94608307036.例6.用Monte-Carlo求积方法求定积分的近似值,用随机点个数分别取为50，60，70，80，200的计算值进行计算实验来观察计算结果特点并提出一种从一组Monte-Carlo求积计算值推断的较好积分近似值的方法。解：修改Monte-Carlo求积方法求二重积分程序为：fx_,y_:=Exp-(x+y);Dos=SumfRandomReal,1,2,y=RandomReal,2,3, i,1,n,j,1,n/n2/N;Printn=,n, 积分=,Ns,10;,n,50,200,10执行后，计算机在屏幕上给出的计算结果如下：n=50 积分=0.0200165674n=60 积分=0.02009240135n=70 积分=0.0197868345n=80 积分=0.01983216355n=90 积分=0.01993588289n=100 积分=0.0199346424n=110 积分=0.01979465142n=120 积分=0.01987866262n=130 积分=0.01985729255n=140 积分=0.01984628761n=150 积分=0.01991800975n=160 积分=0.01985555608n=170 积分=0.01992147037n=180 积分=0.01982700398n=190 积分=0.01988543485n=200 积分=0.01989785087 本题积分准确值为0.0198937375894.,观察计算结果看到积分值并不是严格按照实验点数的增加而单调收敛的，因此有时对某个n值计算得出的结果可能没有比n小的值精确。但它的收敛性是肯定的，不过收敛速度很慢。 观察这16个计算结果值与准确值做比较，看到这样的特点：积分结果从小数点第二位后出现不同，但这16个数在小数点第二位有14个为1、在小数点第三位有14个为9，在小数点第四位有8个为8，它们在相应位数上取值最多，而这些值正好与准确值的相应数相同。因此，我们提出用统计一组计算数相应位数出现的数字规律确定积分结果的较好近似值，即，如果这组树的同一数字位中的某个数字出现的最多，则该数字应为准确值的有效数字，把所有有效数字组合在一起得到的树就是该积分值较好的近似值。7.5思考与提高1. 什么样的求积公式可以求出具有给定精度的定积分值？2.在给定求积精度e后，用复合梯形求积公式的余项来事先估计复合梯形求积公式的小区间数n与用判别式|T2n-Tn|e得出的小区间数n是否一样？用实验说明你的解答。3.分析n点Newto-Cotes求积程序的编程优缺点，你能给出哪些改进？4. 分析复合梯形求积程序的编程优缺点，你能给出哪些改进？5.分析Romber方法程序的编程优缺点，你能给出哪些改进？6.分析Monte-Carlo方法的编程优缺点，你能给出哪些改进？7.6 练习1分别用n=3，4，8的Newton-Cotes求积公式求定积分2求定积分的近似值,要求误差小于10-7。1） 用复合梯形公式2） 用复合抛物线求积公式3. 用Romberg求积方法求定积分的近似值,要求误差小于10-12。4 用Monte-Carlo求积方法求定积分求定积分的近似值。用随机点个数分别取为20，30，40，50进行计算。第8章 常微分方程初值问题数值解8.1 实验目的 了解常微分方程初值问题的数值解概念，掌握解常微分方程初值问题的Euler方法及改进的Euler方法和Runge-Kutta方法解常微分方程初值问题的算法构造和计算。能用程序实现解常微分方程初值问题的Euler方法、改进的Euler方法和经典的Runge-Kutta方法，学习用计算机求常微分方程初值问题数值解的一些科学计算方法和简单的编程技术。8.2 概念与结论1. 常微分方程初值问题 常微分方程特解问题称为初值问题，通常其形式为 2常微分方程初值问题数值解 常微分方程初值问题的解y(x)在a,b上的有限个值y(xk)的近似值yk称为常微分方程初值问题数值解，其中 xk= xk-1 + hk-1 ，k=1,2,3,N。xk称为节点， hk 称为步长。 通常，步长h不变，取为等距步长 hk=h=（b-a）/N，N为等分区间a,b分割数。3.常微分方程初值问题数值解法 求常微分方程初值问题数值解yk的方法称为微分方程初值问题数值解法。4.单步法 在计算yk+1之时只用到yk 的方法，其计算公式有： 显式单步法计算公式: yk+1=yk+h(xk,yk;h） 隐式单步法计算公式: yk+1 = yk +h(xk ,yk,yk+1,h） 式中函数是连续函数，称为增量函数。5.多步法 在计算yk+1之时不仅用到yk ，还要用yk-1，yk-2，,一般m 步法要用到yk,yk-1,yk-2， yk-m+1， 多步法也有显式方法和隐式方法之分。6.数值解法的局部截断误差 假设某常微分方程初值问题数值解法在x= xk没有误差，即y(xk)= yk，称Tk+1 =y(xk+1)- yk+1为该数值方法的局部截断误差。 显式单步法有其局部截断误差为Tk+1 =y(xk+1)- y(xk)- h(xk,y(xk),h）7.数值解法的阶 如果某常微分方程初值问题数值解法的局部截断误差为Tk+1 =O(h p+1)则称该数值方法的阶为P,或称该数值方法为P阶方法。 数值方法的阶越高，方法越好。8.3 程序中Mathematica语句解释 1. k+,+k,k-,-k k+ 表示赋值关系 k=k+1 。 例如, k=1;Table+k,5获得表2,3,4,5,6+k 表示先处理k的值,再做赋值 k=k+1。 例如, k=1;Tablek+,5获得表1,2,3,4,5k- 表示赋值关系 k = k-1。 例如，k=1;Tablek-,5获得表1, 0, -1, -2, -3-k 表示先处理k的值,再做赋值 k=k-1。 例如，k=1;Table-k,4获得表0,-1,-2,-32. x+=k, x*=kx+=k 表示 x = x + kx*=k 表示 x = x * k3.Forstat,test,incr,body 以stat为初值,重复计算incr和body直到test为False终止 。这里start为初始值,test为条件,incr为循环变量修正式,body为循环体,通常由incr项控制test的变化。 注意: 上述命令形式中的start可以是由复合表达式提供的多个初值,如果循环体生成 Break 语句,则退出For循环; 如果循环体生成Continue 语句,则由incr的增量进入For语句的下一次循环。 8.4 方法、程序、实验在实际问题中，常常会遇到一些常微分方程初值问题。对这些问题如果只求助于高等数学中介绍的用求通解加确定边界条件的方法去求解往往是无能为力的。因为一方面通解可能求不出来，另一方面所求的解可能是不能用初等函数表达的形式。人们处理这类问题主要采用常微分方程初值问题数值解的方法来求解。这类方法有很多，这里主要介绍解常微分方程初值问题的Euler方法、改进的Euler方法和Runge-Kutta方法的构造过程和算法程序。1. Euler方法Euler方法是最简单的求微分方程数值解的方法。这个方法由于精度不高，实用中较少使用。人们常用Euler方法来说明求微分方程数值解涉及到的一些问题。1) Euler方法的构造过程 Euler方法涉及的计算公式有很多方法，这里介绍在求微分方程数值解用的较多的Taylor展开法。 设 f(x,y) 充分光滑, xn+1=xn+h,将y(xn+1)在x n点作Taylor展开，得:y (x n+1)=y(xn)+hy(xn)+(h2/2!)y(xn)+O(h3)取其关于h 的线性部分，有：y (x n+1)y(xn)+hy(xn)注意到y(xn)= f(xn,y(xn),用yk代替y(xk)并将“”换为等号“=”,得到Euler公式：y n+1= yn+h f(xn,yn) 于是，由初始条件y0=y(a),借助Euler公式就可以依次计算出微分方程初值问题的数值解：y1 ,y2 ,yk称由Euler公式求数值解的方法为Euler方法。显然Euler方法是单步显式方法。 因为对Euler公式有局部截断误差为Tn+1= O(h2)因此Euler方法是一阶方法。2) Euler方法算法 1 输入函数f(x,y)、初值y0、自变量区间端点a,b步长h2 计算节点数n和节点x k 3 用Euler公式y n+1= yn+h f(xn,yn) 求数值解3) Euler方法程序Clearx,y,ffx_,y_= Input函数f(x,y)=y0=Input初值y0 =a=Input左端点a=b=Input右端点b=h=Input步长h=n=(b-a)/h;Fori=0,iInterpolatingFunctionrange, 的形式给出的,其中的InterpolatingFunctionrange, 是所求的插值函数表示的数值解, range就是所求数值解的自变量范围。4) 例题与实验例1. 用Euler方法求初值问题的数值解。取步长h=0.1,并在一个坐标系中画出数值解与准确解的图形。解：执行Euler方法程序后，在输入的窗口中按提示分别输入x-2x/y、1、0、1、0.1每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出如下数值解结果：y(0.1)=1.1y(0.2)=1.19182y(0.3)=1.27744y(0.4)=1.35821y(0.5)=1.43513y(0.6)=1.50897y(0.7)=1.58034y(0.8)=1.64978y(0.9)=1.71778y(1.)=1.78477本题的准确解为y(x)= (1 + 2 x)1/2，在一个坐标系中画出数值解与准确解的图形如下：图中点图是数值解，曲线为准确解。从图中可以知道数值解与准确解比较接近，但还是有误差的。2. 改进的Euler方法 改进的Euler方法比Euler方法精度高，其中把微分方程初值问题数值隐式解法计算公式变为预测、校正公式的做法很有代表性。1) 改进的Euler方法构造过程 将微分方程y=f(x,y) 在xk,xk+1积分，得对右端的定积分用数值积分的梯形公式，有用yk代替y(xk)并将“”换为等号“=”可得求微分方程初值问题的数值解的梯形公式：这个公式是单步隐式公式。可以证明它是二阶方法，比Euler方法阶高。不过，在给定初始条件y0=y(a)后,要求出数值解，每一次计算yk的值都要进行非线性方程求根的迭代解法来完成，因此计算量很大。为减少计算量，通常采用迭代一两次就转入下一步计算的方法，特别，如果先用Euler公式计算一次以对要计算的下一步值解进行预测，然后再用梯形公式对其进行校正的方法得到下一步的值，它的计算格式为：这个计算格式称为改进的Euler公式，而称由如上计算格式求数值解的方法为改进的Euler方法。2) 改进的Euler方法算法1.输入函数f(x,y)、初值y0、自变量区间端点a,b步长h2.计算节点数n和节点x k 3.用改进的Euler公式求数值解3) 改进的Euler方法程序Clearx,y,ffx_,y_= Input函数f(x,y)=y0=Input初值y0 =a=Input左端点a=b=Input右端点b=h=Input步长h=n=(b-a)/h;Fori=0,in,i+, xk=a+i*h; y1=y0+h*fxk,y0; y1=y0+h/2*(fxk,y0+fxk+h,y1); Printy(,xk+h/N,)=,y1/N; y0=y1注：语句h=Input步长h=的步长输入应该用带小数点的数输入，这样可以加快计算速度，否则，由于Mathematica软件本身的原因，可能会出现计算时间过长等计算不出结果的情况。说明 本程序用改进的Euler公式求常微方程初值问题数值解。程序执行后，按要求通过键盘依次输入输入函数f(x,y)、初值y0、自变量区间端点a,b步长h后，计算机则给出常微方程初值问题数值解。程序中变量说明fx,y: 存放函