Modified Euler vs Runge Kutta vs Euler

소스코드

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

double f(double t, double y)
{
	return(y-(t*t)+1);
}
void meuler2()
{
	int i=1;
	float x, y, x1=0, y1=0.5, xf=2, y2, h=0.05;
	x = x1; 
	y = y1;
	printf("\nModified Euler\n");
	printf("%f\t%f\n", x, y);
	while (x<xf)
	{
		y2 = y + h *f(x, y);
		y1 = y2;
		y2 = y + (h / 2)*(f(x, y) + f(x + h, y2));
		y = y2;
		x = x + h;
		if((i%2)==0 && i<11)
		printf("%f\t%f\n", x, y);
		i++;
	}
}

float runge2()
{
	float x0 = 0, y0 = 0.5,	x = 2, h = 0.1;
	printf("\nRunge_Kutta\n");
	float n = (x - x0) / h;
	float k1, k2, k3, k4, k5;

	float y = y0;
	printf("%f\t%f\n", x0, y);
	for (int i = 1; i <= n; i++)
	{
		k1 = h*f(x0, y);
		k2 = h*f(x0 + 0.5*h, y + 0.5*k1);
		k3 = h*f(x0 + 0.5*h, y + 0.5*k2);
		k4 = h*f(x0 + h, y + k3);
		y = y + (1.0 / 6.0)*(k1 + 2 * k2 + 2 * k3 + k4);;
		x0 = x0 + h;
		if(x0<0.6)
			printf("%f\t%f\n", x0, y);
	}
	return y;
}

void euler()
{
	printf("\nEuler\n");
	float a=0, b=.5, x, y, h=0.025, t=0.5, k;
	int tmp=1, i=4;
	x = a;
	y = b;
	printf("%f\t%f\n", x, y);
	while (x <= t)
	{
		k = h*f(x, y);
		y = y + k;
		x = x + h;
		if (tmp == i)
		{
			printf("%f\t%f\n", x, y);
			i += 4;
		}
		tmp++;
	}
}
int main(void)
{
	euler();
	meuler2();
	runge2();
}

결과

소스코드

#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);
}

결과