我正在解决数值分析中的一个问题（第 16 页是英文），这就是解决方案：
>> z=[0:500:4000 5000:1000:12000];
>> data=[5050 4980 4930 4890 4870 4865 4860 4860 4865 4875 4885 4905 4920 4935 4950 4970 4990];
>> fun=@(p)(4800 + p(1))*ones(size(z)) +p(2)/1000*z+p(3)*exp(p(4)/1000*z)-data;
>> x0=[0 0 0 -1];
>> opt = optimset('MaxFunEvals',1000);
>> p=lsqnonlin(fun,x0,[],[],opt);
Local minimum possible.
lsqnonlin stopped because the final change in the sum of squares relative to
its initial value is less than the default value of the function tolerance.
<stopping criteria details>
>> fitf=@(t)(4800 + p(1))*ones(size(t)) + p(2)/1000*t+ p(3)*exp(p(4)/1000*t);
>> tt=linspace(0,12000,1000);
>> plot(z,data,'r-',tt,fitf(tt),'b-');
>> p
p =
-20.2090 17.3368 272.9057 -0.7528
>>
在情节上看起来不错：
但是当我询问它时，答案中多项式的系数是其他一些值：
在我的程序运行中，coeff. 似乎是 p4=-0.7528，另一个是 p1=-20.2090，p2=17.3368，p3=272.9057。
所以 c(z) 变成
c(z)=4800-20.2090+(17.3368)z/1000+(272.9057)exp(0.7528z/1000)
和
c’(z)=17.3368/1000-0.2054exp(-z/1000)
然后我可以进行替换 z'=u 和 z''=u' 所以继续使用的 matlab 代码将是：
function dZ=sys(x,Z)
c=@(z)4800 -20.2090 + (17.3368)*z/1000+ (272.9057)*exp(0.7528z/1000); % c(z)
c=c(2000);
% Z(1):=z
% Z(2):=u
dZ=zeros(2,1); % a column vector
dZ(1)=Z(2);
dZ(2)=(c/cosd(7.8))^2*((1287*exp(Z(1)/1000))/5000 + 729/50000)/...
(4800 + 5.66 + (14.58)*Z(1)/1000+ (257.4)*exp(-Z(1)/1000))^3;
end
你同意吗？