학교수업, CS/수치해석

상미분 방정식(오일러 방정식, 룬게쿤타방법, 투사법, 유한차분법)

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

상미분을 분류하는 기준

  • 2차 혹은 고차
  • 선형 혹은 비선형
  • 초기 조건 혹은 경계 조건
  1. 첫째 기준인 1차 혹은 고차는 단순히 방정식 중 미분향의 최고차항의 차수가 1차인지 그 이상인지에 따라 나눈 기준으로 보통 그 차수의 개수만큼의 조건식이 주어진다
  2. 선형 혹은 비현형의 분류는 방정식 중 미분향의 계수를 보고 따지는 것인데 그 계수가 y에 관한 함수이면 비선형으로 보고 상수나 혹은 x에 관한 함수며 선형으로 볼 수 있다. 선형 비선형 모두 식을 푸는 알고리즘 면에서는 큰 차이가 없으나 계산량과 복잡도 면에서 비선형 미방을 푸는 것이 더 복잡하다
  3. 초기 조건과 경계 조건의 분류인데 고차 미방에서의 값이 각각주어지면 초기 조건의 문제라고 불 수 있다. 일반적으로 경계조건의 문제가 더 복잡하며 또 초기 조건의 문제를 알아야만 풀 수 있다.

 

오일러 방법

  1. 구하고자 하는 구간을 x0부터 x1으로 하여 받아들인다
  2. 초기 조건 y0을 받아들인다
  3. 인쇄할 간격 tp로 받아들인다. 결과적으로 (x1-x0)/tp개 만큼만 계산한다
  4. x축의 계산 간격을 h에 받아들인다. 

<수식>

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

#define N 500

float func(float x) {
	float a;
	a = 3.5 * x * x - 4.7 * x + 10.5;
	return a;
}

float FX(float x, float ic) {
	float a;
	a = 3.5 / 3.0 * x * x * x - 4.7 / 2.0 * x * x + 10.5 * x + ic;
	return a;
}

int main(void) {
	int i, j, n;
	float xp[N], yp[N], x, y, x0, x1, y0, h, tp, nc, np, fr;

	printf("*start value in x: ");
	scanf_s("%f", &x0);

	printf("*end value in x : ");
	scanf_s("%f", &x1);

	printf("*initial value in y : ");
	scanf_s("%f", &y0);

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

	np = (x1 - x0) / tp;
	if (np > 500.0) {
		printf("reduce integral interval or print step interval\n");
		exit(1);
	}

	printf("\n *step size in x : ");
	scanf_s("%f", &h);

	nc = tp / h;

	x = x0; y = y0; 
	xp[0] = x;	yp[0] = y;

	for (i = 0; i < (int)np; i++) {
		for (j = 0; j < (int)nc; j++) {
			y = y + func(x) * h;
			x = x + h;
		}
		xp[i + 1] = x;
		yp[i + 1] = y;
	}
	printf("result \n");
	printf("X | Y | REAL VALUE | ERROR \n");
	printf("----------------------------------\n");

	for (i = 0; i <= (int)np; i++) {
		fr = FX(x0 + tp * (float)i, y0);
		printf("%12.5f | %12.5f | %12.5f | %12.5f \n",xp[i],yp[i],fr,fabs((yp[i]-fr)/fr));
	}
	_getch();
}

 

변형된 오일러 방정식(룬게 쿠타 방법)

기본적인 개념과 사용하는 항들은 이전의 오일러 방식과 같다. 다만 기울기 S를 추정하는 데 있어서 더욱 정밀한 방법을 사용하여 정밀도를 높인다. 

룬게 쿠타 방법

룬게룬게 쿠타 방법은 테일러 전개식에서 몇 번째 항까지 고려하느냐에 따라서 2차 룬게 쿠타, 3차 룬게 쿠타 방법으로 유도할 수 있다.

<수식>

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

#define N 500

float func(float x) {
	float a;
	a = 3.5 * x * x - 4.7 * x + 10.5;
	return a;
}

float FX(float x, float ic) {
	float a;
	a = 3.5 / 3.0 * x * x * x - 4.7 / 2.0 * x * x + 10.5 * x + ic;
	return a;
}

