// y' = x + y

double f1(double x,double y)

{

return x + y

}

// y' = 3y/(1 + x)

double f2(double x,double y)

{

return 3*y/(1 + x)

}

// y' = y * y

double ftest(double x,double y)

{

return y*y

}

void solve( double (*func)(double x,double y),

double minX,double maxX,

double y0,

double h,

double result[][7],int* resultNum

)

{

double K1,K2,K3,K4

double Xn_1,Yn_1

int n = 0

result[n][0] = n

result[n][1] = minX

result[n][2] = 0

result[n][3] = 0

result[n][4] = 0

result[n][5] = 0

result[n][6] = y0

for(n = 1 n * h <= maxX n ++ )

{

Xn_1 = result[n-1][1]

Yn_1 = result[n-1][6]

K1 = (*func)(Xn_1 , Yn_1)

K2 = (*func)(Xn_1 + h/2 , Yn_1 + h/2*K1)

K3 = (*func)(Xn_1 + h/2 , Yn_1 + h/2*K2)

K4 = (*func)(Xn_1 + h , Yn_1 + h*K3)

result[n][0] = n

result[n][1] = minX + n*h

result[n][2] = K1

result[n][3] = K2

result[n][4] = K3

result[n][5] = K4

result[n][6] = Yn_1 + h*(K1 + 2*K2 + 2*K3 + K4)/6

}

*resultNum = n

}

void print(double result[][7],int resultNum)

{

int i

double x

printf("%5s%15s%15s%15s%15s%15s%15s\n","n","Xn","K1","K2","K3","k4","Y")

for(i = 0 i <95i ++ )

printf("-")

printf("\n")

for(i = 0 i <resultNum i ++)

{

printf("%5d%15f%15f%15f%15f%15f%15f\n",

(int)result[i][0],

result[i][1],

result[i][2],

result[i][3],

result[i][4],

result[i][5],

result[i][6]

)

}

for(i = 0 i <95i ++ )

printf("-")

printf("\n\n")

}

int main(int argc, char *argv[])

{

double minX, maxX

double y0

double h

double result[10000][7]

intresultNum

/****************************************

* y'=x+y

* y(0)=1 (0＜x＜1)

****************************************/

printf("y'=x+y y(0)=1 (0＜x＜1)\n")

minX = 0.0,maxX = 1

y0 = 1

h = 0.1

solve( ftest,

minX,maxX,

y0,

h,

result,&resultNum

)

print(result,resultNum)

/****************************************

* y'=3y/(1+x)

* y(0)=1 (0＜x＜1)

****************************************/

printf("y'=3y/(1+x) y(0)=1 (0＜x＜1)\n")

minX = 0.0,maxX = 1

y0 = 1

h = 0.1

solve( ftest,

minX,maxX,

y0,

h,

result,&resultNum

)

print(result,resultNum)

/****************************************

* y'=y*y

* y(0)=1 (0＜x＜0.5)

****************************************/

printf("y'=y*y y(0)=1 (0＜x＜0.5)\n")

minX = 0.0,maxX = 0.5

y0 = 1

h = 0.1

solve( ftest,

minX,maxX,

y0,

h,

result,&resultNum

)

print(result,resultNum)

/****************************************/

return 0

}

/*

运行结果：

y'=x+y y(0)=1 (0＜x＜1)

n Xn K1 K2 K3 k4 Y

-----------------------------------------------------------------------------------------------

0 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000

1 0.100000 1.000000 1.102500 1.113289 1.235052 1.111110

2 0.200000 1.234567 1.375551 1.392136 1.563310 1.249998

3 0.300000 1.562495 1.763910 1.790762 2.042253 1.428566

4 0.400000 2.040801 2.342756 2.389201 2.780510 1.666653

5 0.500000 2.777733 3.259974 3.347626 4.005666 1.999963

6 0.600000 3.999853 4.839806 5.026356 6.263001 2.499883

7 0.700000 6.249414 7.909333 8.383049 11.143498 3.332844

8 0.800000 11.107850 15.118384 16.717986 25.046449 4.996628

9 0.900000 24.966293 38.999310 48.255163 96.474521 9.929124

10 1.000000 98.587505 220.775003 439.6517502904.595478 81.996399

-----------------------------------------------------------------------------------------------

y'=3y/(1+x) y(0)=1 (0＜x＜1)

n Xn K1 K2 K3 k4 Y

-----------------------------------------------------------------------------------------------

0 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000

1 0.100000 1.000000 1.102500 1.113289 1.235052 1.111110

2 0.200000 1.234567 1.375551 1.392136 1.563310 1.249998

3 0.300000 1.562495 1.763910 1.790762 2.042253 1.428566

4 0.400000 2.040801 2.342756 2.389201 2.780510 1.666653

5 0.500000 2.777733 3.259974 3.347626 4.005666 1.999963

6 0.600000 3.999853 4.839806 5.026356 6.263001 2.499883

7 0.700000 6.249414 7.909333 8.383049 11.143498 3.332844

8 0.800000 11.107850 15.118384 16.717986 25.046449 4.996628

9 0.900000 24.966293 38.999310 48.255163 96.474521 9.929124

10 1.000000 98.587505 220.775003 439.6517502904.595478 81.996399

-----------------------------------------------------------------------------------------------

y'=y*y y(0)=1 (0＜x＜0.5)

n Xn K1 K2 K3 k4 Y

-----------------------------------------------------------------------------------------------

0 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000

1 0.100000 1.000000 1.102500 1.113289 1.235052 1.111110

2 0.200000 1.234567 1.375551 1.392136 1.563310 1.249998

3 0.300000 1.562495 1.763910 1.790762 2.042253 1.428566

4 0.400000 2.040801 2.342756 2.389201 2.780510 1.666653

5 0.500000 2.777733 3.259974 3.347626 4.005666 1.999963

-----------------------------------------------------------------------------------------------

请按任意键继续. . .

*/

float x=0,y=1//初始值

int i=1

while(i<100)

{

    float k = y-(2*x)/(3*y)//求斜率，也就是y'

    y+=k*dx

    x+=dx

    printf("x=%f,y=%f\n",x,y)//输出

    i++

}