各位版上前輩好
小弟最近使用lsqcurvefit+網路上找的Newton法解聯立方程式
流程大概如下 如果有認知錯誤的地方請指正m(_ _)m
1. 一開始lsqcurvefit的起始猜測值傳入lsqcurvefit的副檔案裡面
2. 將lsqcurvefit的起始猜測值帶入聯立方程式的副檔案裡面進行計算
3. 計算之後的答案再用lsqcurvefit副檔案的方程式計算模擬值
4. lsqcurvefit傳出修正之後的猜測值
5. 再將修正之後的值傳入聯立方程式的檔案裡面進行計算
6. 得到解答
首先 主程式的程式碼如下
ans1 = lsqcurvefit(@GNFun,x0,xdata)
(x0是lsqcurvefit的起始猜測值 xdata用不到)
接著起始猜測值進入lsqcurvefit的副檔案裡面
function F = GNFun(x,xdata)
x0 = [一堆猜測值];
Q = NewtonRaphson(@Fun1,x0,1e-6);
F = 1.2*x(1)+1.5*Q(18)+2*Q(19)+1*Q(20);
(F才是lsqcurvefit所要擬合的式子)
Fun1裡面是一堆管流的方程式 並且呼叫Newton這個副檔案求解
function root = Newton(func,x,tol)
if size(x,1) == 1; % x must be column vector
x = x';
end
for i = 1:1e+5
[jac, f0] = jacobian(func, x);
if sqrt(dot(f0,f0)/length(x)) < tol
root = x; return
end
dx = jac\(-f0);
x = x + dx;
if sqrt(dot(dx,dx) / length(x)) < tol
root = x; return
end
end
error('請增加疊代步數或進行修改')
end
function [jac, f0] = jacobian(func,x)
% Returns the Jacobian matrix and f(x).
h = 1.0e-4;
n = length(x);
jac = zeros(n);
f0 = feval(func,x);
for i = 1:n
temp = x(i);
x(i) = temp + h;
f1 = feval(func,x);
x(i) = temp;
jac(:,i) = (f1 - f0)/h;
end
end
上面這個副檔案是小弟在網路上找到的 已經試過基本功能沒問題
可是我不知道要怎樣修改才能讓lsqcurvefit對於GNFun每一次的疊代值帶入Fun1
畫一張簡單的流程圖如下
http://i.imgur.com/shta1Qa.png
希望大大們不吝賜教或是給予修改的建議
謝謝QAQ