四阶龙格_库塔算法的C语言实现_毋玉芝

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，该函数接受一个函数指针、初始值、步长和迭代次数作为参数，并返回求解得到的积分值。

题目： 设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&lt;stdlib.h&gt;#include&lt;stdio.h&gt;/*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是一种用于算法开发、数据分析、可视化和数值计算的高级技术计算语言和交互式环境。

作为一个强大的工具，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是初始条件。

《数值分析》课程设计常微分方程初值问题数值解的实现和分析—四阶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 方法。

—科教导刊（电子版）·2019年第24期/8月（下）—259四阶龙格-库塔方法的程序设计与应用罗丽珍吴庆军（玉林师范学院数学与统计学院广西·玉林537000）摘要本文通过介绍四阶龙格-库塔方法，通过预报斜率和泰勒展开式推导出龙格—库塔格式。

了解它的基本思想与算法步骤、MATLAB 语言编写的程序。

列举一些例子，运用四阶龙格-库塔方法的MATLAB 程序在软件中运行，求解出常微分方程的数值解，同时将求解出的数值解与精确解进行比较。

关键词龙格-库塔方法常微分方程数值解中图分类号：TP337文献标识码：A 0引言从17世纪以来国内外数学家对常微分方程的研究取得了很多的成果.欧拉在研究中指出常微分方程存在唯一解和无数解，他用近似值求解微分方程，发现用积分因子求解微积分方程的特殊算法。

拉格朗日建立了一阶微分方程理论，他将参数变法应用到四阶非齐次方程的求解。

我们生活中许多问题的解决都运用到常微分方程，常微分方程的数值解法中经常使用的方法是四阶龙格-库塔方法。

各个领域和工程问题中的原理和演变规律都是用常微分方程来描述的，如在物理方面的电路中电流变化的规律、航天航空方面卫星运转问题、经济方面物品供给以及需求与物价的之间的关系、军事方面研究深水炸弹在水下的运动等。

对这些事物、现象变化规律的描述、认知和分析，需要运用常微分方程来解决。

人们使用常微分方程数值解法的四阶龙格-库塔方法去研究这些问题很实用，而且具有很重要的应用价值。

目前，常微分方程在解决我们生活中的问题很实用，许多问题都运用常微分方程来求解。

中国科学技术大学学者倪兴在常微分方程的研究中写了关于欧拉法、方法等几种方法，他运用常微分计算卫星运动的初轨，把方法运用到卫星轨道改进的例子中；扬州大学学者冯建强和孙诗一研究四阶方法的推导，他写出了如何推导的过程。

在高校数值分析、数值计算方法与实验等教材中，许多作者都出版关于常微分方程初值问题数值解法的教材书，欧拉方法、改进欧拉法和方法等，同时在教材书中写入各种实际问题的例子，运用这些方法去解决常微分方程的初值问题。

某一地区的病菌传染，三种人员的人数的状态方程，即可能受传染的人数x 1,已被传染的病的人数x2，及已治愈的人数x3，并假设他们是因接触而传染。

设：α是x1中单位时间内的传染系数，β是x2中单位时间内被治愈的比例系数，则得以下状态方程：f`(x1)=-α*x1*x2;f`(x2)=α*x1*x2-β*x2;f`(x3)=β*x2;其中f`(x)是加速度本题中，α=0.01，β=0.027初始条件：t=0, x1=620人, x2=10人, x3=70人选步长：h=0.025, 用四阶RK法模拟。

