Chem & Coding | 一題雙解：用Matlab和R解方程組

這是一道來自北大化院陳慶德老師的《計算機在化學化工中的應用》課程中的作業題，原題目要求使用Matlab完成。

其Matlab代碼和結果如下：

function[y] = nl(x)y(1)=x(1)^2+2*x(2)^2-1;y(2)=2*x(1)^3-x(2);enda = [1 1 -1 -1;2 1 1 1;4 3 -1 -1];b = [0 0 0];c = [a b];rank(a)% ans = % 2rank(c)% ans = % 2%% rank(a)=rank(c)=2<4，方程有無窮多個解。rref(c)% ans = % 1 0 2 2 0 % 0 1 -3 -3 0 % 0 0 0 0 0%% 方程的通解為：%% x1+2x3+2x4=0%% x2-3x3-4x4=0x0=[0.8;0.6];fsolve(@nl, x0)% Equation solved.% fsolve completed because the vector of function values is near zero% as measured by the default value of the function tolerance, and% the problem appears regular as measured by the gradient.% <stopping criteria details>% ans = % 0.5567 % 0.5874%% 解得 x=0.5567 y=0.5874

使用R也可以完成類似的功能，代碼和運行結果如下：

library(pracma)a <- matrix(c(1,1,-1,-1,2,1,1,1,4,3,-1,-1), nrow = 3, ncol = 4, byrow = TRUE)b <- matrix(c(0,0,0), nrow = 3, ncol = 1, byrow = TRUE)c <- cbind(a,b)qr(a)$rank#> [1] 2qr(c)$rank#> [1] 2rref(c)#> [,1] [,2] [,3] [,4] [,5]#> [1,] 1 0 2 2 0#> [2,] 0 1 -3 -3 0#> [3,] 0 0 0 0 0library("BB")nl <- function(x){ y <- numeric(2) y[1]=x[1]^2+2*x[2]^2-1 y[2]=2*x[1]^3-x[2]^2 y}x <- c(0.8,0.6)result <- dfsane(x,nl) result$par#> [1] 0.5566931 0.5874065

推薦閱讀：