수치미분(테일러시리즈, 리차드슨, 직사각형 수치미분, 사다리꼴 적분, 심프슨의 수치적분)

학교수업, CS/수치해석

수치미분(테일러시리즈, 리차드슨, 직사각형 수치미분, 사다리꼴 적분, 심프슨의 수치적분)

빨대도둑 2022. 11. 30. 15:38

테일러 시리즈를 이용한 수치미분

차분 근사법: i+1번째 항에서 i번째 항을 빼어 미분값을 계산하는 방법

뒤에서 빼는 차분 근사법: i번째 미분값을 계산하는데 i번째 항에서 i-1번째 항을 빼어서 사용하는 방법

중간 차분 근사법: i번째 차분법을 계산하는데 i-1번째 항에서 i-1번째 항을 빼서 계산하는 방법

이러한 3가지 방법의 성능을 배교해 본다면 같은 계산량에 대해서 앞의 두 가지 방법이 h에 비례하는 오차를 갖지만 마지막 방법은 hxh에 비례하는 오차를 가진다. 따라서 중간 차분 근사법이 정확도 면에서 가장 우수하다

<수식> 중간 차분 근사식을 이용하여 함수의 미분 계수를 추정하는 프로그램

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

float eval_f(float x) {
	float y;
	y = x * x + sin(2.0 * x) + 7.0;
	return y;
}

float eval_df(float x) {
	float y;
	y = 2.0 * x + 2.0 * cos(2.0 * x);
	return y;
}

float der1(float x, float h) {
	float fx2, fx1, fx11, fx22, y;

	fx2 = eval_f(x + 2.0 * h);
	fx1 = eval_f(x + h);
	fx11 = eval_f(x - h);
	fx22 = eval_f(x - 2.0 * h);

	y = (-fx2 + 8.0 * fx1 - 8.0 * fx11 + fx22) / (12.0 * h);

	return y;
}

int main(void) {
	float xs, xe, h, tp, x, d, df, err, dfr, f;
	int i, N;

	printf("*type left bound of the interval : ");
	scanf_s("%f", &xs);
	
	printf("*type right bound of the interval : ");
	scanf_s("%f", &xe);

	printf("*type print step interval in x: ");
	scanf_s("%f", &tp);

	printf("*type calculation steop size <h> : ");
	scanf_s("%f", &h);

	N = (int)((xe - xs) / tp);
	x = xs;

	printf("------------------------------------------\n");
	printf("X: F(X) : Esti.Value : Real Value : Error \n");
	printf("------------------------------------------\n");

	for (i = 0; i < N + 1; i++) {
		df = der1(x, h);
		dfr = eval_df(x);
		f = eval_f(x);
		err = fabs((df - dfr) * 100.0 / dfr);
		printf("%10.3f : %11.3f : %11.3f : %11.3f : %8.5f%%\n", x, f, df, dfr, err);
	
		x += tp;
	}printf("------------------------------------------\n");

}

리차드슨 방법

수치 미분에서 오차를 줄이려면 스템 간격 h를 줄이거나, 더 높은 차수까지 고려한 복잡한 식을 사용하여야 한다. 또 다른 오차를 줄이는 방법은 리차드슨 방법이다. 수직 적분에 사용하는 리차드슨 외삽법의 식으로부터 구해진다.

실험하여 수집한 데이터는 경우에 따라서 x축으로 등간격이 아닐 수 도 있다. 이때 데이터를 적절히 곡선 적합화 하거나, 보간법을 사용하여 균등한 간격의 데이터로 만들 수 있다.

<수식>

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

float eval_f(float x) {
	float y;
	y = x * x + sin(2.0 * x) + 7.0;
	return y;
}

float eval_df(float x) {
	float y;
	y = 2.0 * x + 2.0 * cos(2.0 * x);
	return y;
}

float der1(float x, float h) {
	float fx2, fx1, fx11, fx22, dh2,dh1,y;

	fx2 = eval_f(x + 2.0 * h);
	fx1 = eval_f(x + h);
	fx11 = eval_f(x - h);
	fx22 = eval_f(x - 2.0 * h);
	dh1 = (-fx2 + 8.0 * fx1 - 8.0 * fx11 + fx22) / (12.0 * h);

	fx2 = eval_f(x + h);
	fx1 = eval_f(x + 0.5 * h);
	fx11 = eval_f(x - 0.5 * h);
	fx22 = eval_f(x - h);
	dh2 = (-fx2 + 8.0 * fx1 - 8.0 * fx11 + fx22) / (12.0 * h*0.5);


	y = 4.0 * dh2 / 3.0 - dh1 / 3.0;

	return y;
}

