用gsl計算非方陣矩陣除法 解線性方程

2021-10-10 03:16:15 字數 1788 閱讀 6565

問題描述:

實際開發過程中遇到的乙個問題axispt=a\b(matlab表示式),求axispt.實際上就是求inv(a)*b,但是a不是乙個方陣,沒有找到非方陣如何求逆,遇到瓶頸。

解決辦法

以上等式等價於a*axispt = b,實際是解線性方程組的問題,並且係數矩陣不是方陣

實現過程

// gsl_matrix* matrixa = gsl_matrix_alloc(2*nrow+1, 3);

// gsl_matrix* matrixb = gsl_matrix_alloc(2*nrow+1, ncol);

// gsl_matrix* axispt = gsl_matrix_alloc(3, ncol);

gsl_vector* x =

gsl_vector_alloc(3

);gsl_vector* b =

gsl_vector_alloc(2

*nrow+1)

; gsl_vector* residual =

gsl_vector_calloc

(b->size)

; gsl_matrix* qr =

gsl_matrix_alloc

(matrixa-

>size1, matrixa-

>size2)

; gsl_vector* tau =

gsl_vector_calloc

(min

(matrixa-

>size1, matrixa-

>size2));

// 利用列迴圈把b的每一列依次當做方程等號右邊的結果,依次求得axispt的每一列值

for(

int n=

0; n)// 二者行數必須相等

if(matrixa-

>size1 != b-

>size)

// 複製係數矩陣

gsl_matrix_memcpy

(qr, matrixa)

;gsl_linalg_qr_decomp

(qr, tau)

;gsl_linalg_qr_lssolve

(qr, tau, b, x, residual)

;// 方程的解儲存在x中,將每一列得到的解x組合起來得到完整的axispt

for(

int k=

0; k<

3; k++)}

// 其他操作

// 物件釋放

gsl_vector_free

(tau)

;gsl_matrix_free

(qr)

;gsl_vector_free

(residual)

;gsl_vector_free

(b);

gsl_vector_free

(x);

誤區

①剛開始看到 「axispt=a\b」,一股腦地想著非方陣矩陣的逆該如何求,陷入了瓶頸。

②計算非線性方程求解,一開始用了lu分解(並且資料上有現成的例子),發現lu只能用來解係數矩陣是方陣的方程

③當時看了qr分解,被表示式a=qr勸退,心想我也沒有要把b分解成q和r啊,也沒接著細細往下看函式說明,實際上這個分解方法是把係數矩陣a分解成qr。正確的方法明晃晃擺在眼前,不靜下心深入研究,只知道惶惶不可終日覺得自己沒有能力解決這個問題,值得反省。

計算n階行列式和方陣逆矩陣

輸入n和乙個n階行列式,求結果 行列式就是化為上三角或下三角之後模擬手算 逆矩陣就是按這種方法做 1 2 3 1 0 0 兩邊做相同的初等行變換直到把左邊化為單位矩陣,右邊就是原矩陣的逆矩陣 2 5 6 0 1 0 4 1 3 0 0 1 寫完後只測試了幾組資料。include include in...

用Matlab計算jacobian矩陣解析解

做擴充套件卡爾曼濾波 ekf 的時候需要用到jacobian矩陣。有時手工求解難度較大這時可以用matlab自動求出jacobian矩陣的解析解。以雷達觀測矩陣為例為例 syms x y vx vy 定義符號變數 jacobian sqrt x2 y2 atan y x xvx yvy sqrt x...

非平衡電橋電阻計算 用非平衡電橋測量電阻

rrr rx2.非平衡電橋 非平衡電橋也稱不平衡電橋或微差電橋。圖 為非平衡電橋的原理圖,bd 之間為一 負載電阻rg 用非平衡電橋測量電阻時,是使rr 和r保持不變,rx 即r 變化時則 u變化。再根據u與 rx的函式關係,通過檢測 u的變化從而測得rx 由於可以檢測連續 變化的u 所以可以檢測連...