分割数を求めるプログラム

[maimai]

クッタ法(3次)にて分割数を求めるプログラムを作りたいのですが、うまく作ることができません。
区間[0,1]をｎ等分して求めるx=1.0の時の解析解と数値解の誤差がε=10^-3以下になる分割数を求めたいです。

以下のようにつくってはみたものの、それっぽい感じにはなったのですが、ちゃんとε=10^-3以下にはなっていません。
分かる方どうかお願い致します。
(使用ソフト:EasyIDEC)

コード:

#include<stdio.h>
#include<math.h>
#define EPS 1.0e-5


double k_func(double ,double ,double);
double func(double ,double);
double func_y(double);

int main(void)
{
	double a=0.0, b=1.0;
	double x=1.0,y,h,k,z,v,s,j;
	int n=10,i;
	printf(" i  x         y        z\n");
	for(n=2.0;n<256;n*=2){
	h=(b-a)/n;
	x=0.0;y=1.0;i=0;

	while(i<n){
		z=func_y(x);
		k=k_func(x,y,h);
		x+=h;y+=h*k;i++;
	
	}
printf("%2.d, %f, %f, %f\n",i,x,y,z);

s=z-y;
		v=fabs(s);
	printf("v=%f\n",v);
	if (v>=EPS) {
	printf("%f,n=%d\n",v,n);
}	
	}return 0;
}

double k_func(double x,double y,double h)
{
	double k1,k2,k3,k;
	k1=func(x,y);
	k2=func(x + h/2.0,y + k1*h/2.0);
	k3=func(x + h,y + (2.0*k2-k1)*h);
	k=(k1 + 4.0*k2 + k3)/6.0;
	return k;
	
}
double func(double x ,double y)
{
	return (3*y)/(x+1);//１階常微分方程式
	
}
double func_y(double x)
{
	return pow(x+1, 3.0);//特殊解
}

[かずま]

ループの中で z を計算しても、直前の値を上書きするから無意味です。
ループ完了後の x の値で z を求めましょう。

コード:

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

#define EPS 1.0e-5

double k_func(double, double, double);
double func(double, double);
double func_y(double);

int main(void)
{
	double a = 0.0, b = 1.0;
	double x = 1.0, y, h, k, z, v, s, j;
	int n, i;

	printf(" i  x         y        z\n");
	for (n = 2; n < 256; n *= 2) {
		h = (b - a) / n;
		x = 0.0;
		y = 1.0;
		i = 0;
		while (i < n) {
			k = k_func(x, y, h);
			x += h;
			y += h * k;
			i++;
		}
		z = func_y(x);  // ★ z = func_y(b) と同じ
		printf("%2.d, %f, %f, %f\n", i, x, y, z);

		s = z - y;
		v = fabs(s);
		printf("v=%f\n", v);
		if (v >= EPS) {
			printf("%f,n=%d\n", v, n);
		}
	}
	return 0;
}

double k_func(double x, double y, double h)
{
	double k1, k2, k3, k;
	k1 = func(x, y);
	k2 = func(x + h / 2.0, y + k1 * h / 2.0);
	k3 = func(x + h, y + (2.0 * k2 - k1) * h);
	k = (k1 + 4.0 * k2 + k3) / 6.0;
	return k;
}

double func(double x, double y)
{
	return (3 * y) / (x + 1);  //１階常微分方程式
}

double func_y(double x)
{
	return pow(x + 1, 3.0);  //特殊解
}

実行結果

コード:

 i  x         y        z
 2, 1.000000, 7.661111, 8.000000
v=0.338889
0.338889,n=2
 4, 1.000000, 7.932150, 8.000000
v=0.067850
0.067850,n=4
 8, 1.000000, 7.989216, 8.000000
v=0.010784
0.010784,n=8
16, 1.000000, 7.998481, 8.000000
v=0.001519
0.001519,n=16
32, 1.000000, 7.999799, 8.000000
v=0.000201
0.000201,n=32
64, 1.000000, 7.999974, 8.000000
v=0.000026
0.000026,n=64
128, 1.000000, 7.999997, 8.000000
v=0.000003
続行するには何かキーを押してください . . .

EPS が 1.0e-5 ですから、分割数 128 で OK です。
ε=10^-3 でいいのなら、分割数は 32 ですね。