-------------程序清单-------------------#include <stdio.h>#include <conio.h>#define A 0.01#define B 0.027void myFun(double *x,double *f){f[0]=-A*x[0]*x[1];f[1]=A*x[0]*x[1]-B*x[1];f[2]=B*x[1] ;}int rht4(double *rs,double h,double x[],int n,int Num){int j,i,ii;double ps[4]={0,0.5,0.5,1};double bs[4]={1.0/6,1.0/3,1.0/3,1.0/6};double K[3];double y[3];double f[3];for(i=0;i<n;i++){y[i]=x[i];rs[i]=y[i];}for(j=1;j<Num;j++){for(ii=0;ii<4;ii++){ /*龙格-库塔*/for(i=0;i<n;i++){x[i]=rs[(j-1)*n+i]+ps[ii]*f[i]*h; /*特别注意rs[(j-1)*n+i] 而不是x[i] 不是前次的*/}myFun(x,f);/*算每次的加速度*/for(i=0;i<n;i++){y[i]+=f[i]*h*bs[ii];}} /*end 龙格-库塔*/for(i=0;i<n;i++){x[i]=y[i];rs[j*n+i]=y[i];}}}int main(){double t=0,x[3]={620.0,10.0,70.0},h=0.025,*rh;int Num,i,j;double f[3];printf("Input Num: ");scanf("%d",&Num);rh=(double *)malloc(Num*3*sizeof(double));rht4(rh,h,x,3,Num) ;printf("x1 x2 x3-----------------------------------------------------"); for(i=0;i<Num;i++){printf(" ");for(j=0;j<3;j++)printf("%f ",rh[i*3+j]);}getch();return 0;}－－－－－－－－－－end－－－－－－－－－－－－－－－。

实验四 龙格库塔法【实验内容】1. 用标准四阶龙格－库塔方法设计算法2. 编程解微分方程初值问题；【实验方法与步骤】1． 龙格－库塔方法的基本思路设法计算(),f x y 在某些点上的函数值，然后对这些函数值做线性组合，构造近似计算公式；再把近似公式和解的泰勒展开式相比较，使前面的若干项吻合，从而达到较高的精度。

2． 四阶龙格－库塔方法的计算步骤求解'0()(,)()()y x f x y a x b y a y ⎧=≤≤⎨=⎩ 对上述给定的(,)f x y ，用四阶龙格-库塔法求解常微分方程初值问题112341213243(22)6(,)11(,)2211(,)22(,)n n n n n n n n n n h y y k k k k k hf x y k hf x h y k k hf x h y k k hf x h y k +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ 程序如下#include<stdio.h>void Runge_Kutta(float( * f)(float x,float y),float a,float b,float y0,int N){ float x=a,y=y0,K1,K2,K3,K4;float h=(b-a)/N;int i;printf("x[0]=%f\ty[0]=%f\n",x,y);for(i=1;i<=N;i++){ K1=( * f)(x,y);K2=( * f)(x+h/2,y+h*K1/2);K3=( * f)(x+h/2,y+h*K2/2);K4=( * f)(x+h,y+h*K3);y=y+h*(K1+2*K2+2*K3+K4)/6;x=a+i*h;printf("x[%d]=%f\ty[%d]=%f\n",i,x,i,y); }}float f(float x,float y){ return y-2*x/y;}void main(){float a=0,b=1,y0=1;Runge_Kutta(f,a,b,y0,5);}[实验内容]课本Ｐ１８２页例７．３．２【实验结果】。

龙格库塔法一、基本原理：“龙格-库塔(Runge-Kutta)方法”是一种在工程上应用广泛的“高精度单步算法”。

由于此算法精度高，采取措施对误差进行抑制，可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法，即：yi＋１＝yi＋h*( K1+ 2*K2 +2*K3+ K4)/6K1=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<stdio.h>#include<math.h>#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);}运行结果：input a,b,x0,y0,n:0 5 0 2 20x0y0k1k2k3 k40.000000 2.000000-0.000000-0.500000-0.469238-0.8861310.250000 1.882308-0.885771-1.176945-1.129082-1.2800600.500000 1.599896-1.279834-1.295851-1.292250-1.2227280.750000 1.279948-1.228700-1.110102-1.139515-0.9901621.000000 1.000027-1.000054-0.861368-0.895837-0.7528521.2500000.780556-0.761584-0.645858-0.673410-0.5621891.5000000.615459-0.568185-0.481668-0.500993-0.4205371.7500000.492374-0.424257-0.361915-0.374868-0.3178552.0000000.400054-0.320087-0.275466-0.284067-0.2435982.2500000.329940-0.244935-0.212786-0.218538-0.1894822.5000000.275895-0.190295-0.166841-0.170744-0.1495632.7500000.233602-0.150068-0.132704-0.135399-0.1197033.0000000.200020-0.120024-0.106973-0.108868-0.0970483.2500000.172989-0.097256-0.087300-0.088657-0.0796183.5000000.150956-0.079757-0.072054-0.073042-0.0660303.7500000.132790-0.066124-0.060087-0.060818-0.0553054.0000000.117655-0.055371-0.050580-0.051129-0.0467434.2500000.104924-0.046789-0.042945-0.043363-0.0398334.5000000.094123-0.039866-0.036750-0.037072-0.0342024.7500000.084885-0.034226-0.031675-0.031926-0.029571xn=5.000000yn=0.076927。

