고유방정식 repeat(),uniroot.all(),뉴튼랩슨법으로 고유값 구하기

https://pastryofjsmath.tistory.com/81

LA14) 고유값,고유벡터(eigen value,eigen vector)-2

위에 대한 내용의 문제풀이입니다.

 

f <- function(x) {-x^3 +4*x^2 +3*x -14}
plot(f,xlim = c(-5,5))
abline(h=0,v=0,col='black')
abline(v=-3,col='red')

> A <- matrix(c(0,0,7,
+               4,2,5,
+               1,0,2),byrow=T,ncol=3)
> eigen(A)$values
[1]  3.828427  2.000000 -1.828427

위 고유값과 실제로 같은지 구해보자.

> #초기값지정
> x <- -3
> repeat{
+   if(-x^3 +4*x^2 +3*x -14<0) break
+   x <- x + 0.00001
+ }
> x
[1] -1.82842
> repeat{
+   if(-x^3 +4*x^2 +3*x -14>0) break
+   x <- x + 0.00001
+ }
> x
[1] 2
> 
> repeat{
+   if(-x^3 +4*x^2 +3*x -14<0) break
+   x <- x + 0.00001
+ }
> x
[1] 3.82843

> #패키지이용
> library(rootSolve)
> f <- function(x) {-x^3 +4*x^2 +3*x -14}
> uniroot.all(f,interval = c(-5,5))
[1]  2.000000 -1.828428  3.828430

library(plotly)
library(magrittr)

# Newton-Raphson 
nt <- function(fun, x0 = -100, N = 100, tol = 1e-6) {
  h <- 1e-6
  iter_data <- data.frame(step = integer(), x = numeric(), y = numeric())
  i <- 1
  
  while (i <= N) {
    df.dx <- (fun(x0 + h) - fun(x0)) / h
    x1 <- x0 - fun(x0) / df.dx
    
    # 현재 x와 f(x) 저장
    iter_data <- rbind(iter_data, data.frame(step = i, x = x1, y = fun(x1)))
    
    if (abs(x0 - x1) < tol) break
    x0 <- x1
    i <- i + 1
  }
  
  return(iter_data)  # x 값을 저장한 데이터프레임 반환
}

# 함수 정의
fun <- function(x) -x^3 + 4*x^2 + 3*x - 14

# Newton-Raphson 결과
results <- nt(fun, x0 = -10, N = 1000, tol = 1e-20)

# 곡선이 그려질 범위
x_vals <- seq(-10, 10, length.out = 1000)
curve_data <- data.frame(x = x_vals, y = fun(x_vals))

# 애니메이션 생성
fig <- plot_ly() %>%
  # 함수 추가
  add_lines(data = curve_data, x = ~x, y = ~y, name = "Function Curve", line = list(color = 'blue')) %>%
  
  # 점 추가 (애니메이션)
  add_markers(data = results, x = ~x, y = ~y, frame = ~step, name = "Newton-Raphson",
              marker = list(size = 10, color = 'red')) %>%
  
  # 레이아웃 설정
  layout(title = "Newton-Raphson Method",
         xaxis = list(title = "x"),
         yaxis = list(title = "f(x)"))

fig

초기값을 조정해가면서 근을 수치해석적으로 구할 수 있습니다.

library(pracma)
library(matlib)
A <- matrix(c(0,0,7,
              4,2,5,
              1,0,2),byrow=T,ncol=3)

> D1 <- diag(eigen(A)$values[1],3,3)
> gaussianElimination(A-D1,c(0,0,0))
     [,1] [,2]      [,3] [,4]
[1,]    1    0 -1.828427    0
[2,]    0    1 -6.734591    0
[3,]    0    0  0.000000    0
> 
> D2 <- diag(eigen(A)$values[2],3,3)
> gaussianElimination(A-D2,c(0,0,0))
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    0    1    0
[3,]    0    0    0    0
> 
> D3 <- diag(eigen(A)$values[3],3,3)
> gaussianElimination(A-D3,c(0,0,0))
     [,1] [,2]      [,3] [,4]
[1,]    1    0  3.828427    0
[2,]    0    1 -2.693981    0
[3,]    0    0  0.000000    0

행렬 A의 고유값에 대응하는 고유벡터는 위와 같이 구할 수 있고

A^5의 고유값에 대응하는 고유벡터는 A의 고유값에 대응하는 고유벡터와 동일합니다.

> D1 <- diag(eigen(B)$values[1],3,3)
> gaussianElimination(B-D1,c(0,0,0))
     [,1] [,2]      [,3] [,4]
[1,]    1    0 -1.828427    0
[2,]    0    1 -6.734591    0
[3,]    0    0  0.000000    0
> 
> D2 <- diag(eigen(B)$values[2],3,3)
> gaussianElimination(B-D2,c(0,0,0))
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    0    1    0
[3,]    0    0    0    0
> 
> D3 <- diag(eigen(B)$values[3],3,3)
> gaussianElimination(B-D3,c(0,0,0))
     [,1] [,2]      [,3] [,4]
[1,]    1    0  3.828427    0
[2,]    0    1 -2.693981    0
[3,]    0    0  0.000000    0

#A^5의 고유값 = (A의 고유값)^5
> (eigen(A)$values)^5 ; eigen(B)$values
[1] 822.43564  32.00000 -20.43564
[1] 822.43564  32.00000 -20.43564