使用R求解平面3D方程的线性方程组
solving set of linear equations using R for plane 3D equation
为了求解我的线性方程组,我遇到了一些麻烦。
我的示例中有三个 3D 点(A、B、C),我想自动求解我的系统。我想用这 3 个点创建一个平面。
手动(数学上)非常简单,但我不明白为什么我在编码时没有解决我的问题...
我有一个笛卡尔方程组,它是一个平面的方程:ax+by+cz+d=0
xAx + yAy + zA*z +d = 0 #point A
xBx + yBy + zB*z +d = 0 #point B
etc
我使用矩阵,例如A=(0,0,1)
; B=(4,2,3)
和 C=(-3,1,0)
。
通过手动求解,我有这个例子这个解决方案:x+3y-5z+5=0.
为了在 R 中解决它:我想使用 solve()
.
A <- c(0,0,1)
B <- c(4,2,3)
C <- c(-3,1,0)
res0 <- c(-d,-d,-d) #I don't know how having it so I tried c(0,0,0) cause each equation = 0. But I really don't know for that !
#' @param A vector 3x1 with the 3d coordinates of the point A
carteq <- function(A, B, C, res0) {
matrixtest0 <- matrix(c(A[1], A[2], A[3], B[1], B[2], B[3],C[1], C[2], C[3]), ncol=3) #I tried to add the 4th column for solving "d" but that doesn't work.
#checking the invertibility of my matrix
out <- tryCatch(determinant(matrixtest0)$modulus<threshold, error = function(e) e)#or out <- tryCatch(solve(X) %*% X, error = function(e) e)
abcd <- solve(matrixtest0, res0) #returns just 3 values
abcd <- qr.solve(matrixtest0, res0) #returns just 3 values
}
这不是好方法...但我不知道如何在我的问题中添加“d”。
我需要的return
是:return(a, b, c, d)
我觉得我的问题很经典而且简单,但是我没有找到像 solve() 或 qr.solve() 这样的函数可以解决我的问题...
你的解决方案实际上是错误的:
A <- c(0,0,1)
B <- c(4,2,3)
C <- c(-3,1,0)
CrossProduct3D <- function(x, y, i=1:3) {
#
To3D <- function(x) head(c(x, rep(0, 3)), 3)
x <- To3D(x)
y <- To3D(y)
Index3D <- function(i) (i - 1) %% 3 + 1
return (x[Index3D(i + 1)] * y[Index3D(i + 2)] -
x[Index3D(i + 2)] * y[Index3D(i + 1)])
}
N <- CrossProduct3D(A - B, C - B)
#[1] 4 2 -10
d <- -sum(N * B)
#[1] 10
#test it:
crossprod(A, N) + d
# [,1]
#[1,] 0
crossprod(B, N) + d
# [,1]
#[1,] 0
crossprod(C, N) + d
# [,1]
#[1,] 0
为了求解我的线性方程组,我遇到了一些麻烦。
我的示例中有三个 3D 点(A、B、C),我想自动求解我的系统。我想用这 3 个点创建一个平面。 手动(数学上)非常简单,但我不明白为什么我在编码时没有解决我的问题...
我有一个笛卡尔方程组,它是一个平面的方程:ax+by+cz+d=0
xAx + yAy + zA*z +d = 0 #point A
xBx + yBy + zB*z +d = 0 #point B
etc
我使用矩阵,例如A=(0,0,1)
; B=(4,2,3)
和 C=(-3,1,0)
。
通过手动求解,我有这个例子这个解决方案:x+3y-5z+5=0.
为了在 R 中解决它:我想使用 solve()
.
A <- c(0,0,1)
B <- c(4,2,3)
C <- c(-3,1,0)
res0 <- c(-d,-d,-d) #I don't know how having it so I tried c(0,0,0) cause each equation = 0. But I really don't know for that !
#' @param A vector 3x1 with the 3d coordinates of the point A
carteq <- function(A, B, C, res0) {
matrixtest0 <- matrix(c(A[1], A[2], A[3], B[1], B[2], B[3],C[1], C[2], C[3]), ncol=3) #I tried to add the 4th column for solving "d" but that doesn't work.
#checking the invertibility of my matrix
out <- tryCatch(determinant(matrixtest0)$modulus<threshold, error = function(e) e)#or out <- tryCatch(solve(X) %*% X, error = function(e) e)
abcd <- solve(matrixtest0, res0) #returns just 3 values
abcd <- qr.solve(matrixtest0, res0) #returns just 3 values
}
这不是好方法...但我不知道如何在我的问题中添加“d”。
我需要的return
是:return(a, b, c, d)
我觉得我的问题很经典而且简单,但是我没有找到像 solve() 或 qr.solve() 这样的函数可以解决我的问题...
你的解决方案实际上是错误的:
A <- c(0,0,1)
B <- c(4,2,3)
C <- c(-3,1,0)
CrossProduct3D <- function(x, y, i=1:3) {
#
To3D <- function(x) head(c(x, rep(0, 3)), 3)
x <- To3D(x)
y <- To3D(y)
Index3D <- function(i) (i - 1) %% 3 + 1
return (x[Index3D(i + 1)] * y[Index3D(i + 2)] -
x[Index3D(i + 2)] * y[Index3D(i + 1)])
}
N <- CrossProduct3D(A - B, C - B)
#[1] 4 2 -10
d <- -sum(N * B)
#[1] 10
#test it:
crossprod(A, N) + d
# [,1]
#[1,] 0
crossprod(B, N) + d
# [,1]
#[1,] 0
crossprod(C, N) + d
# [,1]
#[1,] 0