무리방정식에는 무연근이 있을 수 있다.
여기서는 무리방정식
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 로도 무리방정식을 풀지 못해 이차반정식으로 고쳐서 풀었다. *)
>>> 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.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 에 저장되는 값이 출력된다, *)
> 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
]]--
'학습 > 수학' 카테고리의 다른 글
| 모든 양의 약수를 구하고 소수(prime number)인지 아닌지 판별하기 (0) | 2011.11.01 |
|---|---|
| 온라인 행렬 계산기 소개 (0) | 2011.10.27 |
| 직각삼각형의 빗변의 길이 구하기 / underflow 와 overflow 방지 (0) | 2011.10.22 |
| 여러가지 도구를 이용한 행렬의 고유값(eigenvalue)과 고유벡터(eigenvector) 구하기 (0) | 2011.10.20 |
| 여러 가지 도구로 연습해본 분수 계산 (0) | 2011.10.06 |
