当前位置：文档视界 › 变步长梯形公式

变步长梯形公式

实验七变步长梯形公式

（一）实验目的

掌握变步长梯形公式的应用。

（二）实验项目内容

1．写出变步长梯形公式步骤和流程图。

2．对算法用C程序或C#、C++实现。

3．调试程序。

（三）主要仪器设备

微机

（四）实验室名称

公共计算机实验室

（五）实验报告撰写

程序，运行结果

流程图：

运用程序：

#include"iostream"

#include"math.h"

using namespace std;

double fx(double x1); //函数声明

void main()

{

double xi,T=0,Sn=0,h=0,a=0,b=1,d=0,R1=0,R2=1; int N=11;

//cout<<"a=,b="<<'\n';

//cin>>a;

//cin>>b;

h=b-a;

T=(fx(a)+fx(b))*h/2;

cout<<"T1="<

for(int i=1;i<=N;i++)

{

h=h/2;//(pow(2.0,i)); //循环一次h减半 double T1=0;

double x=a+h;

while(x

{

T1+=fx(x);

x+=2*h;

}

T=T/2+h*T1;

cout<<"T"<

}

double fx(double x1)

{

double Y;

if(x1!=0)

{

Y=sin(x1)/x1;

}

if(x1==0)

{

Y=1;

}

return Y;

}

云行结果：

复化积分法(复化梯形求积-复化Simpson公式-变步长求积法)MATLAB编程实验报告 (1)

复化积分法（复化梯形求积，复化Simpson 公式，变步长求积法） MATLAB 编程实验报告一、问题描述：编写函数实现复化积分法。二、实验步骤（过程）：（一）复化求积法（1）复化梯形求积：用复化梯形求积公式求解 dx x x ?10sin function [f]=Tn(a,b,n,y) syms t; h=(b-a)/n; f=0; for k=1:n+1 x(k)=a+(k-1)*h z(k)=subs(y,t,x(k)); end for i=2:n f=f+z(i); end q=subs(y,t,a); if y=='sin(t)/t'&&a==0 q=1; end p=subs(y,t,b); T=h/2*(q+p+2*f); T=vpa(T,7) clc,clear; syms t; a=0;b=1; y=sin(t)/t; n=8; Tn(a,b,n,y); （2）复化Simpson 公式：用复化Simpson 公式求解?211dx e x function [f]=simpson(a,b,n,y)

syms t; h=(b-a)/n; f=0;l=0; for k=1:n+1 x(k)=a+(k-1)*h w(k)=0.5*h+x(k) z(k)=subs(y,t,x(k)); end for i=2:n f=f+z(i); end for i=1:n l=l+w(i); end q=subs(y,t,a); if y=='sin(t)/t'&&a==0 q=1; end p=subs(y,t,b); T=h/2*(q+p+2*f); T=vpa(T,7) clc,clear; syms t; a=1;b=2; y=exp(1/t); n=5; simpson(a,b,n,y); （3）变步长求积法：以书本例4.5为例function [f]=TN(a,b,y,R0) syms t; T=[]; f=0; q=subs(y,t,a); if y=='sin(t)/t'&&a==0 q=1; end p=subs(y,t,b); T(1)=(b-a)/2*(q+p); i=2; n=i-1; h=(b-a)/n; z1=a+h/2; z2=subs(y,t,z1);

复化梯形法复化矩形法变步长梯形变步长辛普森

陕西科技大学机械教改班用C++的积分其实积分的思想就是,微分—>求和—>取极限，如果是用纯手工法那就是先对一个函数微分，再求出它的面积，在取极限，因为我们的计算速度和计算量有限，现在有了计算机这个速度很快的机器，我们可以把微分后的每个小的面积加起来，为了满足精度，我们可以加大分区，即使实现不了微分出无限小的极限情况，我们也至少可以用有限次去接近他，下面我分析了四种不同的积分方法，和一个综合通用程序。一．积分的基本思想 1、思路：微分—>求和—>取极限。 2、Newton —Leibniz 公式 ?-=b a a F b F dx x f ) ()()( 其中，)(x F 被积函数)(x f 的原函数。 3、用计算机积分的思路在积分区间内“微分—>求和—>控制精度”。因为计算机求和不可以取极限，也就是不可以无限次的加下去，所以要控制精度。二．现有的理论 1、一阶求积公式---梯形公式 ?=+-=b a T b f a f a b dx x f )]()([2 )( 他只能精确计算被积函数为0、1次多项式时的积分。 2、二阶求积分公式——牛顿、科特斯公式 ?=+++-=b a S b f a b f a f a b dx x f )]()2(4)([6)( 他只能精确计算被积函数为0、1、2、3次多项式时的积分。三．四种实现方法 1.复化矩形法将积分区间[a,b]等分成n 个子区间: ],[],[],[],[],[112322110n n n n x x x x x x x x x x ---、、、则h=(b-a)/n,区间端点值k x =a+kh

龙格库塔积分算法

龙格库塔法龙格库塔法是常用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家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）

变步长的梯形积分方法的应用

CENTRAL SOUTH UNIVERSITY 数值分析实验报告

变步长的梯形积分方法的应用一、问题背景实际问题当中常常需要计算积分，有些数值方法，如微分方程和积分方程的求解，也都和积分计算相关联。依据人们所熟知的微积分基本定理，对于积分 ()dx x f I a ?=b ，只要找到被积分函数()x f 的原函数()x F ，便有下列牛顿-莱布尼茨（Newton-Leibniz ）公式： ()()()a F b F dx x f b -=?a . 但实际使用这种求积分方法往往有困难，因为大量的被积函数，诸如 ()0sin ≠x x x ，2x e -等，其原函数不能用初等函数表达，故不能用上述公式计算。即使能求得原函数的积分，有时计算也十分困难。例如对于被积函数 ()6 11x x f +=，其原函数 ()C x x x x x x x x F ++-+++??? ??-+=1 313ln 3411arctan 61arctan 3122，计算()a F ，()b F 仍然很困难，另外，当()x f 是由测量或数值计算给出的一张数据表时，牛顿-莱布尼茨公式也不能直接运用。因此有必要研究积分的数值计算问题。二、数学模型由于牛顿-科特斯积分公式在8≥n 时不具有稳定性，故不能通过提高阶数的方法来提高求积精度。为了提高精度通常可以把积分区间划分成若干的子区间（通常是等分），再在每个子区间上用低阶求积公式。这种方法称为复合求积法。复合梯形法虽然方法简单，但是却不能估计积分精度，这有时候是很不方便的。要想控制积分精度，可以采用如下的方法，设积分区间已经划分为n 个子区间，这时再把区间划分更细，给出新的积分结果，如果前后两次积分的差比给定的误差容限小的话，则停止细华否则继续增加积分区间。这种方法原理很简单也容易实现，但是实际计算中一般采用的比较少，因为这种方法比较机械效率不是太高，实际上采用比较多的通常是Romberg 方法。三、算法及流程给定义误差容限小量TOL ，对于()dx x f b a ?，有复合梯形公式

数值分析与算法变步长梯形求积法计算定积分

变步长梯形求积法计算定积分 1．原理：变步长求积法的主要思想是利用若干小梯形的面积代替原方程的积分，当精度达不到要求时，可以通过增加点数对已有的区间再次划分，达到所需精度时即可；其中由于新的式子中有原来n点中的部分项，故可以省略一些计算，符合了计算机计算存储的思想。主要公式：T2n=T n/2+(h/2)*Σf(x k+; 2．C++语言实现方式：通过每次的T n值和新增的函数值点计算T2n，再通过判断|T n-T2n|的大小来判断是否达到精度要求。 3．源程序如下： #include"" #include"" double f(double x)//预先输入的待积分函数 { double s; s=log(x*x); return(s); } double ffts(double a,double b,double eps) { int n,k; double fa,fb,h,t1,p,s,x,t; fa=f(a);

fb=f(b); n=1; h=b-a; t1=h*(fa+fb)/2; p=eps+1; while(p>=eps) { s=0; for(k=0;k<=n-1;k++) { x=a+(k+*h; s=s+f(x); } t=t1/2+h*s/2; p=fabs(t1-t); cout<<"步长n为:"<

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 =++

变步长梯形法

#include #include float f(float x) { float s; s=1/(x*x); return(s); } main() { float a,b,c,h,x,T1,T2,S; float T0; printf("变步长梯形法求积分：\n"); printf("需要求解的积分式为f(x)=1/(x*x)\n"); printf("请输入a: "); scanf("%f",&a); printf("请输入b: "); scanf("%f",&b); printf("请输入c: "); scanf("%f",&c); h=b-a; T1=h*(f(a)+f(b))/2; S=0; x=a+h/2; do { S=S+f(x); x=x+h; }while(x=c) { T1=T0; S=0; x=a+h/2; do { S=S+f(x); x=x+h; }while(x

T2=T1/2+h/2*S; printf("步长为%f时\t Tn=%f \t T2n=%f \t 误差变化：%f\n",h,T1,T2,fabs(T1-T2)); h=h/2;T0=T2; } printf("************************************************************** **************\n"); printf("所求的结果为=%f\n",T2); }

变步长梯形公式(C语言)

变步长梯形公式（C语言）程序： // cehngxu.cpp : 定义控制台应用程序的入口点。 // #include"stdafx.h" #include"stdio.h" #include"math.h" double f(double x) //这里自定义了函数，因为出现sin(0)/0的情况，系统无法计算；{ double y; if (x == 0) y = 1; //把x=0作为一种情况，单独拿出来； else y = sin(x) / x; //正常情况下的函数； return(y); } //编写这个自定义函数便于你在作业中的计算，对于不用的题目只需改动一下函数即可计算； void main() //主函数； { double a , b ,h,k,s,t[4998]; //数组的定义中不能出现变量，故对其任意取值，一般取一个很大的数； int n = 1, m = 0; printf("please input a="); scanf_s("%lf", &a); //在VC环境下用scanf输入没有问题，但是在我这编译环境visual studio 2013下需要用scanf_s()输入； printf("please input b="); scanf_s("%lf", &b); t[1] = (b - a)*(f(a) + f(b)) / 2; do { h = (b - a) / n; s = 0; for (k = 0; k < n; k++) { s += f(a + (k + 0.5)*h); } t[2 * n] = t[n] / 2 + h / 2 * s; //梯形的递推公式 n = 2*n; m = m + 1;

龙格库塔法例题

四阶龙格一库塔法通常所说的龙格一库塔法是指四阶而言的．我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格一库塔公式 (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

数值分析变步长求解积分

#include #include # define N 100 double function(double s) { double result; result=sin(s)/s; return result; } int explode2(int k) { int i,s=1; for(i=1;i

} for(i=0;i<2*2*2*N-1;i++) { sum3+=function(a+(i+i/2)*h/8); } d[0]=sum;d[1]=sum1;d[2]=sum2;d[3]=sum3; for(j=1;j<4;j++) T[i][j]=1/2*T[i][j-1]+h/explode2(j)*d[j]; for(i=1;i<4;i++) { for(j=0;j<4-i;j++) T[i][j]=explode4(i)/(explode4(i)-1)*T[i-1][j+1]-1/explode4(i)*T[i-1][j]; } printf("Romberg=%lf",T[3][0]); }

常微分方程组的四阶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 数据:

变步长,辛普森,梯形公式

#include #include #include #include #include #define ESC 27 #define EPS 0.5e-7 #define f(x1) (-2.0/(x1*x1-1)) #define g(x2) (-2.0/(x2*x2-1)) #define s(x0) (-2.0/(x0*x0-1)) char t; void tixing() { double a1=2,b1=3; double T,h1,x1; int n,i; printf("please input n:"); scanf("%d",&n); h1=(b1-a1)/n; x1=a1;T=0; for(i=1;i

实验五变步长的龙哥库塔法

实验五变步长的龙格库塔法一、实验目的 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 #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 #include /*n表示几等分，n+1表示他输出的个数*/ int RungeKutta(double y0,double a,double b,int n,double *x,double *y,int style,double (*function)(double，double)) { double h=(b-a)/n,k1,k2,k3,k4; int i; // x=(double*)malloc((n+1)*sizeof(double)); // y=(double*)malloc((n+1)*sizeof(double)); x[0]=a; y[0]=y0; switch(style) { case 2: for(i=0;i

利用龙格库塔法求解质点运动方程

利用龙格-库塔法求解质点运动常微分方程一、待解问题讨论常微分方程的初值问题，边值问题的数值解法，最常用的基本方法就是龙格-库塔法使一些物理方程的计算简便，所得结果的精确提高。下面是利用龙格-库塔法求解一个质点运动的常微分方程以初速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 - =?