《计算机数值方法》课程设计报告题目四阶Runge-Kutta方法学生姓名班级计科12学号成绩指导教师延安大学计算机学院2014年9月1日目录一、摘要 (5)二、问题重述 (5)三、方法原理及实现 (5)四、计算公式或算法 (5)五、Matlab程序 (6)六、测试数据及结果 (6)七、结果分析 (10)八、方法改进 (10)九、心得体会 (10)十、参考文献 (10)一、摘要本课程设计主要内容是用四阶Runge-Kutta 方法解决常微分方程组初值问题的数值解法，首先分析题目内容和要求，然后使用Matlab 编写程序计算结果并绘图，最后对计算结果进行分析并得出结论。

二、问题描述在计算机上实现用四阶Runge —Kutta 求一阶常微分方程初值问题()()[]()⎩⎨⎧=∈=1,,,,'y a y b a x y x f x y的数值解，并利用最后绘制的图形直观分析近似解与准确解之间的比较.三、方法原理及实现龙格—库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法.由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂.该算法是构建在数学支持的基础之上的。

龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均,用做下一点的依据，从而构造出了精度更高的数值积分计算方法.如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。

经典的R K -方法是一个四阶的方法,它的计算公式是：112341213243(22)6(,)(,)22(,)22(,)n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ 四、计算公式或算法1． 输入()y x f y n b a ,,,,,0（编写或调用计算()y x f ,的函数文件()y x F ,）,2． b -a h n=（） ()00x y y = b x a x n ==,03．For n i :1=()()()()()101111211111311112411312,,,,i i i i i i i i i x x i hhh K f x y K f x h y h K K f x h y h K K f x h y hK ---------=+-===++=++=++ ()43211226K K K K h y y i i ++++=- End4．输出.,21n y y y ⋅⋅⋅ 五、Matlab 程序x=[a ：h:b ］；y(1）=y1；n=(b-a)/h+1;for i=2：nfk1=f(x （i-1),y （i —1））；fk2=f （x(i-1)+h/2,y （i —1)+fk1＊h/2）；fk3=f(x(i —1）+h/2，y （i-1)+fk2*h/2）;fk4=f(x （i —1）+h,y(i —1)+fk3*h ）;y(i ）=y （i-1)+h ＊(fk1+2*fk2+2*fk3+fk4)/6；endy六、测试数据及结果用调试好的程序解决如下问题:应用经典的四阶Runge —Kutta 方法解初值问题21(),13,(1)2,y y y t t y ⎧'=+≤≤⎪⎨⎪=-⎩ 取0.5h = （1） 步骤一：编写函数具体程序。

例题一四阶龙格-库塔法的具体形式为：1.1程序:/*e.g: y'=t-y,t∈[0,1]/*y(0)=0/*使用经典四阶龙格-库塔算法进行高精度求解/* y(i＋1)＝yi＋( K1+ 2*K2 +2*K3+ K4)/6/* K1=h*f(ti,yi)/* K2=h*f(ti+h/2,yi+K1*h/2)/* K3=h*f(ti+h/2,yi+K2*h/2)/* K4=h*f(ti+h,yi+K3*h)*/#include <stdio.h>#include <math.h>#define N 20float f(float d,float p) //要求解的微分方程的右部的函数e.g: t-y {float derivative;derivative=d-p;return(derivative);}void main(){float t0; //范围上限float t; //范围下限float h; //步长float nn; //计算出的点的个数,即迭代步数int n; //对nn取整float k1,k2,k3,k4;float y[N]; //用于存放计算出的常微分方程数值解float yy; //精确值float error;//误差int i=0,j;//以下为函数的接口printf("input t0:");scanf("%f",&t0);printf("input t:");scanf("%f",&t);printf("input y[0]:");scanf("%f",&y[0]);printf("input h:");scanf("%f",&h);// 以下为核心程序nn=(t-t0)/h;printf("nn=%f\n",nn);n=(int)nn;printf("n=%d\n",n);for(j=0;j<=n;j++){yy=t0-1+exp(-t0); //解析解表达式error=y[i]-yy; //误差计算printf("y[%f]=%f yy=%f error=%f\n",t0,y[i],yy,error);//结果输出k1=h*f(t0,y[i]); //求K1k2=h*f((t0+h/2),(y[i]+k1*h/2)); //求K2k3=h*f((t0+h/2),(y[i]+k2*h/2)); //求K3k4=h*f((t0+h),(y[i]+k3*h)); //求K4y[i+1]=y[i]+((k1+2*k2+2*k3+k4)/6); //求y[i+1]t0+=float(h);i++;}}1.2结果:input t0:0input t:1input y[0]:0input h:0.2nn=5.000000n=5ti yi yy error0 0 0 00.2 0.019736 0.018731 0.001005 0.4 0.074813 0.070320 0.004493 0.6 0.158303 0.148812 0.0094910.8 0.264635 0.249329 0.0153061.0 0.389331 0.367879 0.021452图一例题二(见《高等流体力学》P45页例1.6)给定速度场u=x+t, v=-y-t, 求t=1时过（1,1）点的质点的迹线。

matlab经典的4级4阶runge kutta法-回复什么是Matlab经典的4级4阶Runge-Kutta法（RK4）？Matlab经典的4级4阶Runge-Kutta法是一种数值解法，用于求解常微分方程（ODEs）。

该方法是由德国数学家卡尔·雷特克（Carl Runge）和瓦尔特·库塔（Martin Kutta）在1900年至1901年期间独立发现和发展的。

它是一种四阶精度的方法，具有较高的稳定性和所需步骤数较少的优点，因此被广泛应用于科学和工程领域的数值模拟和计算中。

Runge-Kutta法的基本思想是通过逐步逼近来计算给定的常微分方程。

根据初值条件，我们可以设置一个初始点并迭代地计算出下一个点的近似值。

使用RK4方法，我们可以通过四个步骤来计算下一个点的近似值。

那么，RK4的四个步骤是什么呢？首先，我们设定初值条件。

对于一个一阶常微分方程y' = f(x, y)，我们需要给出初始点(x0, y0)。

然后，我们选择一个适当的步长h，以确定我们在每个步骤中应该前进的距离。

其次，我们计算k1，这是在初始点上斜率的近似值。

我们使用初始点的坐标(x0, y0)和方程f(x, y)来计算k1 = f(x0, y0)。

接下来，我们计算k2，这是在前进到下一个点之前我们在中间点上斜率的近似值。

中间点的坐标为(x0 + h/2, y0 + h*k1/2)。

我们使用方程f(x, y)和中间点的坐标来计算k2 = f(x0 + h/2, y0 + h*k1/2)。

然后，我们计算k3，这是在前进到下一个点之前我们在另一个中间点上斜率的近似值。

另一个中间点的坐标为(x0 + h/2, y0 + h*k2/2)。

我们使用方程f(x, y)和另一个中间点的坐标来计算k3 = f(x0 + h/2, y0 + h*k2/2)。

最后，我们计算k4，这是在下一个点上的斜率的近似值。

下一个点的坐标为(x0 + h, y0 + h*k3)。

用VB编写的一个经典程序，工程应用很广泛四阶龙格库塔法解一阶二元微分方程//dxi/dt=c*(xi-xi^3/3+yi)+K*(X-xi)+c*zi//dyi/dt=(xi-b*yi+a)/c//i=1,2,3//X=sum(xi)/N#include <stdio.h>#include <conio.h>#include <math.h>#include <stdlib.h>#define N 1000 //定义运算步数；#define h 0.01 //定义步长；float a,b,c;//定义全局变量常数a，b，c//定义微分方程：double fx(double x[],double dx,double y[],double dy,double z[],int i,double k,double xavg){ int j;double xi,yi;xi=x[i]+dx;yi=y[i]+dy;return c*(xi-pow(xi,3)/3+yi)+k*(xavg-xi)+c*z[i];}double fy(double x[],double dx,double y[],double dy,int i){double xi,yi;xi=x[i]+dx;yi=y[i]+dy;return (xi-b*yi+a)/c;}void main(){double Kx[3][4],Ky[3][4],x[3]={1,2,3},y[3]={2,3,4},xavg,k=0;//定义x,y的初值；double z[3]={0};int i,j,m,n,S;FILE *fp1,*fp;fp=fopen("sjy.txt","w");fp1=fopen("sjykxy.txt","w");fprintf(fp1,"k\tx1\tx2\tx3\ty1\ty2\ty3\n");if(fp==NULL||fp1==NULL){printf("Failed to open file.\n");getch();return;}printf("Input the value of const a,b,c(seperated by ,eg 0.1,0.2,0.3):");scanf("%f,%f,%f",&a,&b,&c);printf("Input the three values of z(seperated by spacekey,eg 0.1 0.2 0.3):");for (m=0;m<3;++m){scanf("%lf",&z[m]);}printf("Input the value of Steps to get different values of xt,yt(S):");scanf("%d",&S);while(k<=1){fprintf(fp,"k=%.3f\n",k);fprintf(fp,"t\tx1\tx2\tx3\ty1\ty2\ty3\n");fprintf(fp1,"\n%.3f",k);for(j=1;j<N;++j){printf("%.3lf",j*h);fprintf(fp,"%.3lf",j*h);xavg=0;for(i=0;i<3;++i){xavg+=x[i];}xavg=xavg/3.0;//四阶龙格库塔法：for(i=0;i<3;++i){Kx[i][0]=fx(x,0,y,0,z,i,k,xavg);Ky[i][0]=fy(x,0,y,0,i);Kx[i][1]=fx(x,h/2*Kx[i][0],y,h/2*Ky[i][0],z,i,k,xavg);Ky[i][1]=fy(x,h/2*Kx[i][0],y,h/2*Ky[i][0],i);Kx[i][2]=fx(x,h/2*Kx[i][1],y,h/2*Ky[i][1],z,i,k,xavg);Ky[i][2]=fy(x,h/2*Kx[i][1],y,h/2*Ky[i][1],i);Kx[i][3]=fx(x,h*Kx[i][2],y,h*Ky[i][2],z,i,k,xavg);Ky[i][3]=fy(x,h*Kx[i][2],y,h*Ky[i][2],i);x[i]=x[i]+(Kx[i][0]+2*Kx[i][1]+2*Kx[i][2]+Kx[i][3])/6*h;y[i]=y[i]+(Ky[i][0]+2*Ky[i][1]+2*Ky[i][2]+Ky[i][3])/6*h;}for(i=0;i<3;++i){printf("\t%.3lf",x[i]);fprintf(fp,"\t%.3lf",x[i]);}for(i=0;i<3;++i){printf("\t%.3lf",y[i]);fprintf(fp,"\t%.3lf",y[i]);}printf("\n");fprintf(fp,"\n");//取第S步，即时间为S*h的x1，x2，x3，y1，y2，y3随k值的变化；while(j==S){for(n=0;n<3;++n){fprintf(fp1,"\t%.3lf",x[n]);}for(n=0;n<3;++n){fprintf(fp1,"\t%.3lf",y[n]);}break;}continue;}k+=0.1;}fclose(fp1);fclose(fp);printf("\nFinished.\nDatas have been saved to \"sjy.txt,sjykxy.txt\".\n");getch();}。