Re: [問題] 聯立方程式求解

看板MATLAB作者evanzxcv (鐵扶龍)時間17年前 (2007/04/21 05:18)推噓1(1推 0噓 0→)

留言1則, 1人參與討論串1/4 (看更多)

※ 引述《pto2 (蔡少)》之銘言： : 小弟我有小問題想請教各位先進，就是我有三個方程式分別為： : a1*x+b1*y+c1*z+d1 = 0 ..............(1) : a2*x^2+b2*y^2+c2*z^2+d2 = 0 ........(2) : a3*x^2+b3*y^2+c3*z^2+d3 = 0 ........(3) : 其中a1 a2 a3 b1 b2 b3 c1 c2 c3 d1 d2 d3 為常數 : 每次我用 solve 解就會當掉，請問有什麼好方法來求解這樣的問題嗎?? 我寫了個最佳化的程式solve1.m來解，需要猜一個起始值，測試程式的結果是，起始值假如和解差太遠(例如解是1,2,3，起始值猜一千之類的)，或方程組本身無實數解，就無法得解。此程式只能用來解以上方程式係數都為實數的實數解，假如係數有複數，或是需要複數解的話，可能要另外想辦法了... 先解釋程式中出現的代號好了，輸入的A,B是係數矩陣，X0是起始值，亦即： A=[a1 b1 c1 a2 b2 c2 a3 b3 c3] B=[d1 d2 d3] X0=[x起始值 y起始值 z起始值] 輸出的X是求得的近似解，sum是三個等式左邊的平方和，亦即 X=[x y z] sum=((1)式等號左方)^2+((2)式等號左方)^2+((3)式等號左方)^2 跑程式時把A,B,X0設好，然後直接用[X sum]=solve1(A,B,X0)來跑，我選擇把sum輸出的原因是檢驗程式有沒有跑成功，假如sum很接近零就沒錯，假如sum很大的話，換個X0試試看吧！至於程式裡面的myfun是solve1的副函數，他的目的是算出三個等式左方的平方和，因為係數和解都是實數，平方一定大於等於零，因此把三個平方相加，所得的數值也大於等於零，只要利用fminunc求最小值，那個最小值就是零了，就可以得到解。 myfun的輸出中fun是三個等式左方的平方和，而grad是fun的梯度，因為用梯度會讓fminunc跑得更快(我測試的結果是一眨眼就好了) %%%%%以下是程式solve1.m內容%%%%% function [X sum] = solve1(A,B,X0) options = optimset('GradObj','on','LargeScale','off','TolFun',eps,'TolX',eps); [X sum] = fminunc(@(X) myfun(X,A,B), X0, options); return function [fun grad] = myfun(X,A,B) sum1 = A(1,:)*X+B(1); sum2 = A(2,:)*(X.^2)+B(2); sum3 = A(3,:)*(X.^2)+B(3); fun = sum1^2 + sum2^2 + sum3^2; grad = zeros(3,1); for i = 1:3, grad(i) = 2*A(1,i)*sum1 + 4*A(2,i)*X(i)*sum2 + 4*A(3,i)*X(i)*sum3; end return %%%%%程式結束%%%%% 不知道這樣有沒有解決原po的問題，我是初學MATLAB的人啊，寫的程式可能也不是很好... -- ※ 發信站: 批踢踢實業坊(ptt.cc) ◆ From: 128.42.158.2 ※ 編輯: evanzxcv 來自: 128.42.158.26 (04/21 22:42)