对于方程y'=f(x,y),初始条件:y(x0)=y0,
#include<stdio.h>FILE *fp=fopen("ex的值.dat","w");double func(double x){return x;}void tworder(double x_0, double h, double t_0, double t_n){double n=(t_n-t_0)/h;double t=t_0, x=x_0;double k1, k2, k3, k4;int i;printf("%f\t%f\n", t, x);fprintf(fp,"%f\t%f\n", t, x);for(i=1;i<=n;i++){k1=func(x);k2=func(x+k1*h/2);k3=func(x+k2*h/2);k4=func(x+k3*h);x=x+(k1+2*k2+2*k3+k4)*h/6;t=t+h;printf("%f\t%f\n", t, x);fprintf(fp,"%f\t%f\n", t, x);}}int main(){ double h=0.01;double x_0=1;double t_0=0, t_n=10;tworder(x_0, h, t_0, t_n);fclose(fp);return 0;}
运行后得到x(t=10)=e^10=22026.465777
用龙格-库塔法计算的结果并不是绝对精确。但是，只要时间 t取值不要太大，误差的范围仍然是可以接受的。
实践应用：姿态四元数微分方程
% 四元数微分方程的4阶龙格库塔法% q0:4*1,rotation vector from body-frame to world-frame% gyro：陀螺仪数据% T：更新周期function [ q ] = Quaternion_RungeKutta4( q0,gyro,T)q0=Norm_Quaternion(q0); %归一化K1= Quaternion_Diff( gyro,q0);q1=Norm_Quaternion(q0+T/2*K1);K2 = Quaternion_Diff(gyro,q1);q2=Norm_Quaternion(q0+T/2*K2);K3 = Quaternion_Diff(gyro,q2);q3=Norm_Quaternion(q0+T*K3);K4 = Quaternion_Diff(gyro,q3);q = q0 + T/6*(K1+2*K2+2*K3+K4);q = Norm_Quaternion(q);end% 函数功能：四元数微分方程% 输 出：四元数的一阶导数% 备 注：连续域function [ q_diff ] = Quaternion_Diff( gyro,q)A = [ 0, -gyro(1)/2, -gyro(2)/2, -gyro(3)/2;gyro(1)/2,0, gyro(3)/2, -gyro(2)/2;gyro(2)/2, -gyro(3)/2,0, gyro(1)/2;gyro(3)/2, gyro(2)/2, -gyro(1)/2, 0];q_diff = A*q;end
