各位版上前辈好
小弟最近使用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