간단한 무리방정식의 풀이

무리방정식에는 무연근이 있을 수 있다.
여기서는 무리방정식

              sqrt(x + 1) = x - 2

의 풀이를 Mactimatica, Maxima 등에서는 어떻게 하는지 알아 본다.

우선 손계산으로 풀어보자.
무리방정시에는 무연근(계산 상으로는 근으로 나오지만 원래 방정식의 의마에는 맞지 얺는 근)이 나올 수 있다. 먼저 근이 존재할 범위를 알아 본다. 제곱근 기호 속이 0 이상이라야 하므로 x + 1 >= 0 즉 x >= -1 이다. 또 제곱근의 값 자체가 0 이상이므로 주어진 방정식의 우변도 0 이상이라야 한다. 그러므로 x - 2 >= 0 즉, x >= 2 이다. 그러므로 근의 존재 범위는 x >= -1 과 x >= 2 의 공통범위 x >= 2 이다. 이 범위에 있지 않는 근은 모두 무연근이므로 버린다.

이제 주어진 방정식의 양변을 제곱하면
              x + 1 = (x - 2)^2
        즉,  x + 1 = x^2 - 4x + 4
이것은 단순한 이차방정식이므로 내림차 순으로 잘 정리하여 근의 공식을 써서 풀면 된다.
             x^2 - 5x + 3 = 0
따라서    x = (5 ± sqrt(13))/2    이다.
이 중에 (5 -sqrt(13))/2 는 2 보다 작으므로 (근의 존재 범위에 맞지 않는) 무연근이라서 버린다.
그러므로 주어진 무리방정식의 해는  x = (5 + sqrt(13))/2  뿐이다.


* Mathematica 를 이용한 무리방정식의 풀이
  (Solve 는 심볼 계산으로 해를 구해주고, NSolve 는 수치적 계산으로 해를 구래준다.)

* Maxima 를 이용한 무리방정식의 풀이
  (Maxima 로는 무리반정식의 해를 심볼로 구하는 방법이 없다.
   그대신 find_root 를 사용하면 수치적 해를 구해준다.
   solve 를 써서 다항방정시의 해를 구한 다음 무연근을 버리는 방법으로 그 해를 구하였다.
   numer 는 심볼 값을 수치적 값으로 계산해 주는 명령이다. *)

* Python 에서 sympy 를 이용한 방정식의 풀이
  (* 아래에서 진한 부분만 입력하면 된다. 
      sympy 로도 무리방정식을 풀지 못해 이차반정식으로 고쳐서 풀었다. *)

>>> from sympy import *
>>> x = symbols('x')
>>> y = solve((x + 1) - (x - 2)**2, x)
>>> y
[-13**(1/2)/2 + 5/2, 13**(1/2)/2 + 5/2]
>>> y[0] >= 2
False
>>> y[1] >= 2
True
>>> y[1]
13**(1/2)/2 + 5/2
>>> N(y[1])
4.30277563773199

* Python 에서 scipy 를 이용하여 무리방정식의 수치적 해 구하기
  (* numpy 와 scipy 는 수치적 계산을 위한 파이썬 전용 라이브러리이다. *)

>>> from scipy import *
>>> from scipy.optimize import fsolve
>>> dir(fsolve)

>>> help (fsolve)

>>> def f(x): return sqrt(x + 1) - x + 2
...
>>> fsolve(f, 0)
array([ 4.30277564])
>>> fsolve(f, 9)
array([ 4.30277564])

또는

>>> from scipy import sqrt
>>> from scipy.optimize import fsolve
>>> fsolve(lambda x: sqrt(x + 1) - x + 2, 1)
array([ 4.30277564])
>>> sol = fsolve(lambda x: sqrt(x + 1) - x + 2, 1)
>>> sol
array([ 4.30277564])
>>> sol[0]
4.3027756377319948