int main(void) {
	int i, j;
	float xp[N], yp[N], x, y, x0, x1, y0, h, tp, nc, np, fr;
	float k1, k2, p1, a1, a2, err;

	printf("-------------------------------------\n");
	printf("*program for second order funge kutta method* \n");
	printf("-------------------------------------\n");

	printf("*start value in x: ");
	scanf_s("%f", &x0);

	printf("*end value in x : ");
	scanf_s("%f", &x1);

	printf("*initial value in y : ");
	scanf_s("%f", &y0);

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

	np = (x1 - x0) / tp;

	if (np > 500.0) {
		printf("reduce integral interval or print step interval\n");
		exit(1);
	}

	printf("\n *step size in x : ");
	scanf_s("%f", &h);

	printf("-------------------------------------\n");
	printf("select solving method\n");
	printf("1. heun' method with a single corrector \n");
	printf("2. the improved polygon method\n");
	printf("3. ralston's method \n");
	printf("type the number that you want to solve with...");
	scanf_s("%d", &i);

	printf("-------------------------------------\n");

	switch (i)
	{case(1):
		a1 = 0.5;
		a2 = 0.5;
		p1 = 1.0;
		break;
	case(2):
		a1 = 0.0;
		a2 = 1.0;
		p1 = 0.5;
		break;
	case 3:
		a1 = 1.0 / 3.0;
		a2 = 2.0 / 3.0;
		p1 = 3.0 / 4.0;
		break;
	default: printf("type the number one of above...");
		exit(1);
	}


	nc = tp / h;

	x = x0; y = y0; 
	xp[0] = x;	yp[0] = y;

	for (i = 0; i < (int)np; i++) {
		for (j = 0; j < (int)nc; j++) {
			k1 = func(x);
			k2 = func(x + p1 * h);
			y = y + (a1 * k1 + a2 * k2) * h;
			x = x + h;
		}
		xp[i + 1] = x;
		yp[i + 1] = y;
	}

	err = 0.0;

	printf("-------------------------------------\n");
	printf("result \n");
	printf("-------------------------------------\n");
	printf("X | Y | REAL VALUE | ERROR \n");
	printf("----------------------------------\n");

	for (i = 0; i <= (int)np; i++) {
		fr = FX(x0 + tp * (float)i, y0);
		printf("%12.5f | %12.5f | %12.5f | %12.5f \n",xp[i],yp[i],fr,fabs((yp[i]-fr)/fr));
	}
	_getch();
}

 

투사법

2차 이상의 미분 방정식은 그 차수만큼의 조건이 주어져야만 특수해를 구했을 때 나오는 상수의 값을 결정할 수 있었다. 이번에는 초기 조건이 아닌 경계 조건이 주어진 경우 해를 구하는 방법에 대해 알아볼 것이다. 

만일 어떤 미방을 만족시키는 x에 관한 2개의 개별적인 함수가 있다면 그것들의 선형적인 조합이 또한 해가 된다.

<수식>

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

void clrscr(void)
{
	COORD Cur = { 0, 0 };
	unsigned long dwLen;
	FillConsoleOutputCharacter(GetStdHandle(STD_OUTPUT_HANDLE), ' ', 80 * 25, Cur, &dwLen);
}

float funct(float t,float y, float dy) {
	float out;
	out = -2. / t * dy + 2. / (t * t) * y + sin(log(t)) / (t * t);
	return out;
}

