推 pto2:感謝~~我試試看 04/22 14:05
※ 引述《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)