* Octave 로 풀어보는 무리방정식
  (* Octave 로는 수치적 해를 구할 수 있다. fsolve 는 실수 범위내의 수치적 연산으로
     해를 구하고, fzero 는 복수수 범위 까지의 수치적 연산으로 해를 구해준다. 
     optimset 은 오차의 범위를 조절하기 위해 사용된다. 함수 f 의 구현에는 함수 f2 의
     구현과 달리 수식 부분 y = sqrt(x + 1) - x + 2 의 끝에 세미콜론이 없으므로 이 함수 f 가
     호출될 때 마다 이 세미콜론이 없는 구문을 만나면 y 에 저장되는 값이 출력된다, *)

octave-3.2.4.exe:1> function y = f(x)
> y = sqrt(x + 1) - x + 2
> endfunction
octave-3.2.4.exe:2> function y = f2(x)
> y = sqrt(x + 1) - x + 2;
> endfunction
octave-3.2.4.exe:3> opt = optimset("TolX", 1e-2)
opt =
{
  TolX =  0.010000
}

octave-3.2.4.exe:4> fsolve(@f, 0, opt)
y =  3
y =  3.0000
y =  3.0000
y =  1.7321
y =  1.7320
y =  1.7321
y =  0.24529
y = -0.0053636
y = 1.1321e-004
ans =  4.3026
octave-3.2.4.exe:5> fzero(@f, 0, opt)
y =  3
y =  2.4784
y =  2.3491
y =  3
y =  1.7321
y =  2.7247
y =  2.0811
y =  3
y =  1.7321
y =  12 +  3i
y = -4.6834
y =  0.31008
y = -2.1323
y = -0.0058277
y =  0.0051326
ans =  4.2962
octave-3.2.4.exe:6> fsolve(@f2, 0, opt)
ans =  4.3026
octave-3.2.4.exe:7> fzero(@f2, 0, opt)
ans =  4.2962
octave-3.2.4.exe:8> opt = optimset("TolX", 1e-12)
opt =
{
  TolX = 1.0000e-012
}

octave-3.2.4.exe:9> fsolve(@f, 0, opt)
y =  3
y =  3.0000
y =  3.0000
y =  1.7321
y =  1.7320
y =  1.7321
y =  0.24529
y = -0.0053636
y = 1.1321e-004
y = -2.3671e-006
y = 4.9493e-008
y = -1.0348e-009
ans =  4.3028
octave-3.2.4.exe:10> fzero(@f, 0, opt)
y =  3
y =  2.4784
y =  2.3491
y =  3
y =  1.7321
y =  2.7247
y =  2.0811
y =  3
y =  1.7321
y =  12 +  3i
y = -4.6834
y =  0.31008
y = -2.1323
y = -0.0058277
y = 3.3801e-006
y = -3.3794e-006
y = 4.4409e-016
y = -1.0973e-012
ans =  4.3028
octave-3.2.4.exe:11> fsolve(@f2, 0, opt)
ans =  4.3028
octave-3.2.4.exe:12> fzero(@f2, 0, opt)
ans =  4.3028

수치적 방법으로 근의 근사값을 구하는 방법은 여러 가지가 소개되어 있다.
다음은 방정식을 그 중에 가장 많이 쓰이는 (쉬운 방법) Newton 반복법을 사용하여
여러 가지 프로그램 언어로 구현해 보았다.

* 무리방정식을 풀기 위해 C 언어로 구현한 Newton 방법

/**
 * Filename solveSimpEqC_01.c
 *
 *  Purpose:
 *          Find roots of an equation by using the Newton method.
 *
 *  With VC++
 *      Compile: cl solveSimpEqC_01.c
 *      Execute: solveSimpEqC_01
 *
 *  With g++
 *      Compile: gcc solveSimpEqC_01.c -o solveSimpEqC_01
 *      Execute: ./solveSimpEqC_01
 *
 *  Date: 2011. 10. 30 (Sun)
 */

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

double f(double x) {
 return sqrt(x + 1) - x + 2;
}

double df(double x) {
 return 1/(2*sqrt(x + 1)) - 1;
}

double newtonIter(double (*func)(double), double (*dfunc)(double), double x0, double err) {
  double x = x0;
  double x1;
  int i;
  for (i = 0; i < 100; i++) {
   x1 = x;
   x = x - func(x)/dfunc(x);
   if (abs(x - x1) < err) {
       printf("Break loop after the %d-th iteration\n", i + 1);
       break;
   }
  }
  return x;
}

