Microsoft MVP성태의 닷넷 이야기
글쓴 사람
홈페이지
첨부 파일

그래프 그리기로 알아보는 뉴턴-랩슨(Newton-Raphson's method)법과 제곱근 구하기 - C#

뉴턴-랩슨법은 방정식의 근사해를 구하는 방법인데요,

뉴턴법/뉴턴-랩슨법의 이해와 활용(Newton's method)
; http://darkpgmr.tistory.com/58

이론이 꽤나 직관적입니다.

간단하게 프로그램을 만들어 가면서 예를 들어볼까요? ^^ 다음의 2차 방정식으로 시작해 보겠습니다.

y = x2 - 10

그래프를 그리기 위해 우선 x, y 기준선을 그려보겠습니다.

using System;
using System.Drawing;
using System.Windows.Forms;

namespace WindowsFormsApplication1
{
    public partial class Form1 : Form
    {
        Point _ptCenter;

        public Form1()
        {
            InitializeComponent();
            SetStyle(ControlStyles.ResizeRedraw | ControlStyles.UserPaint | ControlStyles.AllPaintingInWmPaint, true);
            UpdateStyles();
        }

        private void Form1_Load(object sender, EventArgs e)
        {
            this.Width = 800;
            this.Height = 800;
        }

        protected override void OnPaintBackground(PaintEventArgs e)
        {
            e.Graphics.Clear(Color.White);
        }

        protected override void OnPaint(PaintEventArgs e)
        {
            Graphics g = e.Graphics;

            _ptCenter = DrawXYBase(g);
        }
        
        private Point DrawXYBase(Graphics g)
        {
            int centerX = this.Width / 2;
            int centerY = this.Height / 2;

            g.DrawLine(Pens.Black, 0, centerY, this.Width, centerY);
            g.DrawLine(Pens.Black, centerX, 0, centerX, this.Height);

            return new Point(centerX, centerY);
        }
    }
}

우리가 원하는 2차 방정식은 다음과 같이 간단하게 람다식으로 표현할 수 있습니다.

Func<float, float> orgFunc = (x) => (float)(x * x - 10);

x의 증가에 따라 구해지는 y 값에 따라 (x, y) 좌표에 점을 찍어나가는 방법으로 그래프를 표현할 건데요. 아쉽게도 x2의 큰 차이로 y 값이 크게 증가하기 때문에 x 변위로 그래프를 그리기에는 좀 맞지 않습니다. 따라서 다음의 역함수를 구해,

Func<float, float> inversedF = (y) => (float)Math.Sqrt(y + 10);

거꾸로 y 값을 바꿔가면서 x 좌표를 구할 것입니다. 이 함수를 이용해 그래프를 그리면,

protected override void OnPaint(PaintEventArgs e)
{
    Graphics g = e.Graphics;

    _ptCenter = DrawXYBase(g);

    Func<float, float> orgFunc = (x) => (float)(x * x - 10);
    Func<float, float> inversedF = (y) => (float)Math.Sqrt(y + 10);

    DrawGraph(inversedF, g);
}

private void DrawGraph(Func<float, float> inversedF, Graphics g)
{
    for (float y = -_ptCenter.Y; y < 0; y ++)
    {
        float x = inversedF(y);

        if (float.IsNaN(x) == true)
        {
            continue;
        }

        DrawPixel(g, Brushes.Black, x, y);
        DrawPixel(g, Brushes.Black, -x, y);
    }

    for (float y = 0; y < _ptCenter.Y; y ++)
    {
        float x = inversedF(y);

        if (float.IsNaN(x) == true)
        {
            continue;
        }

        DrawPixel(g, Brushes.Black, x, y);
        DrawPixel(g, Brushes.Black, -x, y);
    }
}

private void DrawPixel(Graphics g, Brush brush, float x, float y)
{
    PointF pt = Projection(x, y);
    g.FillRectangle(brush, pt.X, pt.Y, 1, 1);
}

private PointF Projection(float x, float y)
{
    return new PointF(x + _ptCenter.X, -y + _ptCenter.Y);
}

실행 결과가 이렇게 구해집니다.

newton_raphson_1.png

너무 멋이 없군요. ^^; 그래서 좀더 극적인 효과를 높이기 위해 그래프의 기울기를 완만하게 할 수 있도록 1보다 작은 임의의 값을 곱하겠습니다.

float _slopeForGraph = 70.0f;

Func<float, float> orgFunc = (x) => (float)(1 /_slopeForGraph * x * x - 10);
Func<float, float> inversedF = (y) => (float)Math.Sqrt((y + 10) * _slopeForGraph);

이렇게 하면 좀더 볼만하게 바뀌는데,

newton_raphson_2.png

여전히 아쉽게도 하단의 영역이 뚫린 영역이 좀 거슬립니다. 이를 보완하기 위해 y의 증감을 1에서 0.1 정도로 낮춰주면,

float _incOffset = 0.1f;

private float DrawTangentLine(Func<float, float> slopeFunc, Graphics g)
{
    float minY = 100;
    float minX = 0;

    for (float x = -_ptCenter.X; x < 0; x += _incOffset)
    {
        // ...[생략]...
    }

    for (float x = 0; x < _ptCenter.X; x += _incOffset)
    {
        // ...[생략]...
    }

    return minX;
}

이제야 좀 그래프다운 모습이 나옵니다. ^^

newton_raphson_3.png




자, 이제 x 좌표의 아무 값이나 (하지만 극적인 효과를 위해 화면에 보일 정도의 값으로) 선택해서 (x * x - 10) 방정식에 대입해 나오는 y 좌표까지의 선을 연결합니다.

protected override void OnPaint(PaintEventArgs e)
{
    // ...[생략]...

    DrawGraph(inversedF, g);

    float startX = 150.0f;

    DrawVerticalLine(orgFunc, g, startX);
}

private void DrawVerticalLine(Func<float, float> orgFunc, Graphics g, float x)
{
    float y = orgFunc(x);

    for (float y1 = 0; y1 < y; y1++)
    {
        DrawPixel(g, Brushes.Blue, x, y1);
    }
}

그럼, 이렇게 나오겠죠?!

newton_raphson_4.png

바로 저렇게 수직선과 2차 방정식 그래프가 만난 점의 접선을 구해 그려보겠습니다. 접선은 1차 방정식이므로 다음과 같은 유형으로 나와야 합니다.

y = a * x + b

여기서 접선의 기울기 a 값은 수직선이 만난 2차 방정식 그래프 지점의 기울기와 같습니다. 따라서 원래의 2차 방정식을 미분해 "y = 2 * x" 식을 구할 수 있고 여기에 x 값을 수직선의 x 값으로 대체해 주면 됩니다. 따라서 "y = a * x + b"의 a는 다음의 식으로 구할 수 있습니다.

float startX = 150.0f; // 임의로 정한 x == 150.0f;

float lineSlope = (1 / _slopeForGraph) * 2.0f * startX;

a 값을 구했으니 "y = a * x + b"에서의 b 값을 구해야 하는데요. 역시 2차 방정식에 x를 대입했을 때 나오는 y 값의 (x, y) 쌍이 접선 방정식에서도 그대로 구해져야 한다는 점을 이용해 다음과 같이 구할 수 있습니다.

Func<float, float> orgFunc = (x) => (float)(1 /_slopeForGraph * x * x - 10);

float startX = 150.0f;

float lineSlope = (1 / _slopeForGraph) * 2.0f * startX;
float slopeBParam = (orgFunc(startX) - (lineSlope * startX));

자, 그럼 이제 접선의 방정식을 다음과 같이 구성할 수 있고,

protected override void OnPaint(PaintEventArgs e)
{
    // ...[생략]...
    Func<float, float> slopeLineFunc = (x) => (lineSlope * x + slopeBParam);

    DrawTangentLine(slopeLineFunc, g);
}

void DrawTangentLine(Func<float, float> slopeFunc, Graphics g)
{
    for (float x = -_ptCenter.X; x < 0; x += _incOffset)
    {
        float y = slopeFunc(x);
        DrawPixel(g, Brushes.Red, x, y);
    }

    for (float x = 0; x < _ptCenter.X; x += _incOffset)
    {
        float y = slopeFunc(x);
        DrawPixel(g, Brushes.Red, x, y);
    }
}

이를 그려서 나온 직선을 통해 다시 y = 0이 되는 지점의 x 값을 구할 수 있습니다.

newton_raphson_5.png

이 x 값을 구하기 위해 slopeLineFunc 접선의 방정식을 역함수로 바꾸면 되는데,

Func<float, float> slopeLineFunc = (x) => (lineSlope * x + slopeBParam);
Func<float> inverseSlopeLineFunc = (y) => ((y - slopeBParam) / lineSlope);

어차피 y 좌표의 값이 0인 순간에 대한 x 값을 구하는 것이기 때문에 다음과 같이 간략화할 수 있습니다.

Func<float> inverseSlopeLineFunc = () => (-slopeBParam / lineSlope);

float newX = inverseSlopeLineFunc();

자... 그런데, 재미있게도 새롭게 구한 newX를 가장 처음에 잡았던 임의의 x 값(150.0f) 대신 사용해 다시 수직선을 긋고 그것의 접선을 다시 구한 다음 또 newX 값을 구할 수 있게 됩니다. 프로그램도 간단합니다. 단지 위의 코드를 2번만 반복해 주면 됩니다.

float startX = 150.0f;

for (int i = 0; i < 2; i++) // 2번 반복
{
    DrawVerticalLine(orgFunc, g, startX);

    float lineSlope = (1 / _slopeForGraph) * 2.0f * startX;

    float slopeBParam = (orgFunc(startX) - (lineSlope * startX));

    Func<float, float> slopeLineFunc = (x) => (lineSlope * x + slopeBParam);
    DrawTangentLine(slopeLineFunc, g);

    Func<float> inverseSlopeLineFunc = () => (-slopeBParam / lineSlope);

    float newX = inverseSlopeLineFunc();

    startX = newX; // 마지막 x 값을 다시 피드백받아 진행
}

그래프 모양은 이렇습니다.

newton_raphson_6.png

for 문의 반복을 늘려보면 자연스럽게 y = 0인 지점의 접선으로 결과가 수렴하는 것을 볼 수 있습니다.

newton_raphson_7.png




이런 성격을 가진 뉴턴-랩슨법이 제곱근을 구하는 데 사용하는 것도 역시 직관적입니다. 예를 들어 2의 제곱근의 값이 x라고 할 때 이는 x를 제곱해서 2가 나온다는 것입니다. 따라서 다음의 방정식의 해를 구하는 것과 같고,

y = x2 - 2

0 = x2 - 2 // (이때의 x 값은 2의 제곱근!)

y = 0일 때이니, 바로 위에서 "for 문의 반복을 늘려보면 자연스럽게 y = 0인 지점의 접선으로 결과가 수렴하는 것을 볼 수 있습니다." 라고 말한 내용이 결국 제곱근과 연관되어 있음을 알 수 있습니다.

간단하지요? ^^

그래서 결국 위에서 그래프를 그리기 위해 사용한 공식을 (멋진 그래프를 위해 보정한 수치를 빼면) 그대로 사용할 수 있습니다. 이를 이용해 원하는 수의 제곱근을 구하기 위한 일반적인 함수를 다음과 같이 만들 수 있습니다.

static double newtonRaphson3(double xsquare, double x0, int depth = 1)
{
    Func<double, double> orgFunc = (x) => (double)(x * x - xsquare);

    for (int i = 0; i < depth; i++)
    {
        double lineSlope = 2.0f * x0;

        double slopeBParam = (orgFunc(x0) - (lineSlope * x0));

        Func<double> inverseSlopeLineFunc = () => (-slopeBParam / lineSlope);

        double newX = inverseSlopeLineFunc();

        x0 = newX; // 마지막 x 값을 다시 피드백받아 진행
    }

    return x0;
}

위의 코드를 이용해 닷넷에서 기본제공되는 Math.Sqrt와 값을 비교해 테스트할 수 있습니다.

double x = 10; // 10의 제곱근을 구한다.
double x0 = x / 2; // 초기값은 5로.

Console.WriteLine(string.Format("{0}", newtonRaphson3(x, x0, 25)));
Console.WriteLine(string.Format("{0}", Math.Sqrt(x)));

x = 2; // 2의 제곱근을 구한다.
x0 = x / 2; // 초기값은 1로.

Console.WriteLine(string.Format("{0}", newtonRaphson3(x, x0, 25)));
Console.WriteLine(string.Format("{0}", Math.Sqrt(x)));

// 출력 결과
3.16227766016838
3.16227766016838
1.41421356237309
1.4142135623731

마지막으로, newtonRaphson3 메서드의 내부 코드를 좀더 치환해서 정리하면 다음과 같이 간단한 함수 하나로 남습니다.

static double newtonRaphson4(double xsquare, double x0, int depth = 1)
{
    for (int i = 0; i < depth; i++)
    {
        x0 = (x0 * x0 + xsquare) / (2.0f * x0);
    }

    return x0;
}

오~~~ 어디서 많이 보던 공식 아닌가요? ^^

뉴턴-랩슨 법(Newton-Raphson method)
; http://tro.kr/34




끝맺기 전에, 최근에 본 글에서 재미있는 공식을 봤습니다.

왜 함수형 프로그래밍이 중요한가
; https://drive.google.com/file/d/0B01vnjtP-BZsSXBnLXlfSW4zYU9zS2pXZFB4YXpOa0t0N293/view

위의 글에 보면 뉴튼-랩슨 제곱근 구하는 공식으로 다음을 소개하고 있습니다.

ai+1 = (ai + n / ai) / 2

그런데, 이 공식의 유추 과정이 재미있습니다. ^^ 미분 및 접선의 방정식을 구할 필요도 없이 그냥 다음과 같은 사고로 접근하고 있습니다.







간단하죠? ^^ 이렇게 제곱근을 구하는 방식을 "바빌로니아 법"이라 한다고! ^^

바빌로니아 법
; https://ko.wikipedia.org/wiki/%EB%B0%94%EB%B9%8C%EB%A1%9C%EB%8B%88%EC%95%84_%EB%B2%95

위의 공식에 초기값인 a를 분모/분자에 곱해주면 그것이 바로 뉴턴-랩슨법의 공식과 바로 일치합니다.

끝!

(첨부 파일은 위의 그래프를 그리는 C# 윈폼 프로젝트와 newtonRaphson4 메서드를 테스트하는 콘솔 프로젝트를 포함합니다.)




[이 글에 대해서 여러분들과 의견을 공유하고 싶습니다. 틀리거나 미흡한 부분 또는 의문 사항이 있으시면 언제든 댓글 남겨주십시오.]

[연관 글]





[최초 등록일: ]
[최종 수정일: 2/29/2016 ]

Creative Commons License
이 저작물은 크리에이티브 커먼즈 코리아 저작자표시-비영리-변경금지 2.0 대한민국 라이센스에 따라 이용하실 수 있습니다.
by SeongTae Jeong, mailto:techsharer@outlook.com

비밀번호

댓글 쓴 사람
 




[1]  2  3  4  5  6  7  8  9  10  11  12  13  14  15  ...
NoWriterDateCnt.TitleFile(s)
12266정성태7/14/202028오류 유형: 630. 사용자 계정을 지정해 CreateService API로 서비스를 등록한 경우 "Error 1069: The service did not start due to a logon failure." 오류발생
12265정성태7/10/202050오류 유형: 629. Visual Studio - 웹 애플리케이션 실행 시 "Unable to connect to web server 'IIS Express'." 오류 발생
12264정성태7/9/202040오류 유형: 628. docker: Error response from daemon: Conflict. The container name "..." is already in use by container "...".
12261정성태7/9/2020198VS.NET IDE: 148. 윈도우 10에서 .NET Core 응용 프로그램을 리눅스 환경에서 실행하는 2가지 방법 - docker, WSL 2 [3]
12260정성태7/8/202084.NET Framework: 926. C# - ETW를 이용한 ThreadPool 스레드 감시파일 다운로드1
12259정성태7/8/202041오류 유형: 627. nvlddmkm.sys의 BAD_POOL_HEADER BSOD 문제
12258정성태7/8/2020106기타: 77. DataDog APM 간략 소개
12257정성태7/7/202076.NET Framework: 925. C# - ETW를 이용한 Monitor Enter/Exit 감시파일 다운로드1
12256정성태7/7/2020121.NET Framework: 924. C# - Reflection으로 변경할 수 없는 readonly 정적 필드 [4]
12255정성태7/6/202086.NET Framework: 923. C# - ETW(Event Tracing for Windows)를 이용한 Finalizer 실행 감시파일 다운로드1
12254정성태7/2/202055오류 유형: 626. git - REMOTE HOST IDENTIFICATION HAS CHANGED!
12253정성태7/2/2020135.NET Framework: 922. C# - .NET ThreadPool의 Local/Global Queue파일 다운로드1
12252정성태7/2/2020113.NET Framework: 921. C# - I/O 스레드를 사용한 비동기 소켓 서버/클라이언트파일 다운로드2
12251정성태7/1/2020137.NET Framework: 920. C# - 파일의 비동기 처리 유무에 따른 스레드 상황파일 다운로드2
12250정성태7/1/2020368.NET Framework: 919. C# - 닷넷에서의 진정한 비동기 호출을 가능케 하는 I/O 스레드 사용법 [1]파일 다운로드1
12249정성태6/29/202055오류 유형: 625. Microsoft SQL Server 2019 RC1 Setup - 설치 제거 시 Warning 26003 오류 발생
12248정성태6/29/202057오류 유형: 624. SQL 서버 오류 - service-specific error code 17051
12247정성태6/29/2020157.NET Framework: 918. C# - 불린 형 상수를 반환값으로 포함하는 3항 연산자 사용 시 단축 표현 권장(IDE0075) [2]파일 다운로드1
12246정성태6/29/202095.NET Framework: 917. C# - USB 관련 ETW(Event Tracing for Windows)를 이용한 키보드 입력을 감지하는 방법
12245정성태6/25/2020265.NET Framework: 916. C# - Task.Yield 사용법 (2) [2]파일 다운로드1
12244정성태6/29/2020115.NET Framework: 915. ETW(Event Tracing for Windows)를 이용한 닷넷 프로그램의 내부 이벤트 활용파일 다운로드1
12243정성태6/23/202069VS.NET IDE: 147. Visual C++ 프로젝트 - .NET Core EXE를 "Debugger Type"으로 지원하는 기능 추가
12242정성태6/24/202041오류 유형: 623. AADSTS90072 - User account '...' from identity provider 'live.com' does not exist in tenant 'Microsoft Services'
12241정성태6/26/2020131.NET Framework: 914. C# - Task.Yield 사용법파일 다운로드1
12240정성태6/23/2020109오류 유형: 622. 소켓 바인딩 시 "System.Net.Sockets.SocketException: An attempt was made to access a socket in a way forbidden by its access permissions" 오류 발생
[1]  2  3  4  5  6  7  8  9  10  11  12  13  14  15  ...