自适应变步长占空比扰动法在光伏发电MPPT中的应用_薛阳
短路电流结合扰动观察法在光伏发电最大功率点跟踪控制中的应用
第26卷第20期中国电机工程学报V ol.26 No.20 Oct. 2006 2006年10月Proceedings of the CSEE ?2006 Chin.Soc.for Elec.Eng. 文章编号：0258-8013 (2006) 20-0098-05 中图分类号：TM51 文献标识码：A 学科分类号：480?60 短路电流结合扰动观察法在光伏发电 最大功率点跟踪控制中的应用 张超，何湘宁 （浙江大学电力电子研究所，浙江省杭州市310027） Short-current Combined With Perturbation and Observation Maximum-power-point Tracking Method for Photovoltaic Power Systems ZHANG Chao, HE Xiang-ning (Power Electronics Institute of Zhejiang University, Hangzhou 310027, Zhejiang Province, China) ABSTRACT: The output power of PV module varies with module temperature, solar insolation and loads, so it is necessary to track MPP of the PV array all the time. In past years, many MPP control algorithms were presented to draw maximum power from the solar array. A novel online short circuit current method is presented. This method can track MPP changes rapidly without disturbing PV system. On the basis of this method, P&O(perturbation and observation) method with optimized perturbation step was proposed to reduce the power oscillation around MPP. Simulations and experimental results show that the PV generation system has good steady state and transient characteristics with the proposed MPPT control method. KEY WORDS：photovoltaic; maximum power point tracking; short circuit control; perturbation and observation control 摘要：光伏电池输出功率随外部环境和负载的变化而变化，为充分发挥光伏器件的效能，需采用最大功率点跟踪电路。对于最大功率点跟踪电路的控制已经提出了许多方法，其中短路电流法和扰动观察法因其具有简单有效的优点而得到广泛应用。针对短路电流法的缺点，该文提出一种新的在线短路电流控制方法。该方法在不干扰系统正常工作的情况下，能迅速感知外部环境变化，但该方法效率不高。为充分发挥光伏电池的效能，在线短路电流控制方法的基础上再引入扰动观察法。该文扰动观察法的扰动步长针对最大功率点处稳态特性进行优化，优化后，扰动观察法可有效消除光伏器件输出功率在最大功率点的振荡现象，从而提高系统效率。仿真和实验研究证明，该方法可以快速跟踪外部环境变化，并消除系统在最大功率点的振荡现象。 基金项目：国家教育部博士点基金项目(20050335059)．关键词：光伏；最大功率点跟踪；短路电流法；扰动观察法0 引言 光伏发电作为一种具有广阔前景的绿色能源已成为国、内外学术界和工业界研究的热点[1-2]。光伏电池输出功率与外界环境和负载情况有关，为充分发挥光伏电池的功效，需在光伏器件和负载之间串联最大功率点跟踪(maximum power point tracking，MPPT)电路[3-4]。MPPT电路常用的控制方法有固定参数法(如固定电压法、固定电流法)、扰动观察法及增加电导法等[5]。 固定参数法利用在最大功率点工作时光伏器件工作电压、电流与器件开路电压、短路电流的近似比例关系进行控制，此方法只需一个检测参数，控制简单易行，但获取开路电压或短路电流要中断系统正常工作，对系统运行存在干扰，此外所采用的控制关系是近似关系，不能实现最优控制，因此该方法控制精度低，仅适用于小功率场合。扰动观察法根据光伏器件在最大功率点处?P/?U=0的特性进行最大功率点跟踪控制，以左侧为例说明该方法具体工作过程：在系统稳定工作情况下，假设增大最大功率点跟踪电路功率器件的占空比，控制器对占空比调节前后的光伏器件输出功率、输出电压进行采样计算，如果输出功率与输出电压为?P/?U >0，则表明系统工作在最大功率点左侧，应继续增加占空比，直到?P/?U =0，具体工作中，由于扰动观测法需要比较占空比变化前后的功率、电压，因此光伏输出功率会在最大功率点两次反复变化，严重时
龙格库塔积分算法
龙格库塔法 龙格库塔法是常用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家C. Runge和M.W. Kutta于1900年左右发明。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。 龙格库塔法是一种在工程上应用广泛的高精度单步算法，可以应用在物理、工程、控制、动力学中，如模糊控制、弹道分析以及分析光纤特性等，在系统仿真中得到广泛应用。 龙格库塔法源自于相应的泰勒级数方法，在每一插值节点用泰勒级数展开,其截断误差阶数也是，根据可省略更高阶的导数计算, 这种方法可构造任意阶数的龙格库塔法。其中4 阶龙格库塔法是最常用的一种方法。因为它相当精确、稳定、容易编程。在计算中一般不必使用高阶方法, 因为附加的计算误差可由增加精度来弥补。如果需要较高的精度, 可采取减小步长的方法即可。4 阶龙格库塔法的精度类似4 阶泰勒级数法的精度。 1、初值问题 对于一阶常微分方程的初值问题 根据常微分方程的理论可知，此初值问题的解在区间[a，b]上存在，且唯一。 2、离散化
取步长h=(b-a)/n，将区间[a , b]分成n个子区间： a=<=b 在其中任意两点的曲线段上，根据积分中值定理，一段光滑曲 线上至少有一点，它的斜率与整段曲线的平均斜率相同， 得=y’() (0<<1) 其中，= 可以将上式改写成y()=y()+h*K (2.1) 其中K为平均斜率，K=f() 公式（2.1）表明，如果能够确定平均斜率K，就可以根据（2.1）式得到y()的值。 欧拉法和龙格库塔法就是用不同方法确定不同精度的平均斜率K，从而求得y()的近似值。 3、Euler法 欧拉法虽然精度低，但它是最简单的一种显式单步法，也是龙 格库塔法的基础。 首先，令、为y() 及y()的近似值，并且令平均斜 率K=f()，即以点的斜率作为平均斜率K，便得到欧拉公式=+h* f() （3.1） 4、改进的欧拉法 此种方法是取、两点的斜率的平均值作为平均斜率K， 即K= ,其中、均为y()以及y()的近似值，就得到 改进后的欧拉公式（4.1）
扰动观察法
function [sys,x0,str,ts,simStateCompliance] = PO_MPPT_Boost(t,x,u,flag) switch flag, case 0, [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes; case 1, sys=mdlDerivatives(t,x,u); case 2, sys=mdlUpdate(t,x,u); case 3, sys=mdlOutputs(t,x,u); case 4, sys=mdlGetTimeOfNextVarHit(t,x,u); case 9, sys=mdlTerminate(t,x,u); otherwise DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag)); end function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes sizes = simsizes; sizes.NumContStates = 0; sizes.NumDiscStates = 3; sizes.NumOutputs = 1; sizes.NumInputs = 2; sizes.DirFeedthrough = 1; sizes.NumSampleTimes = 1; % at least one sample time is needed
sys = simsizes(sizes); x0 = [0,0,0.5];%[U I Uref] str = []; ts = 1e-6; simStateCompliance = 'UnknownSimState'; % end mdlInitializeSizes function sys=mdlDerivatives(t,x,u) sys = []; % end mdlDerivatives function sys=mdlUpdate(t,x,u) DU=0.001; %步长 dU=u(1)-x(1); dP=u(1)*u(2)-x(1)*x(2); if dP>0 if dU>0 Uref=x(3)+DU;
一种应用于光伏系统MPPT的变步长扰动观察法
万方数据
万方数据
万方数据
一种应用于光伏系统MPPT的变步长扰动观察法 作者：朱铭炼， 李臣松， 陈新， 龚春英， ZHU Ming-lian， LI Chen-song， CHEN Xin，GONG Chun-ying 作者单位：南京航空航天大学,江苏,南京,210016 刊名： 电力电子技术 英文刊名：POWER ELECTRONICS 年，卷(期)：2010,44(1) 参考文献(5条) 1.欧阳名三;余世杰;沈玉樑一种太阳能电池MPPT控制器实现及测试方法的研究[期刊论文]-电子测量与仪器学报2004(02) 2.Wen-Jung Chiang;Hurng-Liahng Jou;Jinn-Chang Wu Maximum Power Point Tracking Method for the Voltage-mode Grid-connected Inverter of Photovoltaic Generation System 2008 3.LIU Fang-rui;DUAN Shan-xu;LIU Fei A Variable Step Size INC MPPT Method for PV Systems 2008(07) 4.张超;何湘宁短路电流结合扰动观察法在光伏发电最大功率点跟踪控制中的应用[期刊论文]-中国电机工程学报2006(20) 5.徐鹏威;刘飞;段善旭几种光伏系统MPPT方法的分析比较及改进[期刊论文]-电力电子技术 2007(05) 本文链接：https://www.wendangku.net/doc/e717424817.html,/Periodical_dldzjs201001008.aspx
固定电流法与变步长扰动观察法结合的M算法研究
2012～2013学年度电气与电子工程学院研究生课程 《太阳能光伏并网发电系统》课程报告 固定电流法与变步长扰动观察法结合的 MPPT算法研究 院系：电气与电子工程学院 专业：应用电子工程系 任课教师： 学生姓名： 学号： 指导教师： 二○一三年五月
ABSTRACT:The output power of PV module is a non-linear function of temperature, solar insolation and loads, so it is necessary to track MPP of the PV array all the time. When the external environment rapidly changes，in order to regulate the PV array operating point near the maximum power point quickly, the online constant-current method to track MPP is utilized at first, and then the variable step P&O is adopted to adjust the PV array operating point to be at the maximum-power point. The simulation of the improved method and common methods show that, contrast to the latter, the former can trace the maximum-power point more quickly, efficiently and accurately, which is concise and easy to implement, and also can reduce the energy loss caused by the oscillation of the operating point about the maximum-power point, thus enhancing the PV system efficiency. KEYWORDS:Photovoltaic panel characteristics; maximum power point tracking; constant current method; variable perturb step perturbation and observation method 摘要：光伏器件输出功率是外部环境、负载的非线性函数，为了充分发挥光伏器件的效能，需采用最大功率点跟踪方法进行控制。当外界环境突变时，采用在线固定电流法进行初步跟踪，调整光伏阵列的工作点到最大功率点附近。在此基础上再使用变步长扰动观察法，使得工作点进一步调节到最大功率点，并有效减少了光伏阵列输出功率在最大功率点的振荡。对该结合方法及相关的MPPT算法分别仿真，结果表明，该方法可以在外界环境剧烈变化下快速、有效、准确地跟踪最大功率点，简明易于实现，同时有效减少在光伏阵列最大功率点附近振荡所带来的能量损失，提高了光伏发电系统的效率。 关键词：光伏电池特性；最大功率点跟踪；固定电流法；变步长扰动观察法 1引言 光伏电池是光伏发电系统电能的来源，光伏电池输出功率是其所受日照强度、器件结温的非线性函数。即使在外部环境稳定的情况下，光伏电池的输出功率也会随着外部负载的变化而变化，只有当外部负载与光伏器件达到阻抗匹配时，光伏器件才会输出最大功率。为了实时从光伏阵列获得最大输出功率，需要在光伏发电系统中实现最大功率点的跟踪控制。通常的实现方法是需根据外部环境和负载情况不断调节光伏器件的工作点使其输出最大功率，我们将此功率调节过程称为最大功率点跟踪(Maximum Power Point Tracking，MPPT)。 MPPT电路常用的控制方法有固定参数法（如固定电压法、固定电流法）、扰动观察法及增加电导法等。各种MPPT控制方法中，固定电流法和扰动观察法因简单有效而较常用，但各自也存在缺陷。本文在分析固定电流法和扰动观察法的基础上，采用了一种固定电流法结合扰动观察法的MPPT控制。在外界环境或负载突变时，采用在线固定电流法将光伏阵列的工作点调整到最大功率点附近，以保证跟踪的快速性。在此基础上，为进一步提高对光伏阵列的利用效率，在最大功率点附近，采用变步长扰动观察法，从而减小系统在最大功率点附近的振荡。 2光伏电池特性 光伏电池单元是非线性器件，它的等效电路模型如图2.1所示。图2.1中的电流源产生光生电流I ph，它的数值由光照强度与温度共同决定。串联电阻R S与并联电阻R P用来表征太阳能电池内部的功率损耗。由于太阳能电池表面的材料的电阻率，当电流经过太阳能板流向外部内路时，其必然为产生串联损耗，故引入串联电阻R S。并联电阻R P表征由漏电流引起的损耗。
光伏系统变步长扰动观察法MPPT算法研究
光伏系统变步长扰动观察法MPPT算法研究 王小昆，胡贤新 (中国矿业大学信电学院，江苏，徐州221116) 摘要：光伏电池阵列输出功率受光照强度和温度变化的影响，因此最大功率点跟踪（MPPT）技术广泛应用于光伏系统中。在所有最大功率点(MPP)控制策略中，扰动观察（P&O）MPPT算法因易实现被广泛应用，然而它的缺点是在稳定工作状态下工作点通过MPP时会导致能量振荡损耗，并且在光照强度或温度发生突变时表现较差的动态响应。在本文中，提出一种改进型变步长扰动观察MPPT算法，此方法依据工作点动态调整步长变化，与传统固定步长方法比较，本文提出的方法能有效地提高MPPT速度和转换效率，通过仿真和实验结果分析，验证了此改进算法的可行性。 关键词：最大功率点跟踪；扰动观察法；变步长 Research on Variable Step Size P&O MPPT Algorithm for PV Systems WANG Xiao-kun，HU Xian-xin (School of Information and Electrical Engineering, China University of Mining and Technology, Xuzhou, Jiangsu 221116, China) Abstract:Maximum power point tracking (MPPT) techniques are employed in photovoltaic (PV) systems to maximum the PV array output power which depends on solar irradiance and temperature. Among all the MPPT strategies,The P&O Maximum Power Point Tracking algorithm is mostly used, due to its ease of implementation. However, its main drawbacks are the waste of energy in steady conditions, when the working point moves across the MPP and the poor dynamic performances exhibited when a step change in solar irradiance or in temperature occurs. In this paper, a modified variable step size P&O MPPT algorithm is proposed, the step size is automatically tuned according to the operating point. Compared with the conventional fixed step size method, the proposed approach can effectively improve the MPPT speed and efficiency simultaneously.A theoretical analysis and the design principle of the proposed algorithm are provided and its feasibility is also verified by simulation and experimental results. Keywords：Maximum power point tracking (MPPT); the Perturb and Observe (P&O); Variable Step Size 1 引言 随着世界能源消耗的不断增加，石油、煤炭等不可再生能源不断枯竭，可再生能源在能源消费中的地位越来越突出，其取之不尽和对自然环境影响较小的特点，得到了广泛的应用，其中太阳能是主要清洁能源之一，主要是光伏发电系统投资成本越来越低，并且技术不断进步、系统运行稳定、维护简单、对环境无污染。光伏发电系统现在主要应用于蓄电池充电、扬水系统、居民生活供电、卫星系统供电等场所[1-2]。 定稿日期：2010-12-11 作者简介：王小昆（1983-），男，安徽巢湖人，硕士研究生，研究方向光伏并网系统及电力电子技术 光伏系统存在两个重要不足，其一电池板发电效率非常低，在低光照强度下表现尤为突出，其二电能的产生随着天气状况变化而产生明显的变化，比如光照强度和温度变化等。 最大功率点跟踪（Maximum power point tracking，简称MPPT）变成了光伏发电系统中的重要组成部分，它能将电池板阵列产生的最大功率传输到逆变系统中，实现效率的最优化。基于控制器复杂程度，所用传感器数量，制造成本和性能有多种不同的MPPT控制方法被相继提出。MPPT算法实现目标是跟踪的快速性和准确性，即在光照强度和温度发生突变时快速响应以及稳定工作时较小振荡。目前常用的方法有恒压法
5.2龙格库塔法
第五章 常微分方程的数值解法 5.2 龙格-库塔法 一、教学目标及基本要求 通过对本节课的学习，使学生掌握常微分方程、常微分方程方程组的数值解法。 二、教学内容及学时分配 本节课主要介绍常微分方程的数值解法。具体内容如下： 讲授内容：龙格库塔方法、收敛性与稳定性 三、教学重点难点 1．教学重点：龙格库塔方法、收敛性与稳定性。 2. 教学难点：收敛性与稳定性。 四、教学中应注意的问题 多媒体课堂教学为主。适当提问，加深学生对概念的理解。 五、正文 龙格－库塔方法 引言 龙格－库塔方法的基本思想 '11()() ()()()(,())n n n n y x y x y y x y x hf y h ξξξ++-=?=+ 令*(,())K f y ξξ=称为区间1[,]n n x x +上的平均斜率，只要对*K 提供一种算法，即可推导出一种计算公式。 欧拉公式只是取点n x 的斜率作为区间1[,]n n x x +的平均斜率*K ，精度自然很低。 考察改进的欧拉公式 11211 ()22 n n y y h k k +=++1(,)n n k f x y =21(,)n n k f x h y hk =++ 它利用n x ，1n x +两点的斜率取算术平均，1n x +处斜率通过已知信息n y 用欧拉公式预报得到。
可以考虑设法在1[,]n n x x +多取几个点的斜率值，将它们加权平均作为区间 1[,]n n x x +的平均斜率*K 。这就是龙格－库塔方法的基本思想。 1、二阶龙格－库塔方法 考察1[,]n n x x +内一点,01n p n x x hp p +=+<≤，希望通过,n n p x x +两个点的斜率 12,K K 加权平均得到*K ，即令112((1))n n y y h K K λλ+=+-+ 取1(,)n n K f x y =，如何预报n p x +处斜率2K ？ 仿照改进的欧拉公式，先用欧拉公式预测()n p y x +的值n p y +： 1n p n y y phK +=+ 然后用n p y +计算2(,)n p n p K f x y ++=，从而得 112121[(1)](,)(,) n n n n n p n y y h K K K f x y K f x y phK λλ++=+-+==+ 适当选取,p λ，使上述公式具有较高得精度。假定()n n y y x =，分别将12 ,K K 泰勒展开： '1(,)()n n n K f x y y x == 212 ' '' 2 (,) (,)[(,)(,)(,)]()()()() n p n n n x n n n n y n n n n K f x y pkK f x y ph f x y f x y f x y O h y x phy x O h +=+=+++=++ 代入得 '2''31()()()()n n n n y y x hy x ph y x O h λ+=+++ 按泰勒展开2' '' 31()()()()2 n n n n h y y x hy x y x O h +=+++ 比较得，只要1 2 p λ= ，公式截断误差为3()O h 特别，当1,1/2p λ==，就是改进的欧拉公式， 改取1,1/2p λ==， 12n n y y hk +=+，1(,)n n k f x y =，21(,)22 n n h h k f x y k =++
龙格库塔法例题
四阶龙格一库塔法 通常所说的龙格一库塔法是指四阶而言的．我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格一库塔公式 (9.22) 公式(9.22)的局部截断误差的阶为． 龙格一库塔法具有精度高，收敛，稳定(在一定的条件下)，计算过程中可以改变步长，不需要计算高阶导数值等优点．但仍需计算在一些点上的值，如四阶龙格－库塔法每计算一步需要算四次 的值，这给实际计算带来一定的复杂性，因此，多用来计算“表头”．(即开始几点的近似值)．例3.用标准龙格一库塔法求初值问题 在处的解． 解因与．若应用标准龙格一库塔方法公式(9.22)计算，对于n=0时，则有
于是得 这个值与准确解在处的值已十分接近．再对n=1，2，3，4应用式(9.22)计算，具体计算结果如表3所示：
例3写出用四阶龙格――库塔法求解初值问题 的计算公式，取步长h＝0.2计算y(0.4)的近似值。至少保留四位小数。 解此处f(x，y)＝8－3y，四阶龙格――库塔法公式为 其中κ1＝f(x k，y k)；κ2＝f(x k+0.5h，y k+0.5hκ1)； κ3＝f(x k+0.5h，y k+0.5hκ2)；κ4＝f(x k+h，y k+hκ3) 本例计算公式为： 其中κ1＝8－3y k；κ2＝5.6－2.1y k； κ3＝6.32－2.37y k；κ4＝4.208－1.578y k ＝1.2016＋0.5494y k (k＝0，1，2，…) 当x0＝0，y0＝2， y(0.2)≈y1＝1.2016＋0.5494y0＝1.2016＋0.5494×2＝2.3004 y(0.4)≈y2＝1.2016＋0.5494y1＝1.2016＋0.5494×2.3004＝2.4654
常微分方程组的四阶RungeKutta龙格库塔法matlab实现
常微分方程组的四阶Runge-Kutta方法1.问题： 1.1若用普通方法-----仅适用于两个方程组成的方程组 编程实现： 创建M 文件： function R = rk4(f,g,a,b,xa,ya,N) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here %x'=f(t,x,y) y'=g(t,x,y) %N为迭代次数 %h为步长 %ya,xa为初值 f=@(t,x,y)(2*x-0.02*x*y);
g=@(t,x,y)(0.0002*x*y-0.8*y); h=(b-a)/N; T=zeros(1,N+1); X=zeros(1,N+1); Y=zeros(1,N+1); T=a:h:b; X(1)=xa; Y(1)=ya; for j=1:N f1=feval(f,T(j),X(j),Y(j)); g1=feval(g,T(j),X(j),Y(j)); f2=feval(f,T(j)+h/2,X(j)+h/2*f1,Y(j)+g1/2); g2=feval(g,T(j)+h/2,X(j)+h/2*f1,Y(j)+h/2*g1); f3=feval(f,T(j)+h/2,X(j)+h/2*f2,Y(j)+h*g2/2); g3=feval(g,T(j)+h/2,X(j)+h/2*f2,Y(j)+h/2*g2); f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3); g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3); X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6; Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6; R=[T' X' Y']; end 情况一：对于x0=3000，y0=120 控制台中输入： >> rk4('f','g',0,10,3000,120,10) 运行结果： ans = 1.0e+003 * 0 3.0000 0.1200 0.0010 2.6637 0.0926 0.0020 3.7120 0.0774 0.0030 5.5033 0.0886 0.0040 4.9866 0.1193 0.0050 3.1930 0.1195 0.0060 2.7665 0.0951 0.0070 3.6543 0.0799 0.0080 5.2582 0.0884 0.0090 4.9942 0.1157 0.0100 3.3541 0.1185 数据:
实验五 变步长的龙哥库塔法
实验五变步长的龙格库塔法 一、实验目的 1、了解变步长的龙哥库塔法的特点及具体实现过程。 2、运用vc6.0软件编写龙哥库塔法的算法程序。 3、采用编制的程序，对实际问题进行数值模拟计算。 4、将数值计算结果与精确解进行对比分析，讨论计算误差。 二、实验内容 1、根据变步长的龙哥库塔法法算法的特点，设计程序的流程 本实验采用经典的变步长的龙哥库塔法法格式（四阶变步长的龙哥库塔法法格式） 2、用vc6.0语言编写变步长的龙哥库塔法法算法程序 在编写程序时，充分考虑到程序的交互性和实用性。本程序可以实现计算机主动提示操作步骤和输出结果； 输入待计算的信息后，计算机自动生成图像，并且将得到的结果也显示在图像中，便于直观观察； 3、算法的源程序 #include "stdio.h" #include "math.h" double f(double,double); main() { double y0=1,k1,k2,k3,k4,x0=0,h=0.1,h1=0.05,q,tol=0.1,y00,y01; int i,n=10; k1=f(x0,y0); k2=f(x0+h/2,y0+h/2*k1); k3=f(x0+h/2,y0+h/2*k2); k4=f(x0+h/2,y0+h*k3); y0=y0+h/6*(k1+2*k2+2*k3+k4); y00=y0; printf("the y00 is: %lf\n",y00); /*步长为h时从x0出发求一步得y00*/ k1=f(x0,y0); k2=f(x0+h1/2,y0+h1/2*k1); k3=f(x0+h1/2,y0+h1/2*k2); k4=f(x0+h1/2,y0+h1*k3); y0=y0+h1/6*(k1+2*k2+2*k3+k4); x0=x0+h1; k1=f(x0,y0); k2=f(x0+h1/2,y0+h1/2*k1); k3=f(x0+h1/2,y0+h1/2*k2);
龙格库塔法-原理及程序实现
龙格库塔法一、基本原理：
可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法，即： yi＋１＝yi＋h*( K1+ 2*K2 +2*K3+ K4)/6 K1=f(xi,yi) K2=f(xi+h/2,yi+h*K1/2) K3=f(xi+h/2,yi+h*K2/2) K4=f(xi+h,yi+h*K3) 通常所说的龙格-库塔法就是指四阶——龙格库塔法，我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格-库塔法公式。 （1） 计算公式（1）的局部截断误差是。 龙格-库塔法具有精度高，收敛，稳定（在一定条件下），计算过程中可以改变步长，不需要计算高阶导数等优点，但仍需计算 在一些点上的值，如四阶龙格-库塔法每计算一步需要计算四次 的值，这给实际计算带来一定的复杂性，因此，多用来计算“表头”。 二、小程序 #include
#define f(x,y) (-1*(x)*(y)*(y)) void main(void) { double a,b,x0,y0,k1,k2,k3,k4,h; int n,i; printf("input a,b,x0,y0,n:"); scanf("%lf%lf%lf%lf%d",&a,&b,&x0,&y0,&n); printf("x0\ty0\tk1\tk2\tk3\tk4\n"); for(h=(b-a)/n,i=0;i!=n;i++) { k1=f(x0,y0); k2=f(x0+h/2,y0+k1*h/2); k3=f(x0+h/2,y0+k2*h/2); k4=f(x0+h,y0+h*k3); printf("%lf\t%lf\t",x0,y0); printf("%lf\t%lf\t",k1,k2); printf("%lf\t%lf\n",k3,k4); y0+=h*(k1+2*k2+2*k3+k4)/6; x0+=h; } printf("xn=%lf\tyn=%lf\n",x0,y0); }
龙格库塔方法解题
龙格库塔方法 从常微分方程数值解法的几何意义看，欧拉方法取一点n t 处的斜率() n n Q t f k ,1 =作为平均斜率，因此欧拉方法近似公式为 t k Q Q n n ?+=+11 向后的欧拉方法则采用点1+n t 处的斜率 ()112,++=n n Q t f k 作为平均斜率，即 t k Q Q n n ?+=+21 所以这两种方法也称作矩形法。改进的欧拉方法则取点n t 处和点1+n t 处斜率 1k 和2k 的平均值作为平均斜率，即 ()t k k Q Q n n ?++=+2112 1 因此改进的欧拉方法又称为梯形方法。可以预见，若去多点处斜率的加权平均值作为平均斜率，误差会更小，这就是龙格-库塔方法。 最常用的是四阶龙格库卡近似计算公式，即： ()t k k k k Q Q n n ?++++=+43211226 1 式中 () () 314221312 121,2,2,,tk Q t f k k t Q t f k k t Q t f k Q t f k n n n n n n n n ?+=??? ? ???+=???? ???+==+++ 使用标准四阶龙格——库塔法求解初值问题
()(){}10210≤≤='=x xy y y 计算所用程序如下所示 #include "stdio.h" #include "conio.h" float func(float x,float y) { return(2*x*y); } float runge_kutta(float x0,float xn,float y0,int n) { float x,y,y1,y2,h,xh; float d1,d2,d3,d4; int i; x=x0; y=y0; h=(xn-x0)/n; for(i=1;i<=n;i++) { xh=x+h/2; d1=func(x,y); d2=func(xh,y+h*d1/2.0); d3=func(xh,y+h*d2/2.0); d4=func(xh,y+h*d3); y=y+h*(d1+2*d2+2*d3+d4)/6.0; x=x0+i*h; } return(y); } void main() { float x0,xn,y0,e; int n; printf("\ninput n:\n"); scanf("%d",&n); printf("input x0,xn:\n"); scanf("%f%f",&x0,&xn); printf("input y0:\n"); scanf("%f",&y0); e=runge_kutta(x0,xn,y0,n);
龙格-库塔方法
8.2 龙格-库塔方法 8.2.1 二阶龙格-库塔方法 常微分方程初值问题： 做在点的泰勒展开： 这里。取，就有 (8.11) 截断可得到近似值的计算公式，即欧拉公式： 若取，式（8.11）可写成：
或 （8.12） 截断可得到近似值的计算公式： 或 上式为二阶方法，一般优于一阶的欧拉公式（8.2），但是在计算时，需要计算 在点的值，因此，此法不可取。 龙格-库塔设想用在点和值的线性组合逼近式（8.12）的主体，即用 （8.13） 逼近
得到数值公式： （8.14） 或更一般地写成 对式（8.13）在点泰勒展开得到： 将上式与式（8.12）比较，知当满足 时有最好的逼近效果，此时式（8.13）－式（8.14）。这是4个未知数的3个方程，显然方程组有无数组解。
若取，则有二阶龙格-库塔公式，也称为改进欧拉公式： （8.15） 若取，则得另一种形式的二阶龙格-库塔公式，也称中点公式： （8.16） 从公式建立过程中可看到，二阶龙格-库塔公式的局部截断误差仍为，是二阶精度的计算公式。类似地，可建立高阶的龙格-库塔公式，同时可知四阶龙格-库塔公式的局部 截断误差为，是四阶精度的计算公式。 欧拉法是低精度的方法，适合于方程的解或其导数有间断的情况以及精度要求不高的情况，当解需要高精度时，必须用高阶的龙格-库塔等方法。 四阶龙格-库塔方法应用面较广，具有自动起步和便于改变步长的优点，但计算量比一般方法略大。为了保证方法的收敛性，有时需要步长取得较小，因此，不适于解病态方程。 8.2.2 四阶龙格-库塔公式 下面列出常用的三阶、四阶龙格-库塔计算公式。 三阶龙格-库塔公式
龙格库塔法的编程
龙格库塔法的编程 #include 利用龙格-库塔法求解质点运动常微分方程 一、待解问题 讨论常微分方程的初值问题，边值问题的数值解法，最常用的基本方法就是龙格-库塔法使一些物理方程的计算简便，所得结果的精确提高。 下面是利用龙格-库塔法求解一个质点运动的常微分方程 以初速s km v /80=自地球表面（半径km R 6000=）竖直向上发射一质量为 m 的火箭，如图所示。若不计空气阻力，火箭所受引力 F 大小与它到地心距离的平方成反比，求火箭所能到达的最大高度H 。 火箭为我们分析的对象，对其做受力分析，火箭在任意位置 x 处仅受地球引力F 作用。由题意知，F 的大小与 x 2 成反比，设μ为比例系数，则有： （1） 当火箭处于地面时，即R x =时F = m g ，由式(1)可得2mgR =μ，于 是火箭在任意位置 x 处所受地球引力 F 的大小为 22 x mgR F = 由于火箭作直线运动，火箭的直线运动微分方程式为： 2x F μ =2 222x mgR dt x d m -= 分离变量积分 22 vx gR dx dv -= 二、物理机理 1.龙格-库塔法的基本思想及一般形式 设初值问题],[)(b a x y y ∈=，由微分中值定理可知，必存在],[1+∈n n x x ξ，使 设)(n n x y y =，并记))(,(*ξξy f K =，则 *1)(hK y x y n n +=+ 其中*K 称为)(x y 在][1,+n n x x 上的平均斜率，只要对平均斜率*K 提供一种算法，上式就给出了一种数值解公式，例如，用),(1n n y x f K =代替*K ，就得到欧拉公式，用),(112++=n n y x f K 代替*K ，就得到向后欧拉公式，如果用21,K K 的平均值来代替*K ，则可得到二阶精度的梯形公式。可以设想，如果在],[1+n n x x 上能多预测几个点的斜率值，用它们的加权平均值代替*K ，就有望的到具有较高精度的数值解公式，这就是龙格-库塔（Runge-Kutta ）法的基本思想。 龙格-库塔公式的一般形式： )) (,()()()()(1ξξξy hf x y y h x y x y n n n +='+=+dx dv v dt dx dx dv dt dv dt x d =?==2 22 2 x mgR dx dv mv - =?利用龙格库塔法求解质点运动方程