常微分方程的数值解法

#include <stdio.h>
#include <math.h>

typedef double (*Fun)(double x, double y);

double RK4(Fun f, double x0, double y0, double xn, double h)
{
double k1 = 0.0, k2 = 0.0, k3 = 0.0, k4 = 0.0;
double x = 0.0, y = y0, ns;
int nsteps, i;

// Calculate nsteps
ns = (xn - x0)/h;
ns = ceil(ns);
nsteps = (int) ns; // number of steps
h = (xn - x0)/ns;

for(i = 0; i < nsteps; i++)
{
x = x0 + i * h;
k1 = h * f(x, y);
k2 = h * f(x + h / 2, y + k1 / 2);
k3 = h * f(x + h / 2, y + k2 / 2);
k4 = h * f(x + h, y + k3);
y += k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6;
}
return y;
}

/*y'=-y-y^2sin(alpha)
*assume sin(alpha) = 0.5
*/
double f(double x, double y)
{
return -y - y * y * 0.5;
}

int main()
{
//[1,2]
double ret = RK4(f,