Runge-Kutta vs Trapezoidal (Method)

소스코드

#include<stdio.h>
#include<math.h>
#define TOL 0.0001
double f(double t, double w);
double y(double t);
double absol(double x);
void main()
{
	int i, k, N = 5, n = 0;
	double t, k1, k2, k3, k4, p1, p2, h = 0.2;
	double w[6], v[5];
	for (k = 1; k <= 2; k++)
	{
		w[0] = -1;
		v[0] = -1;
		p2 = -1;
		for (i = 0; i<N; i++) {
			t = h*i;
			k1 = h*f(t, w[i]);
			k2 = h*f(t + h / 2, w[i] + k1 / 2);
			k3 = h*f(t + h / 2, w[i] + k2 / 2);
			k4 = h*f(t + h, w[i] + k3);
			w[i + 1] = w[i] + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
		}
		for (i = 0; i<N; i++) {
			t = h*i;
			p1 = 100;
			while (absol(p2 - p1) >= TOL)
			{
				p1 = p2;
				p2 = p1 - (p1 - v[i] - h / 2 * (f(t, v[i]) + f(t + h, p1))) /
					(1 - h / 2 * (5 * exp(5 * (t + h))*(2 * p1 - 2 * (t + h))));
			}
			v[i + 1] = p2;
		}
		if (k == 1) {
			printf(" \tRunge-Kutta Method \t\tTrapezoidal Method\n");
			printf(" \t\t h=0.2 \t\t\t\t h=0.2\n");
			printf("t \tw   \t\terror\t\t w\t\terror\n");
			for (i = 0; i <= 5; i++)
			{
				t = 0.2*i;
				printf("%.1f\t%f\t%e\t %f\t%e\n",t, w[i], absol(y(t) - w[i]), v[i], absol(y(t) - v[i]));
			}
		}
		h = 0.25;
		N = 4;
	}
	printf("\n\tRunge-Kutta Method \t\tTrapezoidal Method\n");
	printf(" \t\t h=0.2 \t\t\t\t h=0.2\n");
	printf("t \tw   \t\terror\t\t w\t\terror\n");
	for (i = 0; i <= 4; i++)
	{
		t = 0.25*i;
		if (w[i] == INFINITY) //OVERFLOW 체크
			printf("%.2f\t OVERFLOW\t OVERFLOW\t %f\t%e\n", t, v[i], absol(y(t) - v[i]));
		else	//w[i] = atof("stackoverflow");
			printf("%.2f\t%e\t%e\t %f\t%e\n",t, w[i], absol(y(t) - w[i]), v[i], absol(y(t) - v[i]));
	}
}
double f(double t, double w)
{
	return 5 * exp(5 * t)*(w - t)*(w - t) + 1;
}
double absol(double x)
{
	if (x >= 0)return x;
	else return (-1) * x;
}
double y(double t)
{
	return t - exp(-5 * t);
}

결과