int main() {
 double x = newtonIter(&f, &df, 9, 1e-8);
    printf("x = %g\n", x);
    return 0;
}

/*
Execution result after compiling with VC++:
Break loop after the 2-th iteration
x = 4.30302

Execution result after compiling with gcc:
Break loop after the 4-th iteration
x = 4.30278
*/

* 무리방정식을 풀기 위해 C++ 언어로 구현한 Newton 방법

/**
 * Filename solveSimpEqCPP_01.cpp
 *
 *  Purpose:
 *          Find roots of an equation by using the Newton method.
 *
 *  With VC++
 *      Compile: cl solveSimpEqCPP_01.cpp /EHsc
 *      Execute: solveSimpEqCPP_01
 *
 *  With g++
 *      Compile: g++ solveSimpEqCPP_01.cpp -o solveSimpEqCPP_01
 *      Execute: ./solveSimpEqCPP_01
 *
 *  Date: 2011. 10. 30 (Sun)
 */

#include <iostream>
#include <cstring>
#include <cmath>

using namespace std;

double f(double x) {
 return sqrt(x + 1) - x + 2;
}

double df(double x) {
 return 1/(2*sqrt(x + 1)) - 1;
}

double newtonIter(double(*func)(double), double(*dfunc)(double), double x0, double err) {
  double x = x0;
  double x1;
  for (int i = 0; i < 100; i++) {
   x1 = x;
   x = x - func(x)/dfunc(x);
   if (abs(x - x1) < err) {
       cout << "Break loop after the " << i + 1 << "-th iteration" << endl;
       break;
   }
  }
  return x;
}

int main() {
 double sol = newtonIter(&f, &df, 9, 1e-8);
    cout << "x = " <<  sol << endl;
    return 0;
}

/*
Execution Result:
Break loop after the 4-th iteration
x = 4.30278
*/

* 무리방정식을 풀기 위해 Java 언어로 구현한 Newton 방법

/**
 *  Filename SolveSimpEqJava_01.java
 *
 *  Purpose:
 *          Find roots of an equation by using the Newton method.
 *
 *  Compile: javac SolveSimpEqJava_01.java
 *  Execute: java SolveSimpEqJava_01
 *
 *  Date: 2011. 10. 30 (Sun)
 */

import java.lang.*;

public class SolveSimpEqJava_01 {

 private static double f(double x) {
  return Math.sqrt(x + 1) - x + 2;
 }

 private static double df(double x) {
  return 1/(2*Math.sqrt(x + 1)) - 1;
 }

 private static double newtonIter(double x0, double err) {
  double x = x0;
  double x1;
  for (int i = 0; i < 100; i++) {
   x1 = x;
   x = x - f(x)/df(x);
   if (Math.abs(x - x1) < err) {
        System.out.printf("Break loop after the %d-th iteration\n", i + 1);
    break;
   }
   }
  return x;
 }

 public static void main(String[] args) {
  double sol = newtonIter(9, 1e-8);
  System.out.printf("x = %g\n", sol);
 }
}

/*
Execution Result:
Break loop after the 4-th iteration
x = 4.30278
*/

* 무리방정식을 풀기 위해 Groovy 언어로 구현한 Newton 방법
   (소스 내용은 Java 의 것과 같으나 파일 확장명만 groovy 로 바꾼 것)

/**
 *  Filename solveSimpEqGroovy_01.groovy
 *
 *  Purpose:
 *          Find roots of an equation by using the Newton method.
 *
 *  Execute: groovy solveSimpEqGroovy_01.groovy
 *
 *  Date: 2011. 10. 30 (Sun)
 */

import java.lang.*;

public class SolveSimpEqJava_01 {

 private static double f(double x) {
  return Math.sqrt(x + 1) - x + 2;
 }

 private static double df(double x) {
  return 1/(2*Math.sqrt(x + 1)) - 1;
 }

 private static double newtonIter(double x0, double err) {
  double x = x0;
  double x1;
  for (int i = 0; i < 100; i++) {
   x1 = x;
   x = x - f(x)/df(x);
   if (Math.abs(x - x1) < err) {
         System.out.printf("Break loop after the %d-th iteration\n", i + 1);
    break;
   }
   }
  return x;
 }