int main(void) {
	float xs, xe, h, tp, x, d, df, err, dfr, f;
	int i, N;

	printf("*type left bound of the interval : ");
	scanf_s("%f", &xs);
	
	printf("*type right bound of the interval : ");
	scanf_s("%f", &xe);

	printf("*type print step interval in x: ");
	scanf_s("%f", &tp);

	printf("*type calculation steop size <h> : ");
	scanf_s("%f", &h);

	N = (int)((xe - xs) / tp);
	x = xs;

	printf("------------------------------------------\n");
	printf("X: F(X) : Esti.Value : Real Value : Error \n");
	printf("------------------------------------------\n");

	for (i = 0; i < N + 1; i++) {
		df = der1(x, h);
		dfr = eval_df(x);
		f = eval_f(x);
		err = fabs((df - dfr) * 100.0 / dfr);
		printf("%10.3f : %11.3f : %11.3f : %11.3f : %8.5f%%\n", x, f, df, dfr, err);
	
		x += tp;
	}printf("------------------------------------------\n");

}

직사각형 수치 미분법

수치 적분의 가장 초보적인 형태이며 가장 에러율이 큰 방법이다. 이 방법으로 보다 정확한 값을 구하려면 h의 간격을 매우 좁히는것, 즉 소수 간의 수 n을 매우 늘리는 방법을 사용할 수 있다.

<수식>

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

#define SQ(X) ((X)*(X))

float funct(float x) {
	float out;
	out = exp(-SQ(x));
	return out;
}

int main(void) {
	int i;
	float xl, xh, f_out, h, integ,num;

	printf("\n enter the lower value of x-domain : ");
	scanf_s("%f", &xl);

	while (1)
	{
		printf("\n enter the higher value of x-domain : ");
		scanf_s("%f", &xh);

		if (xh > xl)break;
		else {
			printf("\n Retype");
			continue;
		}
	}

	while (1)
	{
		printf("\n enter the number of division : ");
		scanf_s("%f", &num);

		if (num > 0)break;
		else {
			printf("\n Retype");
			continue;
		}
	}
	h = (xh - xl) / (float)num;

	integ = 0.0;
	for (i = 0; i < num; i++) {
		f_out = funct(xl + (float)i * h);
		integ += h * f_out;
	}
	printf("\n the integrated region =%f\n", integ);
	_getch();
}

사다리꼴 적분법

사다리꼴은 소구간을 직사각형법보다 적분하려는 함수에 더 근접하게 접근함으로써 비교적 정확하게 적분치를 구할 수 있다.

<수식>

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

#define SQ(X) ((X)*(X))

float funct(float x) {
	float out;
	out = exp(-SQ(x));
	return out;
}

int main(void) {
	int i;
	float xl, xh, f_out, h, integ,num;

	printf("\n enter the lower value of x-domain : ");
	scanf_s("%f", &xl);

	while (1)
	{
		printf("\n enter the higher value of x-domain : ");
		scanf_s("%f", &xh);

		if (xh > xl)break;
		else {
			printf("\n Retype");
			continue;
		}
	}

	while (1)
	{
		printf("\n enter the number of division : ");
		scanf_s("%f", &num);

		if (num > 0)break;
		else {
			printf("\n Retype");
			continue;
		}
	}
	h = (xh - xl) / (float)num;

	integ = 0.0;
	for (i = 0; i < num; i++) {
		f_out = funct(xl + (float)i * h);
		integ += h * f_out;
	}

	integ += h * (funct(xl) + funct(xh)) / 2.;

	printf("\n the integrated region =%f\n", integ);
	_getch();
}

심프슨의 수치 적분 1/3

기본적인 알고리즘은 사다리꼴법과 유사하다. 주의할 점은 n은 반드시 2의 배수로 정해야 한다는 것이다. n의 수를 증가시키면 그만큼 정확도도 증가한다.

<수식>

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

