問題描述:
實際開發過程中遇到的乙個問題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 所以可以檢測連...