 public static void main(String[] args) {
  double sol = newtonIter(9, 1e-8);
  System.out.printf("x = %g\n", sol);
 }
}

/*
Execution Result:
Break loop after the 4-th iteration
x = 4.30278
*/

* 무리방정식을 풀기 위해 Groovy 언어로 구현한 Newton 방법
  (Groovy 언어의 특징을 살려 아주 간단하게 고친 것)

/**
 *  Filename solveSimpEq_01.groovy
 *
 *  Purpose:
 *          Find roots of an equation by using the Newton method.
 *
 *  Execute: groovy solveSimpEq_01.groovy
 *
 *  Date: 2011. 10. 30 (Sun)
 */

def f(x) {
 return Math.sqrt(x + 1) - x + 2
}

def df(x) {
 return 1/(2*Math.sqrt(x + 1)) - 1
}

def newtonIter(x0, err) {
 def x = x0
 def x1
 for (int i = 0; i < 100; i++) {
  x1 = x
  x = x - f(x)/df(x)
      if (Math.abs(x - x1) < err) {
           println("Break loop after the " + (i + 1) + "-th iteration")
           break
      }
  }
  return x
}

sol = newtonIter(9, 1e-8)
println("x = " +  sol)

/*
Execution Result:
Break loop after the 4-th iteration
x = 4.302775637731995
*/

* 무리방정식을 풀기 위해 Python 언어로 구현한 Newton 방법

#!/env/python
#
#  Filename solveSimpEq_01.py
#
#  Purpose:
#          Find roots of an equation by using the Newton method.
#
#  Execute: python solveSimpEq_01.py
#
#  Date: 2011. 10. 30 (Sun)

import math

def f(x):
 return math.sqrt(x + 1) - x + 2

def df(x):
 return 1/(2*math.sqrt(x + 1)) - 1

def newtonIter(x0, err):
 x = x0
 for i in range(0, 100):
   x1 = x
   x = x - f(x)/df(x)
   if abs(x - x1) < err:
       print "Break loop after the %d-th iteration" % (i + 1)
       break
 return x

sol = newtonIter(9, 1e-8)
print "x = %g" % sol

"""
Break loop after the 4-th iteration
x = 4.30278
"""

* 무리방정식을 풀기 위해 Ruby 언어로 구현한 Newton 방법

#!/env/ruby
#
#  Filename solveSimpEq_01.rb
#
#  Purpose:
#          Find roots of an equation by using the Newton method.
#
#  Execute: ruby solveSimpEq_01.rb
#
#  Date: 2011. 10. 30 (Sun)

def f(x)
 return Math::sqrt(x + 1) - x + 2
end

def df(x)
 return 1/(2*Math::sqrt(x + 1)) - 1
end

def newtonIter(x0, err)
 x = x0
 for i in 0..100 do
   x1 = x
   x = x - f(x)/df(x)
   if (x - x1).abs < err
       print "Break loop after the %d-th iteration\n" % (i + 1)
       break
      end
     end
 return x
end

sol = newtonIter(9, 1e-8)
print "x = %g\n" % sol

=begin
Execution Result:
Break loop after the 4-th iteration
x = 4.30278
=end

* 무리방정식을 풀기 위해 Lua 언어로 구현한 Newton 방법

--[[
  Filename solveSimpEq_01.lua

  Purpose:
          Find roots of an equation by using the Newton method.

  Execute: lua solveSimpEq_01.lua

  Date: 2011. 10. 30 (Sun)
]]--

function f(x)
 return math.sqrt(x + 1) - x + 2
end

function df(x)
 return 1/(2*math.sqrt(x + 1)) - 1
end

function newtonIter(x0, err)
 x = x0
 for i = 1,100 do
   x1 = x
   x = x - f(x)/df(x)
   if math.abs(x - x1) < err then
       print( "Break loop after the " .. (i + 1) .. "-th iteration" )
       break
      end
     end
 return x
end

sol = newtonIter(9, 1e-8)
print( "x = " .. sol )

--[[
Execution Result:
Break loop after the 5-th iteration
x = 4.302775637732
]]--