四阶龙格_库塔算法的C语言实现_毋玉芝
c语言中用runge-kutta 4阶算法求积分
Runge-Kutta 4阶算法是一种用于求解常微分方程数值解的方法。
以下是一个使用C语言实现Runge-Kutta 4阶算法的示例代码，用于求解函数f(x) = x^2在区间[0, 1]上的定积分：c复制代码#include<stdio.h>#include<math.h>double f(double x) {return x * x;}double rungeKutta4(double (*f)(double), double x0, double h, int n) {double k1, k2, k3, k4;double x = x0;for (int i = 0; i < n; i++) {k1 = h * f(x);k2 = h * f(x + 0.5 * k1);k3 = h * f(x + 0.5 * k2);k4 = h * f(x + k3);x += (k1 + 2 * k2 + 2 * k3 + k4) / 6;}return x;}int main() {double x0 = 0.0; // 初始值double h = 0.01; // 步长int n = 1000; // 迭代次数double integral = rungeKutta4(f, x0, h, n);printf("The integral of %f from %f to %f is %f\n", f(x0), x0, x0 + h * n,integral);return0;}在上面的代码中，我们定义了一个函数f(x) = x^2，然后使用Runge-Kutta 4阶算法求解该函数在区间[0, 1]上的定积分。
具体来说，我们首先定义了Runge-Kutta 4阶算法的实现函数rungeKutta4，该函数接受一个函数指针、初始值、步长和迭代次数作为参数，并返回求解得到的积分值。
四阶龙格库塔求解边界层问题(C语言)
题目： 设52/ 1.110/m s μρ-=⨯的气体以10/v m s ∞=的速度以零攻角定常扰流长度为L ＝1m 的大平板，用数值解讨论边界层内的流动规律。
如图所示，沿板面方向取x 坐标，板的法线方向取y 坐标。
在流体力学中，介绍了存在相似性解的二维平面不可压缩流体的边界层方程：C.E ：0=∂∂+∂∂yx u υM.E : 221y udx dp y u x u u ∂∂+-=∂∂+∂∂νρυ0p y∂=∂ 边界条件为：y = 0；u = v = 0y =∞；u =v ∞ ( u = u (x) 为x 点处壁面势流流速 )在外边界上221122e e p v p v c ρρ∞+=+=所以 00e dpdx+=对于平板扰流，则平板扰流的边界层方程可简化为 C.E ：0=∂∂+∂∂yx u υ M.E : 22u u uu x y yυν∂∂∂+=∂∂∂ 其定解的边界条件为y = 0；u = v = 0y =∞；u =v ∞为了求解边界层方程，我们可以利用“相似性解”的方法，对其进行转化，由C.E ，引进流函数),(y x ψ，令xy u ∂∂-=∂∂=ψυψ, 由M.E 式的偏导阶次，可望减少自变量数目令()yg x η==()()f f v g x η∞== 其中2x x ηη∂=-∂1()y g x η∂=∂ 由()v g x f ψ∞=，得()'v g x fu v f y yψ∞∞∂∂===∂∂ ()()(')2v g x f v g x v f f x x xψη∞∞∂∂=-==-∂∂ 所以，(')"2v u v f f x x x η∞∞∂∂==-∂∂(')"()v u v f f y y g x ∞∞∂∂==∂∂222(")"'()()v v u f f y y g x g x ∞∞∂∂==∂∂ 将其都代入M.E ，化简可得"'"0f ff += 这就使原方程变化为："'"0f ff +=此时的边界条件为：η = 0；f(0) = 0 , f’(0) = 0η = ∞；f’(∞) = 1那么，如何利用编辑程序的方法求解这个变化后的边界层微分方程呢？一. 解方程的基本思路为了简化运算，此时边界层微分方程化简成："'"0f ff +=，边界条件不变。
龙格库塔方法
龙格库塔方法#include<stdlib.h>#include<stdio.h>/*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#include<stdlib.h>#include<stdio.h>/*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<n;i++){x[i+1]=x[i]+h;k1=function(x[i],y[i]);k2=function(x[i]+h/2,y[i]+h*k1/2);y[i+1]=y[i]+h*k2;}break;case 3:for(i=0;i<n;i++){x[i+1]=x[i]+h;k1=function(x[i],y[i]);k2=function(x[i]+h/2,y[i]+h*k1/2);k3=function(x[i]+h,y[i]-h*k1+2*h*k2);y[i+1]=y[i]+h*(k1+4*k2+k3)/6;}break;case 4:for(i=0;i<n;i++){x[i+1]=x[i]+h;k1=function(x[i],y[i]);k2=function(x[i]+h/2,y[i]+h*k1/2);k3=function(x[i]+h/2,y[i]+h*k2/2);k4=function(x[i]+h,y[i]+h*k3);y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;}break;default:return 0;}return 1;}double function(double x,double y){return y-2*x/y;}//例子求y'=y-2*x/y(0<x<1);y0=1;/*int main(){double x[6],y[6];printf("用二阶龙格-库塔方法\n");RungeKutta(1,0,1,5,x,y,2,function);for(int i=0;i<6;i++)printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]); printf("用三阶龙格-库塔方法\n");RungeKutta(1,0,1,5,x,y,3,function);for(i=0;i<6;i++)printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]); printf("用四阶龙格-库塔方法\n");RungeKutta(1,0,1,5,x,y,4,function);for(i=0;i<6;i++)printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]); return 1;}。
四阶龙格—库塔法的原理及其应用
《四阶龙格—库塔法的原理及其应用》
龙格—库塔法(又称龙格库塔法)是由一系列有限的、独立的可能解组成的无穷序列,这些解中每个都与原来的数列相差一个常数。
它是20世纪30年代由匈牙利著名数学家龙格和库塔提出的,故得此名。
1.它的基本思想是:在n 阶方阵M 上定义一个函数,使得当n 趋于无穷时,它在m 中所表示的数值为M 的某种特征值,从而构造出一族具有某种特性的可计算函数f (x)= Mx+ C (其中C 为任意正整数)。
例如,若f (x)=(a-1) x+ C,则称之为(a-1) x 的龙格—库塔法。
2.它的应用很广泛,可以求解各类问题,且能将大量的未知数变换成少数几个已知数,因此它是近似计算的一种重要工具。
3.
它的优点主要有:(1)可以将多项式或不等式化成比较简单的形式;(2)对于同一问题可以用不同的方法来解决,并取得同样的结果;(3)适合处理高次多项式或者不等式,尤其适合处理多元函数的二次型。
matlab经典的4级4阶runge kutta法
MATLAB是一种用于算法开发、数据分析、可视化和数值计算的高级技术计算语言和交互式环境。
作为一个强大的工具，MATLAB提供了许多数值计算方法，其中4级4阶Runge-Kutta方法就是其中之一。
1. Runge-Kutta方法简介Runge-Kutta方法是求解常微分方程（ODE）的数值方法之一。
在MATLAB中，用户可以使用内置的ode45函数来调用4级4阶Runge-Kutta方法。
具体来说，4级4阶Runge-Kutta方法是一种单步迭代方法，通过在每个步骤中计算斜率来逐步逼近解析解。
它的优点是数值稳定性好，适用于多种类型的微分方程。
2. Runge-Kutta方法的公式4级4阶Runge-Kutta方法的一般形式如下：$$k_1 = hf(t_n, y_n)$$$$k_2 = hf(t_n + \frac{1}{2}h, y_n + \frac{1}{2}k_1)$$$$k_3 = hf(t_n + \frac{1}{2}h, y_n + \frac{1}{2}k_2)$$$$k_4 = hf(t_n + h, y_n + k_3)$$$$y_{n+1} = y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)$$其中，$t_n$是当前的独立变量值，$y_n$是当前的解向量，h是步长，$f(t_n, y_n)$是给定点$(t_n, y_n)$处的斜率。
通过不断迭代上述公式，可以逐步求解微分方程的数值解。
3. MATLAB中的4级4阶Runge-Kutta方法的应用在MATLAB中，用户可以使用ode45函数调用4级4阶Runge-Kutta方法来求解常微分方程。
使用ode45函数的基本语法如下：```matlab[t, y] = ode45(odefun, tspan, y0)```其中，odefun是用户定义的ODE函数句柄，tspan指定了求解的时间范围，y0是初始条件。
常微分方程初值问题数值解的实现和分析—四阶Rungekutta方法与预估校正算法毕业论文
《数值分析》课程设计常微分方程初值问题数值解的实现和分析—四阶Runge-kutta方法及预估-校正算法常微分方程初值问题数值解的实现和分析—四阶Runge-kutta方法及预估-校正算法摘要求解常微分方程的初值问题，Euler方法，改进的Euler方法及梯形方法精度比较低，所以本文构造高精度单步的四级Runge-kutta方法及高精度的多步预估—校正算法及其Matlab编程来实现对常微分方程初值问题的求解，使在求解常微分方程时，对以前积分方法的收敛速度及精度都有了很高的提高。
关键词：Runge-kutta方法，Adams方法，预估—校正算法，Matlab目录1.前言 (1)2. 几个简单的数值积分法 (2)2.1Runge-kutta方法 (2)2.1.1 Runge-kutta方法的应用 (5)2.2预估—校正算法 (7)2.2.1 Adams数值积分方法简介及预估—校正算法 (7)2.2.2 预估—校正算法的应用 (12)3. 结果分析 (16)总结 (17)参考文献 (18)英文原文和中文翻译 (19)1英文原文 (19)2中文翻译 (20)1.前言常微分方程的初值问题是微分方程定解问题的一个典型代表，以下面的例子介绍常微分方程初值问题数值解的基本思想和原理。
例1.1 一重量垂直作用于弹簧所引起的震荡，当运动阻力与速度的平方成正比时，可借助如下二阶常微分方程描述若令和，则上述二阶常微分方程可化成等价的一阶常微分方程组类似于例1.1，对于m阶常微分方程其中。
若定义可得如下等价的一阶常微分方程组我们知道多数常微分方程主要靠数值解法。
所谓数值解法，就是寻求解在一系列离散节点上的近似值。
相邻两个节点之间的间距称为步长[1]。
2. 几个简单的数值积分法2.1 Runge-kutta方法Runge在1985年提出了一种基于Euler折线法的新的数值方法，此后这种新的数值方法又经过其同胞K.Heun和Kutta的努力[2]，发展完善成为后世所称的Runge-kutta 方法。
Runge-Kutta 四阶
7 1.400000 1.516896
8 1.600000 1.245982
9 1.800000 1.062639
10 2.000000 0.999959
请选择函数
1.Y'(x)=y*sin(PIE*x) Y(0)=1
2.Y'(x)=x+y Y(0)=1.442
3.结束
1
请输入区间[a,b]的端点值
a=0
b=0
输入区间[a,b]端点值有误，请重新输入
a=0
b=2
请输入区间剖分点数n
n=10
n Xn Yn
1 0.200000 1.062680
2 0.400000 1.246016
3 0.600000 1.516916
4 0.800000 1.778619
5 1.000000 1.890103
goto main;
}
else;
}
printf("*****谢谢使用*****\n");
}
请ห้องสมุดไป่ตู้择函数
1.Y'(x)=y*sin(PIE*x) Y(0)=1
2.Y'(x)=x+y Y(0)=1.442
3.结束
4
选择错误，请重新选择
请选择函数
1.Y'(x)=y*sin(PIE*x) Y(0)=1
2.Y'(x)=x+y Y(0)=1.442
2.Y'(x)=x+y Y(0)=1.442
3.结束
3
*****谢谢使用*****
请按任意键继续. . .
for(i=0;i<MAX;i++){
四阶龙格-库塔方法的程序设计与应用
—科教导刊（电子版）·2019年第24期/8月（下）—259四阶龙格-库塔方法的程序设计与应用罗丽珍吴庆军（玉林师范学院数学与统计学院广西·玉林537000）摘要本文通过介绍四阶龙格-库塔方法，通过预报斜率和泰勒展开式推导出龙格—库塔格式。
了解它的基本思想与算法步骤、MATLAB 语言编写的程序。
列举一些例子，运用四阶龙格-库塔方法的MATLAB 程序在软件中运行，求解出常微分方程的数值解，同时将求解出的数值解与精确解进行比较。
关键词龙格-库塔方法常微分方程数值解中图分类号：TP337文献标识码：A 0引言从17世纪以来国内外数学家对常微分方程的研究取得了很多的成果.欧拉在研究中指出常微分方程存在唯一解和无数解，他用近似值求解微分方程，发现用积分因子求解微积分方程的特殊算法。
拉格朗日建立了一阶微分方程理论，他将参数变法应用到四阶非齐次方程的求解。
我们生活中许多问题的解决都运用到常微分方程，常微分方程的数值解法中经常使用的方法是四阶龙格-库塔方法。
各个领域和工程问题中的原理和演变规律都是用常微分方程来描述的，如在物理方面的电路中电流变化的规律、航天航空方面卫星运转问题、经济方面物品供给以及需求与物价的之间的关系、军事方面研究深水炸弹在水下的运动等。
对这些事物、现象变化规律的描述、认知和分析，需要运用常微分方程来解决。
人们使用常微分方程数值解法的四阶龙格-库塔方法去研究这些问题很实用，而且具有很重要的应用价值。
目前，常微分方程在解决我们生活中的问题很实用，许多问题都运用常微分方程来求解。
中国科学技术大学学者倪兴在常微分方程的研究中写了关于欧拉法、方法等几种方法，他运用常微分计算卫星运动的初轨，把方法运用到卫星轨道改进的例子中；扬州大学学者冯建强和孙诗一研究四阶方法的推导，他写出了如何推导的过程。
在高校数值分析、数值计算方法与实验等教材中，许多作者都出版关于常微分方程初值问题数值解法的教材书，欧拉方法、改进欧拉法和方法等，同时在教材书中写入各种实际问题的例子，运用这些方法去解决常微分方程的初值问题。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
四阶龙格 库塔算法的C 语言实现
毋玉芝
(焦作财会学校)
摘 要 本文叙述了四阶龙格 库塔算法的C 语言实现过程、数据存储及其结果的曲线显示,并以具体实例说明了这一过程。
关键词 龙格 库塔算法 数据存储 曲线显示
在科学技术中常常需要求解常微分方程的定解问题,这就需要一种合适的数值解法求出常微分方程的解。
在诸多数值算法中龙格 库塔算法具有较高的精确度,是一种优先选取的算法。
1.四阶龙格 库塔算法简述
龙格 库塔方法实际上是间接地使用台劳级数法的一种技术。
龙格 库塔算法的数学描述如下:
y n+1=y n +h (K 1+2K 2+2K 3+K 4)/6;
K 1=f (x n ,y n );
K 2=f (x n +h/2,y n +K 1 h/2);
K 3=f (x n +h/2,y n +K 2 h/2);
K 4=f (x n +h,y n +K 3 h);
其中: h 表示计算过程中选取的步长;K 1表示x n 点处的斜率;
K 2表示利用K 1求得的(x n +h/2)点处的斜率;
K 3表示利用K 2求得的(x n +h/2)点处的斜率;
K 4表示利用K 3求得的(x n +h)点处的斜率;
2.C 语言的实现过程
2 1 算法实现
龙格 库塔算法关键是选择步长h ,必须根据题目的要求选出合适的步长,这对龙格 库塔算法结果的精确度及其平滑性尤为重要。
在选择了恰当的步长后,利用上述迭代表达式,并根据题目要求的迭代次数,或求解的精度,利用C 语言加以实现。
2 2 数据存储
由于此计算结果数据庞大,程序运行后数据不可能一屏显示,鉴于此首先利用fopen ()函数创建并打开一文本文件,利用fprintf ()函数将数据存储到数据文件,可将此文件打印输出,以检验结果的正确性。
实现过程如下:
if(fp==NU LL)
{
printf( Can t open this file\n );
exit(0);
}
for(i=0;i<=j;i++)
{t=ts+i *h;
2001年3月第1期 焦作大学学报JO UR NAL OF JIAO ZUO U NI VERSIT Y 1Mar.2001
56 焦 作 大 学 学 报 2001年3月
fprintf(fp, t=%3.2f\n ,t);
fprintf(fp, x=%8.5f ,x[i]);
fprintf(fp, y=%8.5f ,y[i]);
fprintf(fp, z=%8.5f \n ,z[i]);
}
fclose(fp);
2 3 曲线显示
根据上述所取得的结果,并根据坐标的适当变换即可以描绘出结果的平滑曲线。
根据此曲线可以进一步了解结果的精确程度。
我们可以利用一子函数graph()实现此功能。
首先在graph()函数中要调用C语言的图形函数库,其次利用C语言中的图形函数:moveto ()、outtextxy()、lineto()、setbkcolor()等,并将所取得的结果进行适当的坐标变换即可实现曲线的显示。
最后必须调用closeg raph()函数关闭图形库。
格式如下:
void graph()
{
int gdrive,gmode=VGAH I
g drive=VGA;
initgraph(&g drive,&gmode, d:\tc20 );
(曲线描绘)
closeg raph();
3.实例分析
本实例是解决在 网络理论 中微分方程组的求解及曲线描绘问题。
所求解的微分方程组如下所示:
x 1(t)=-x1(t)+6x2(t)+2x3(t);
x 2(t)=-x2(t)+3x3(t)+2sin(t);
x 3(t)=-x3(t)+t2 e-t+2cos(t);
条件:0 0 t 5 0 迭代次数:100
利用上述过程,根据循环迭代的方法取步长为h=(5.0-0.0)/100,数据处理编程实现如下:
float k1,k2,k3,k4,m1,m2,m3,m4,n1,n2,n3,n4;
float t,ts,tf,h;
float x[101],y[101],z[101];
x[0]=1.0;y[0]=-1.0;z[0]=0.0;ts=0.0;tf=5.0;
h=(tf-ts)/j;
for(i=1;i<=j;i++)/*数据处理*/
{
t=ts+(i-1)*h
k1=-x[i-1]+2*y[i-1]+6*z[i-1];
m1=-y[i-1]+3*z[i-1]+2*sin(t);
n1=-z[i-1]+pow(t,2)*ex p(-t)+cos(t);
k2=-(x[i-1]+k1*(h/2.0))+2*(y[i-1]+m1*(h/2.0))+6*(z[i-1]+n1*(h/2.0));
m2=-(y[i-1]+m1*(h/2.0))+3*(z[i-1]+n1*(h/2.0))+2*sin(t+h/2.0);
n2=-(z[i-1]+n1*(h/2.0))+pow(t+h2.0,2)*exp(-(t+h/2.0))+cos(t+h/2.0);
(下转第69页)
由前面分析可等效为图5(b)、(c),其中r =r (1+Q 2L ),L =L ,R P =R S r R L 。
在谐振频率 0=1L C
下,容易得出电路有载品质因数 Q =R P
L /C
(6)
进一步分析可知道,谐振回路品质因数愈高,回路谐振特性就愈好。
4.结 论
通过以上分析,可以很清楚地看到品质因数有元件品质因数和回路品质因数,其中回路外接负载时又叫有载品质因数。
品质因数实质上是无功功率和有功功率的比值。
元件品质因数是基础,它决定于本身材料和结构,并且随角频率 变化。
回路品质因数是在谐振频率 0处定义的,对于高Q L 电
感和理想电容组成的谐振回路, 0=1LC
,不同的回路有不同的品质因数。
本文中回路品质因数的分析,以高Q L 电感和理想电容为前提,这也符合实际情况。
参考文献
1.卢淦主编.高频电子电路.北京:中国铁道出版社,1983
2.谢沅清主编.现代电子电路与技术.北京:中央电大出版社,1996
(上接第56页)
k3=-(x [i-1]+k2*(h/2.0))+2*(y[i-1]+m2*(h/2.0))+6*(z[i-1]+n2*(h/2.0));
m3=-(y[i-1]+m2*(h/2.0))+3*(z[i-1]+n2*(h/2.0))+2*sin(t+h/2.0);
n3=-(z[i-1]+n2*(h/2.0))+pow (t+h/2.0,2)*ex p(-(t+h/2.0))+cos(t+h/2.0);
k4=-(x [i-1]+k3*(h/2.0))+2*(y[i-1]+m3*(h/2.0))+6*(z[i-1]+n3*(h/2.0));
m4=-(y[i-1]+m3*(h/2.0))+3*(z[i-1]+n3*(h/2.0))+2*sin(t+h/2.0);
n4=-(z[i-1]+n3*(h/2.0))+pow (t+h/2.0,2)*ex p(-(t+h/2.0))+cos(t+h/2.0);
x[i]=x[i-1]+h *(k1+2*k2+2*k3+k4)/6.0;
y[i]=y[i-1]+h *(m1+2*m2+2*m3+m 4)/6.0;
z[i]=z[i-1]+h *(n1+2*n2+2*n3+n4)
/6.0;
}
其结果存放于data txt 文件中,曲线显示如图1
所示。
4.结束语
利用C 语言可方便地实现四阶龙格 库塔算法,
这对于科学技术中所遇到的常微分方程的求解及曲线
显示问题具有重要意义。
参考文献
1.李庆扬,王能超,易大义著.数值分析.武汉:华中理工大学出版社,1995
2.李仁发著.C 语言程序设计基础.北京:中国矿业大学出版社,199569第1期 靳孝峰等:L C 谐振回路的品质因数。