미분방정식 3가지 기법 비교

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

결과