int main(void) {
	int try, i, num;
	float x[2], fx[2], h, time[101], y[3][101], dy[3][101], f_out, k[3][101], kk[3][101];

	while (1)
	{
		clrscr();
		printf("-------------- shooting method-----------------\n");
		printf("\n enter the value fo division number between defind region (below 100) : ");
		scanf_s("%d", &num);

		if (num > 0 && num <= 100)
			break;
		else {
			printf("retype | press any key");
			_getch();
			continue;
		}
	}

	/*initialization*/

	x[0] = 1.0;
	fx[0] = 1.0;
	x[1] = 2.0;
	fx[1] = 2.0;

	h = (x[1] - x[0]) / (float)num;
	time[0] = x[0];

	/*algorithm implementation*/

	for (try = 0; try < 3; try++) {
		y[try][0] = fx[0];	/*initialization*/
		if (try != 2) {	/*for tiral =0 or 1*/
			printf("\n guess the value of y(%f) : ", x[0]);
			scanf_s("%f", &dy[try][0]);
		}
		else /*for finding solution*/ {
			printf("\n type any key to continue");
			_getch();
			clrscr();

			printf("----------------------------------\n");
			printf(" the correct solution by shotting method is \n");
			printf("----------------------------------\n");

			dy[try][0] = dy[0][0] + (dy[1][0] - dy[0][0]) / (y[1][num] - y[0][num]) * 
            (fx[1] - y[0][num]);
		}
		printf("\n\t time : % 4.2f \t the value of y[%4.2f]=%10.6f", time[0], time[0], 
        y[try][0]);
		/*second order runge-kutta method implementation*/

		for (i = 0; i < num; i++) {
			f_out = funct(time[i], y[try][i], dy[try][i]);
			k[try][i] = h * h / 2. * f_out;
			f_out = funct(time[i] + 2. / 3. * h, y[try][i] + 2. / 3. * h * dy[try][i] + 2. / 
            3. * k[try][i], dy[try][i] + 4. / 3. / h * k[try][i]);
			kk[try][i] = h * h / 2. * f_out;

			y[try][i + 1] = y[try][i] + h * dy[try][i] + (k[try][i] + kk[try][i]) / 2.0;
			dy[try][i + 1] = dy[try][i] + (k[try][i] + 3. * kk[try][i]) / 2. / h;

			time[i + 1] = time[i] + h;
			printf("\n\t time : %4.2f \t the value of y [%4.2f]=%10.6f", time[i + 1], 
            time[i + 1], y[try][i + 1]);
		}
		printf("\n");
	}
	_getch();

}

 

유한 차분법

경계 조건이 주어지는 미방을 푸는 데는 투사법 말고 유한 차분법이 있다. 유한 차분법은 각 미분항을 유한 차분으로 근사화하여 푸는 것으로서 안정도 면에서 투사법보다 좋은 이점을 갖는다. 

<수식>

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

#define P(X) ((-2.)/(X))
#define Q(X) ((-2.)/(X)/(X))
#define R(X) (sin(log(X))/(X)/(X))

void clrscr(void)
{
	COORD Cur = { 0, 0 };
	unsigned long dwLen;
	FillConsoleOutputCharacter(GetStdHandle(STD_OUTPUT_HANDLE), ' ', 80 * 25, Cur, &dwLen);
}

float tridia(float a[], float b[], float c[], float d[], int num) {
	int i;
	float temp;
	for (i = 2; i <= num; i++) {
		temp = c[i] / a[i - 1];
		a[i] -= temp * b[i - 1];
		d[i] -= temp * d[i - 1];
	}
	d[num] /= a[num];
	for (i = num - 1; i >= 1; i--)
		d[i] = (d[i] - b[i] * d[i + 1] / a[i]);
}

int main(void) {
	int num, i;
	float a[100], b[100], c[100], d[100], temp, y[2], h;

	clrscr();
	printf("-----------------FINTE DIFFERANCE METHOD------------------\n");
	printf("enter the number of iteration: ");
	scanf_s("%d", &num);

	num--;
	a[0] = y[0] = 1.;
	b[0] = y[1] = 2.;

	h = (b[0] - a[0]) / ((float)num + 1.0);
	temp = a[0] + h;

	a[1] = 2.0 + h * h * Q(temp);
	b[1] = -1.0 + h * P(temp) / 2.0;
	d[1] = -h * h * R(temp) + (1.0 + h - P(temp) / 2.) * y[0];

	for (i = 2; i < num; i++) {
		temp = a[0] + (float)i * h;
		a[i] = 2.0 + h * h * Q(temp);
		b[i] = -1.0 + h * P(temp) / 2.0;
		c[i] = -1.0 - h * P(temp) / 2.0;
		d[i] = -h * h * R(temp);
	}

	temp = b[0] - h;

	a[num] = 2.0 + h * h * Q(temp);
	c[num] = -1.0 - h * P(temp) / 2.0;
	d[num] = -h * h * P(temp) + (1. - h * P(temp) / 2.) * y[1];
	d[0] = y[0];
	d[num + 1] = y[1];

	tridia(a, b, c, d, num);

	for (i = 0; i <= num + 1; i++) {
		temp = a[0] + (float)i * h;
		printf("\n x: %4.2f \t W[%4.2f] : %f", temp, temp, d[i]);
	}
	_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

'학교수업, CS/수치해석'의 다른글

  • 현재글상미분 방정식(오일러 방정식, 룬게쿤타방법, 투사법, 유한차분법)

관련글

티스토리툴바