1、 1计算物理学练习题及参考解答 1计算物理学的英文表示：computatioal physics 或者computer physics 2什么是计算物理学？它与理论物理、实验物理有什么区别和联系？ 答：计算物理是指以计算机及计算机技术为工具和手段，运用计算数学的方法解决复杂物理问题的一门应用科学。 计算物理方法是除理论方法和实验方法之外的第三种研究手段，计算物理现已成为物理学研究的三大支柱之一，它与实验物理和理论物理的关系如下图： 3计算物理学是物理学、数学、计算机科学三者结合的产物，它也是物理学的一个分支，与理论物理、实验物理有着密切的联系。 4计算机在物理学中有哪些应用？ 答：计算机数值分
2、析、计算机符号处理、计算机模拟、计算机实时控制 5计算机技术有各种各样的算法，可以概括为最基本的两类：串行计算和并行计算。 6理论物理在实际计算中遇到许多困难：非线性问题求解和非对称问题的求解；自变量较多问题求解；非规则界面问题求解等。 7计算物理的优点有：省时省钱；具有更大的自由度和灵活性；能够模拟极端条件下的实验。 8第一原理方法是基于量子力学基本原理建立起来的；分子动力学方法是基于经典力学基本原理建立起来的；蒙特卡罗方法是基于统计力学基本原理建立来的。 9计算机模拟一般有哪两种类型？ 答：随机模拟和确定性模拟，比如蒙特卡罗模拟和分子动力学模拟。 10什么是蒙特卡罗模拟？它的应用一般有哪三
3、种形式？ 答：通过不断产生随机数序列来模拟过程。直接蒙特卡罗模拟、蒙特卡罗积分、Metropolis蒙特卡罗模拟。 11蒙特卡洛方法的理论依据 答： （1）大数法则：人们发现，在一个随机事件中，随着试验次数的增加，事件发生的频率趋于一个稳定值；人们同时也发现，在对物理量的测量实践中，测定值的算术平均也具有稳定性。大数法则反映了大量随机数之和的性质。 （2）中心极限定理：中心极限定理，是概率论中讨论随机变量和的分布以正态分布为极限的一组定理。这组定理是数理统计学和误差分析的理论基础，指出了大量随机变量近似服从正态分布的条件。中心极限定理告诉我们：在有足够大，但又有限的抽样数n的情况下，蒙特卡洛估
4、计值是如何分布的。 12请简述著名的巴夫昂（Buffon）投针实验。并写出用 Matlab 实现的代码。（给出方程、算法框图、程序） 答：在平滑桌面上画一组相距为s的平行线，向此桌面随意地投掷长度 sl 的细针，那么从针与平行线相2交的概率就可以得到 的数值。 NMp ， p2M2N clear S=1; %平行线间距 L=1; %针长 N=1000000; %总投针次数 M=0; for i=1:N x=rand*S/2; %针到最近平行线的距离 a=rand*pi/2; %偏角 if(xsin(a)*L/2) M=M+1; %统计相交次数 end end testpi=2*N/M %pi的
5、实验值 13在考虑蒙特卡罗模拟的精确度时，不能只是简单地减少方差和增加模拟次数，还要同时兼顾计算费用，即机时耗费。 14假定我们研究连续的随机变量，由随机变量的分布可以得到它取某给定值的概率，即 )( duuuuPduug )(ug 称为u 的概率密度分布函数，它表示随机变量u取u 到 duu 之间值的概率。而 u dxxguG )()( 则称为u 的分布函数。 15高斯分布可以由给定的期望值 和方差 2 完全确定下来，通常用 ),N( 2 来表示 2/)(exp21),N( 222 x 比如期望值为1，方差为1的高斯分布表达式为 2/)1(exp21)1,1N( 2 x 16对物理问题的计算
6、机模拟所需要的伪随机数应当满足什么样的标准？有哪些统计检验方法？ 答：良好的统计分布特性；高效率；循环周期长；产生程序可以移植性好；可以重复产生。 统计检验有：均匀性检验；独立性检验；组合规律检验；无连贯性检验；参数检验等等。 17在蒙特卡洛方法应用中减小方差的基本技术：重要抽样法，分层抽样法，控制变量法和对偶变量法。然而，单独使用这四种减小方差的技巧仍然有其局限性。人们发展了一些复合蒙特卡洛计算技术，如适应性蒙特卡洛方法和多道蒙特卡洛抽样方法等。这些蒙特卡洛技巧对于被积函数在积分范围内具有多个尖峰的情况，特别具有实用价值。 18真随机数列 答：真随机数列是不可预测的，因而也是不能重复产生的数
7、据序列。 19伪随机数列 答：通过某些数学公式计算而产生的随机数列 20同余产生器及程序代码 答：一般通过如下的线性同余关系式来产生数列 3)(mod(x 1n mcaxn mxnn / float Random(int n, int m, float seed, float a, float b) int i; float r; r = seed; for (i = 1; i = n; i+) r = (a * r + b) % m; return r / m; 21均匀性检验 答：均匀性检验又称频率检验，它用来检验用随机数（样本值）确定的经验频率和均匀分布频率是否有显著性差异。常用的统计检
8、验方法有 2x 检验和累积频率检验（K-S检验）。 22随机变量抽样 答：指的就是由给定分布函数产生随机数的方法。首先，在0,1区间抽取均匀分布的伪随机数列，再从中抽取满足给定分布密度函数的简单子样，并且各个伪随机数相互独立。 23连续分布的随机变量抽样一般有哪些方法？ 答：直接抽样法；变换抽样法；舍选抽样法；复合抽样法；近似抽样方法 24试述离散型分布的随机变量的直接抽样 答：对于离散序列数 , 21 ixxx 给定每个数的取值概率 , 21 ippp ，则我们可定义其分布函数F(x) 如下： xxiipxF )( 。 在区间0,1上取均匀分布的随机数，判断满足下式的j值： )()( 1 j
9、j xFxF 则抽样值为 jx ，分布符合分布函数F(x)的要求为。 25、试述连续分布的随机变量的变换抽样法。 答：设连续型随机变量的分布密度函数为 )(xf 。要对满足分布密度函数 f(x)的随机变量 抽样较难时可考虑通过其它已知函数的抽样来得到。考虑变换 )(xhy ， )(ygx 选择 )(y ，使得 )()()()( xhxhdxdyyxf 则可对 )(y 抽样得到，通过变换 )( g ，得到满足分布密度函数 )(xf 的随机抽样。 4更为一般的情况，设连续型随机抽样 ),( 的分布密度函数为 ),( yxf 。考虑变换 ),(1 yxhu ， ),(2 yxhv ， Jyxhyxh
10、gJvugyxf ),(),(),(),( 21 ， yvxv yuxuJ / / ， 这样就可以通过抽取满足分布密度函数 ),( vug 随机抽样 ),( 得到待求的满足分布密度函数 ),( yxf 的随机抽样 ),( ，其中 ),(1 h ， ),(2 h 。 26试给出一个用随机数计算的Matlab程序。（10分） 解：物理模型： 如图第一项限中单位正方形内投点在圆内的概率即为单位圆面积的四分之一。 数学方程： 1010222121 )1(4 xxdxdx 算法框图：产生随机点（，）M个；统计其中满足条件 122 的点的个数N；计算值 MN/4 。 Matlab程序：P=4/100000
11、*length(find(sum(rand(2,100000).2)1) 27对指数分布 其他000,xef(x) - 直接抽样 答：第一步：求分布函数 xxx-e1dxef(x)dxF(x) 第二步：抽样 1,0 第三步：计算分布函数的反函数 -e-1)F( )1ln(1 28梅氏游走法计算氦原子中两电子间库伦作用的平均值。（给出动力学方程、数值解方程、算法框图、程序） 答：! 梅氏游走法计算氦原子中两电子间库伦作用的平均值 program main dimension x(6) !输入抽样参数Nt,Ng,Nf,Ns,dx Nt=1000 !热化步数 5Ng=100 !样本组数 Nf=10
12、!抽样间隔 Ns=10000 !每组抽样个数 dx=0.1 !游走步长 !初始化：随机设置初始值x do 10 i=1,6 call random(RND) 10 x(i)=0.01*(rnd-0.5) !两电子在核附近 !热化：消除初始化影响，趋于平衡分布 do 20 j=1,Nt 20 call walk(RND,dx,x) !分组间隔抽样，计算平均值和误差 su=0 su2=0 sdu=0 do 40 ig=1,Ng !样本分成Ng个组 ug=0 ug2=0 do 30 k=1,Ns !Ns间隔Nf抽样，Ns个样本为一组 call walk(RND,dx,x) if(mod(k,Nf).
13、ne.0) goto 30 x12=x(1)-x(4) y12=x(2)-x(5) z12=x(3)-x(6) r12=dist(x12,y12,z12) u=1/r12 ug=ug+u ug2=ug2+u*u !组内求和 30 continue ug=ug/Ns sigmag2=ug2/Ns-ug*ug dug=sqrt(sigmag2/Ns) !组内平均、方差和误差 su=su+ug su2=su2+ug*ug 40 sdu=sdu+dug !组间求和 avu=su/Ng sigma2=su2/Ng-avu*avu du1=sqrt(sigma2/Ng) du2=sdu/Ng del=du
14、1-du2 !组间平均、方差和误差 !输出avu,du1,du2,del 100 open(12,file=out.dat) write(12,1000) Nt,Ng,Nf,Ns,dx,avu,du1,du2,del close(12) 61000 format(4i10,5f15.4) end ! 计算距离的函数子程序 function dist(x,y,z) dist=sqrt(x*x+y*y+z*z) return end ! 计算权重的函数子程序 subroutine weight(x,f) dimension x(6) r1=dist(x(1),x(2),x(3) r2=dist(x
15、(4),x(5),x(6) f=exp(-3.375*(r1+r2) return end ! 梅氏游动一步的子程序 subroutine walk(RND,dx,x) dimension x(6),x0(6) call weight(x,f0) do 10 i=1,6 x0(i)=x(i) call random(RND) ! 存旧 10 x(i)=x(i)+dx*(RND-0.5) ! 生新 call weight(x,f) call random(RND) if(f.ge.f0*RND) goto 30 !游动 do 20 i=1,6 20 x(i)=x0(i) !不动 30 retur
16、n End 29有限差分法 答：微分方程和积分微分方程数值解的方法。基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组 ，解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。 30采用有限差分法求解微分方程时可以用直接法、随机游走法和迭代求解法。其中迭代法被广泛采用，有直接迭代法、高斯赛德尔迭代法和超松弛迭代法
17、。 )(41 ,2)( 1,)( 1,)( ,1)( ,11)(kji, jikjikjik jik ji qh 直接迭代式 )(41 ,2)1( 1,)1( 1,)( ,1)( ,11)(kji, jikjikjik jik ji qh 高斯赛德尔迭代式 7)4(4)1()(,2)1(1,)1(1,)(,1)(,1)(,)(,)1(,1)(kji,kjijikjikjikjikjikjikjikjiqh 超松弛迭代式 31有限元素方法 答：求解微分方程，特别是椭圆型边值问题的一种离散化方法，其基础是变分原理和剖分逼近。有限元方法是传统的里茨加廖金方法的发展，并融会了差分法的优点，处理上统一，
18、适应能力强，已广泛应用于科学与工程中庞大复杂的计算问题。 32采用有限元素法对场域进行三角形元素划分，下图所示划分是否允许（ ） yx0 33有限元素法的一般步骤 答：（1）推导出与给定边界条件的偏微分方程等价的泛函； （2）把求解的区域用三角形元素划分为小的单元，并对每个节点和三角形元素按照约定的规则进行编号； (3) 利用公式计算出各个三角形元素(e)的系数矩阵 eK)( 和 eP)( ； 将各个三角形单元的系数矩阵 eK)( 和 eP)( 装配成总矩阵 )(K 和 )(P ，形成有限元方程组，然后利用强加边界条件法对有限元方程组进行修正； (4) 最后，利用超松弛迭代法求解有限元方程组，
19、则得到区域内各个节点上的函数 值。 34有限元素法与有限差分法的比较 答：（1）数学方法上，有限元素法从泛函的变分开始，而有限差分法从物理模型的偏微分方程开始； （2）区域的离散化方法上，有限元素法采用三角形元素划分，而有限差分法一般采用矩形网络划分； （3）有限元素法采用统一的观点对区域内的节点和边界节点列出计算格式，而有限差分法则是孤立地对微分方程及定解条件分别列差分方程； （4）有限元素法要求计算机的内存量比较大，需要输入的数据量也比较大，而有限差分法的适用范围要广，特别是在边界形状比较规则时。 35关于分子动力学下面的各种说法中正确的有（ ） A分子动力学方法的英文名称是molecul
20、ar dynamics method； B分子动力学方法是按照体系内部的内禀动力学规律来计算并确定位形的转变； C分子动力学模拟方法可以看作是体系在一段时间内的发展过程的模拟，这个过程不存在任何随机因素； D原则上，分子动力学方法所适用的微观物理体系并无什么限制，既可以是少体系统，也可以是多体系统；既可以是点粒子体系，也可以是具有内部结构的体系；处理的微观客体既可以是分子，也可以是其他的微观粒子。 （1）分子动力学是在原子、分子水平上求解多体问题的重要的计算机模拟方法。 （2）通过求解所有粒子的运动方程，分子动力学方法可以用于模拟与原子运动路径相关的基本过程。 （3）在分子动力学中，粒子的运动
21、行为是通过经典的Newton运动方程所描述。 8（4）可以处理非平衡态问题。但是使用该方法的程序较复杂，计算量大，占内存也多。 （5）分子动力学方法中不存在任何随机因素。 36分子动力学 答：分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系综中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。 37采用分子动力学方法，必需对被计算的粒子系统给定适当的边界条件。这些边界条件大致可分成自由表面边界、固定边界、柔性边界、周期性边界四种。 38系统
22、达到平衡所需的时间称为弛豫时间。 39系综调节主要是指在进行分子动力学计算过程中，对温度和压力参数的调节，分为调温技术和调压技术。 40分子间作用力主要有长程的引力和短程的斥力。 411924 年，Lennard-Jones发表了著名的负幂函数式的势函数的解析式为 )r(-)r(4U(r) 612 。 42系综(ensemble)是指在一定的宏观条件下(约束条件)，大量性质和结构完全相同的、处于各种运动状态的、各自独立的系统的集合。 43、请写出分子动力学模拟的实际步骤。 答：分子动力学模拟的实际步骤可以划分为四步： （1）设定模拟所采用的模型模型的设定，也就是势函数的选取。势函数的研究和物理
23、系统上对物质的描述研究息息相关。 （2）给定初始条件。运动方程的求解需要知道粒子的初始位置和速度，不同的算法要求不同的初始条件。如：verlet算法需要两组坐标来启动计算，一组零时刻的坐标，一组是前进一个时间步的坐标或者一组零时刻的速度值。 （3）趋于平衡的计算过程。为使得系统平衡，模拟中设计一个趋衡过程，即在这个过程中，增加或者从系统中移出能量，直到持续给出确定的能量值。称这时的系统已经达到平衡。 （4）宏观物理量的计算。实际计算宏观物理量往往是在分子动力学模拟的最后阶段进行的。它是沿着相空间轨迹求平均来计算得到的。 44实际上，分子动力学模拟方法和随机模拟方法一样都面临着两个基本限制： (
24、1)有限观测时间的限制； (2)有限系统大小的限制。 通常人们感兴趣的是体系在热力学极限下（即粒子数目趋于无穷时）的性质。但是计算机模拟允许的体系大小要比热力学极限小得多，因此可能会出现有限尺寸效应。为了减小有限尺寸效应，人们往往引入周期性、全反射、漫反射等边界条件。当然边界条件的引入显然会影响体系的某些性质。 45请写出实现周期性边界条件的具体操作 答：当有一个粒子穿过基本原胞的六方体表面时；就让这个粒子也以相同的速度穿过此表面对面的表面，重新进入该原胞内。 46在考虑粒子间的相互作用时，通常采用的最小像力约定是什么？。 答：最小像力约定是：设原胞内有 N个粒子，在重复的分子动力学原胞中，每
25、一个粒子只同它所在的原胞内的另外 N-1 个中的每个粒子或其最邻近的影像粒子发生相互作用。 47采用最小像力约定会使得在截断处粒子的受力有一个-函数的奇异性，这会给模拟计算带来误差。怎样减少这种误差？ 答：为减小这种误差，总可以将截断处粒子的相互作用势能换算成 V(r)- Vc ，以保证在截断处相互作用为零。 48EAM 模型是一种半经验多体势模型，其中心思想是什么？ 9答：其中心思想是将原子周围复杂的环境用胶冻模型简化描述（所谓胶冻指的是在均匀正电荷背景上的均匀电子气）。在嵌入原子势方法中，主要的参数都放在原子的电子密度表示及相关形式中。这样就把 “原子对间”的性质主要归结到“原子”的性质，
26、大大的简化了计算。 49分子动力学模拟的时间步长如何选择？ 答：分子动力学计算的基本思想是赋予分子体系初始运动状态之后利用分子的自然运动在相空间中抽取样本进行统计计算，时间步长就是抽样的间隔，因而时间步长的选取 对动力学模拟非常重要。太长的时间步长会造成分子间的激烈碰撞，体系数据溢出；太短的时间步长会降低模拟过程搜索相空间的能力，一般来说，时间步长应该设为分子运动的最小振动周期的1/10 左右为宜。 50Verlet 算法中分子动力学计算的简单步骤是什么？ 答：Verlet 算法，速度项相消，计算坐标时无需速度出现。 简单步骤： 设定粒子的初始位置和速度； 根据粒子的位置计算每个粒子的受力；
27、根据粒子的位置、速度和受力，计算粒子的新位置和新速度； 更新粒子的位置和速度，然后回到 2 步骤。 51简述分子动力学模拟步骤。 答：第一步 模型的设定。 也就是势函数的选取。势函数的研究和物理系统上对物质的描述研究息息相关。最早是硬球势，即小于临界值时无穷大，大于等于临界值时为零。常用的是 LJ势函数，还有EAM势函数，不同的物质状态描述用不同的势函数。模型势函数一旦确定，就可以根据物理学规律求得模拟中的守恒量。 第二步 给定初始条件。 运动方程的求解需要知道粒子的初始位置和速度，不同的算法要求不同的初始条件。如：verlet算法需要两组坐标来启动计算，一组零时刻的坐标，一组是前进一个时间步
28、的坐标或者一组零时刻的速度值。 第三步 趋于平衡计算。 在边界条件和初始条件给定后就可以解运动方程，进行分子动力学模拟。但这样计算出的系统是不会具有所要求的系统的能量，并且这个状态本身也还不是一个平衡态。 为使得系统平衡，模拟中设计一个趋衡过程，即在这个过程中，我们增加或者从系统中移出能量，直到持续给出确定的能量值。我们称这时的系统已经达到平衡。这段达到平衡的时间成为驰豫时间。 第四步 宏观物理量的计算。 实际计算宏观的物理量往往是在模拟的最后阶段进行的。它是沿相空间轨迹求平均来计算得到的。 52神经网络 答：人工神经网络是一种应用类似于大脑神经突触联接的结构进行信息处理的数学模型。在工程与学
29、术界也常直接简称为“神经网络”或类神经网络。神经网络是一种运算模型，由大量的节点（或称“神经元”，或“单元”）和之间相互联接构成。每个节点代表一种特定的输出函数，称为激励函数（activation function）。每两个节点间的连接都代表一个对于通过该连接信号的加权值，称之为权重（weight），这相当于人工神经网络的记忆。网络的输出则依网络的连接方式，权重值和激励函数的不同而不同。而网络自身通常都是对自然界某种算法或者函数的逼近，也可能是对一种逻辑策略的表达。 53计算机仿真 答：利用计算机科学和技术的成果建立被仿真的系统的模型，并在某些实验条件下对模型进行动态实验的一门综合性技术。它具
30、有高效、安全、受环境条件的约束较少、可改变时间比例尺等优点，已成为分析、设计、运行、评价、培训系统（尤其是复杂系统）的重要工具。 54并行计算有什么优点？ 答案要点：(1) 并行计算可以大大加快运算速度，即在更短的时间内完成相同的计算量，或解决原来根本不能计算的非常复杂的问题。 (2) 提高传统的计算机的计算速度一方面受到物理上光速极限和量子效应的限制，另一方面计算机器件产品和材料的生产受到加工工艺的限制，其尺寸不可能做得无限小。因此我们只能转向并行算法。(3) 并行计算对设备的投入较低，既可以节省开支又能完成计算任务。并行计算(Parallel Computing)是指同时使用多台计算机或一
31、台计算机多个处理器协同合作解决计算问题的过程，其主要目的是快速解决大型且复杂的计算问题。 55Matlab中有专门求解偏微分方程的工具箱，打开此工具箱的命令为 pdetool 1056.下面哪一图是边值问题xxyxuyxuyyxuyxuyxuyycos3sin)1,(,0)0,(3sin),1(,0),0(1,0,0uxx的解 A B 57用 simulink 仿真自由落体运动过程 mxgdtdx 11 ，12 xdtdx ，其中 速度:1x ， 位移:2x 。g=9.8，m=2， -1 。 58用 matlab仿真乒乓球弹跳过程（给出动力学方程、数值解方程、算法框图、程序） %pingpan
32、gstatefun.m function xdot=pingpangstatefun(t,x,flag) %乒乓球弹跳模型的标准状态方程 %x(1)为小球速度，x(2)为小球位移 g=9.8; %重力加速度 m=2; %乒乓球质量 beta=-1; %空气阻力系数 rho=1.29; %空气密度 V=1; %落体的体积 xdot=zeros(2,1); %状态变量矩阵初始化 xdot(1)=-g+beta*x(1)/m+rho*V*g/m;%速度加速度方程 11xdot(2)=x(1); %位移速度方程 %main.m clear; v0=0;% 球的初始态 y0=1; x_state=v0,
33、y0;% 球的初始状态赋值到状态变量中 dt=0.01;% 仿真步进 t=0:dt:5;% 仿真时间序列 K=1;% 碰撞衰减系数 x=zeros(length(t),length(x_state);% for k=1:length(t)% 仿真开始，每次循环向前推进一个仿真步进dt x(k,:)=x_state;% 记录并保存当前状态的计算结果 t_out,x_out=ode45(pingpangstatefun,t(k),t(k)+dt,x_state);% 计算下一时刻的新状态 x_state=x_out(length(x_out),:);% 更新状态 if(x_state(2)=0)
34、% 处理碰撞瞬间情况：速度反向并衰减K end %动画作图：显示小球弹跳过程 y=x_state(2);% 小球当前位置 subplot(2,1,1); plot(0,y,.,MarkerSize,40,Color,black); axis(-2 2 0 1);% 坐标范围固定 set(gcf,DoubleBuffer,on);% 双缓冲避免作图闪烁 drawnow;% 立即显示作图 end %仿真结束，最后输出计算时间序列上的结果：随时间变化的速度和位移曲线 subplot(2,1,2); plotyy(t,x(:,1),t,x(:,2),plot); 补充： 1按四舍五入取的近似值3.14
35、，则其相对误差限为 解：按四舍五入取近似值=3.14，其绝对误差限为=0.005， 0.159% 2=3.1415926，取3.14做为近似值，误差为，有3位有效数字；若取3.141做为近似值，误差为，仍然是3位有效数字；若取3.142做为近似值，误差为有4位有效数字。 3误差来源有哪几种？ （1）. 模型误差（Modeling Error）： （2）. 观测误差（Measurement Error）： （3）. 截断误差（Truncation Error）： （4）舍入误差（Roundoff Error） 4线性方程组的数值解法一般有哪几种？ 答：（1）高斯消去法；（2）LU分解法；（3）三
36、对角矩阵追赶法；（4）迭代法 5非线性方程的数值解法一般有：（1）二分法；（2）弦截法；（3）不动点迭代法；（4）牛顿迭代法 6常微分方程 ),()( utftu 对应的显式 Euler 公式为 ),(1 kkkk uthfuu ；对应的隐式 Euler 公式为12),( 111 kkkk uthfuu ；对应的两步法Euler公式为 ),(211 kkkk uthfuu 7数值微分是用邻近一些点上的函数值近似表示函数在一个点上的各阶导数值，Rung-Kutta 方法就是基于这一思路，用 )(,( tutf 在一些特殊点上函数值的线性组合来表示泰勒展开中某点的各阶导数值。最常用的四阶RK公式为：