#define SQ(X) ((X)*(X))

float funct(float x) {
	float out;
	out = 1./SQ(1.+x);
	return out;
}

int main(void) {
	int i;
	float xl, xh, f_out, h, integ_1,integ_2,integ,num;

	printf("\n enter the lower value of x-domain : ");
	scanf_s("%f", &xl);

	while (1)
	{
		printf("\n enter the higher value of x-domain : ");
		scanf_s("%f", &xh);

		if (xh > xl)break;
		else {
			printf("\n Retype");
			continue;
		}
	}

	while (1)
	{
		printf("\n enter the number of division : ");
		scanf_s("%f", &num);

		if ((num > 0)&&((int)num%2==0)) break;
		else {
			printf("\n Retype");
			continue;
		}
	}
	h = (xh - xl) / (float)num;
	integ_1 = integ_2 = 0.0;


	for (i = 1; i <= num/2; i++) {
		f_out = funct(xl +h*(2.0* (float)i -1.));
		integ_1 += 4.0 * f_out;
	}

	for (i = 2; i <= num / 2; i++) {
		f_out = funct(xl + h * (2.0 * (float)i - 2.));
		integ_2 += 2.0 * f_out;
	}

	integ = h * (integ_1 + integ_2 + funct(xl) + funct(xh)) / 3.;

	printf("\n the integrated region =%f\n", integ);
	_getch();
}

심프슨의 3/8 수치 적분

심프슨의 1/3 적분법과 유사하다. 정확도 면에서는 이 방법이 앞선다. 차이점은 3점과 4점의 차이일 뿐이다. 즉 심프슨의 1/3 법은 3개의 점에 대해 라그랑제 보간 다항식을 이용했으나 심프슨의 3/8 법은 4개의 점에 대해 라그랑제 보간 다항식을 이용하게 된다.

<수식>

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

#define SQ(X) ((X)*(X))

float funct(float x) {
	float out;
	out = 1./SQ(1.+x);
	return out;
}

int main(void) {
	int i;
	float xl, xh, f_out, h, integ_1,integ_2,integ_3, integ,num;

	printf("\n enter the lower value of x-domain : ");
	scanf_s("%f", &xl);

	while (1)
	{
		printf("\n enter the higher value of x-domain : ");
		scanf_s("%f", &xh);

		if (xh > xl)break;
		else {
			printf("\n Retype");
			continue;
		}
	}

	while (1)
	{
		printf("\n enter the number of division : ");
		scanf_s("%f", &num);

		if ((num > 0)&&((int)num%2==0)) break;
		else {
			printf("\n Retype");
			continue;
		}
	}
	h = (xh - xl) / (float)num;
	integ_1 = integ_2 = integ_3 = 0.0;


	for (i = 2; i <= num/3; i++) {
		f_out = funct(xl +h*(3.0* (float)i -3.));
		integ_1 += 2.0 * f_out;
	}

	for (i = 1; i <= num / 3; i++) {
		f_out = funct(xl + h * (3.0 * (float)i - 2.));
		integ_2 += 3.0 * f_out;
	}

	for (i = 1; i <= num / 3; i++) {
		f_out = funct(xl + h * (3.0 * (float)i - 1.));
		integ_3 += 3.0 * f_out;
	}

	integ = h * (integ_1 + integ_2 +integ_3+ funct(xl) + funct(xh)) * 3./8.;

	printf("\n the integrated region =%f\n", integ);
	_getch();
}

'학교수업, CS > 수치해석' 카테고리의 다른 글

편미분 방정식(PDE, 라플라스, 포아송) (0)	2022.11.30
상미분 방정식(오일러 방정식, 룬게쿤타방법, 투사법, 유한차분법) (0)	2022.11.30
푸리에 변환(DFT, FFT, 파워스펙트럼 추정, 자기상관도) (0)	2022.11.30
난수(변환법, 거부법, 가우시안 분포, 몬테카를로 적분) (0)	2022.11.30
보간법(선형, 랑그라제, 네빌레, 뉴튼 다항식, 전/후향 차분법, 스플라인) (0)	2022.11.30

현재글수치미분(테일러시리즈, 리차드슨, 직사각형 수치미분, 사다리꼴 적분, 심프슨의 수치적분)

brilliant-star