|
|
【答案】应助回帖
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ... jiqiang: 金币+50, ★有帮助, 还是很谢谢你~~~ 2013-05-31 23:11:14 csgt0: 金币+2, 谢谢 2013-06-03 14:39:36
对于第一组数据,按z0= 0.0015117处理
参数 Beta = 0.0000013251423784
各常数数值都很大,初值不好给出,MATLAB程序列如下,仅供参考,高版本的1stopt跑出来效果应该会好很多。
%--------------------------------------------------------------------------------------------------
function feixianxingnihe314
%clear all;clc
format long
zspan=[-40
-37.9375
-35.875
-33.8125
-31.75
-29.6875
-27.625
-25.5625
-23.5
-21.4375
-19.375
-17.3125
-15.25
-13.1875
-11.125
-9.0625
-7
-4.9375
-2.875
-0.8125
1.25
3.3125
5.375
7.4375
9.5
11.5625
13.625
15.6875
17.75
19.8125
21.875
23.9375
26
28.0625
30.125
32.1875
34.25
36.3125
38.375
40.4375
];
Texp=[
0.99345
0.99
0.98806
0.98477
0.98216
0.98184
0.97636
0.97683
0.96763
0.96138
0.9594
0.9472
0.92474
0.89896
0.8736
0.85005
0.80012
0.7539
0.74086
0.678
0.60586
0.68747
0.7964
0.85913
0.89399
0.90719
0.93361
0.95553
0.96915
0.97888
0.98453
0.98644
0.98746
0.98978
0.99194
0.99773
0.99861
1.00266
1.00016
1.00033
]; %E的数据,在此输入
k0=1e-6;
lb=0;
ub=inf;
%-------------------------------------------------------------------------
% 使用函数lsqnonlin()进行参数估计
OPTIONS=optimset('MaxFunEvals',1000);
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunc,k0,lb,ub,OPTIONS,Texp);
%ci = nlparci(k,residual,jacobian);
%residual;
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\n\t参数 D = %.16f',k(1))
y=KineticsEqs(k);
R2=1-sum((Texp-y').^2)./sum((Texp'-mean(y)).^2);
fprintf('\n\t相关系数之平方R^2 = %.16f',R2);
figure
plot(zspan,KineticsEqs(k),'b',zspan,Texp,'or'),legend('计算值','实验值','Location','Best')
%-------------------------------------------------------------------------
function f = ObjFunc(k,Texp)
f=KineticsEqs(k)'-Texp;
%------------------------------------------------------------------------
function xt = KineticsEqs(k)
zspan=[-40
-37.9375
-35.875
-33.8125
-31.75
-29.6875
-27.625
-25.5625
-23.5
-21.4375
-19.375
-17.3125
-15.25
-13.1875
-11.125
-9.0625
-7
-4.9375
-2.875
-0.8125
1.25
3.3125
5.375
7.4375
9.5
11.5625
13.625
15.6875
17.75
19.8125
21.875
23.9375
26
28.0625
30.125
32.1875
34.25
36.3125
38.375
40.4375
];
I0=1.492e12;
Leff=1;
z0= 0.0015117;
%z0=4.2764e-16;
%q0=k(1)*I0*Leff./(1+zspan.^2./(z0^2));
n=length(zspan);
for i=1:n;
q0=k(1)*I0*Leff./(1+zspan(i).^2./(z0^2));
fun=@(x) log(1+q0.*exp(-x.^2));
integ(i)=quad(fun,-1e5,1e5);
xt(i)=1./sqrt(pi)./q0.* integ(i);
end
![origin 非线性自定义拟合问题~~方程组编写及运行~~-1]()
附图1.jpg
|
|