求助大神!!对已知数据点进行公式参数拟合!。

已知实验数据点(为了简化就先弄十个数据点,反射率关于温度的数据如下)。然后根据一定的理论模型,要对这些点进行拟合,拟合参数是两个未知数,但是我不会拟合,就先凑了两个数字上去(红色行的蓝色数字)。求教大神,如何编程实现拟合出来那两个参数??
temperature=[9689.7973
9746.4442
9766.11479
9636.48478
9489.68108
9732.24874
9490.63608
9332.22171
9399.91128
9173.54988];
reflectivity=[0.06542
0.06639
0.06534
0.06887
0.06933
0.06867
0.06816
0.06632
0.06874
0.06649];
n=0;
n0=1.39;
me=9*1E-31;
gamma=1.05;
Mol=18; %g/mol
for temp=5000:1:25000;
density=(0.03674*1E-3).*temp+2.57706; %根据(舒桦的2016、2017数据+2015Kimura数据)拟合density与temp的线性关系
pressure=density.*357.34-853; %GPa
Ni=3*density./Mol*6*1E29; %单位体积粒子数,粒子密度是h2o三个原子的密度
Eg=6.5-0.2*(density./1.18-1)-0.09*(temp./295-1);
Eg2kT = Eg.*1.6*(1E-19)/2/1.38/(1E-23)./temp; %Eg/2kT
yr = @(x) 2/sqrt(pi)*x.^(1/2)./(1+exp(x+Eg2kT)); %f(-Eg/2kT) fermi-dirac
it = integral(yr,0,inf); %fermi-dirac积分
ne = 2*(me*gamma*1.38*(1E-23).*temp./2/pi/(1.055*1E-34)^2)^(3/2).*it;
wp2 = ne.*1.6*1.6*(1E-38)/1.05/me/8.854*1E12; %wp^2
ve=((1.055*1E-34)*(3.*ne*pi*pi)^(1/3))/(me); %这里的me是否是电子质量?注意这里是ne不是Ni,是电子密度
l=2*(3/4/pi./Ni)^(1/3); %粒子间距
relaxtime=gamma.*l./ve; %弛豫时间计算公式
ns = sqrt(1-(wp2)*(660*1E-9)^2/(4*3.14*3.14*9*1E16)/(1+1i/(2*pi*3*1E8/660/1E-9*relaxtime))); %660nm
R = (abs((ns-n0)/(ns+n0)))^2 ;
n=n+1;
r(n)=R;
te(n)=temp;
hold on
end
plot(temperature,reflectivity,'*');
axis([5000 25000 0 0.5])
xlabel('Temperature K'),ylabel('Reflectivity')
hold on
plot(te,r);
axis([5000 25000 0 0.5])
xlabel('Temperature K'),ylabel('Reflectivity')
hold off

 Accepted Answer

查一下有没有退火算法simulannealbnd或者遗传算法ga的函数,有其中之一就可以代替粒子群算法
具体的options设置可能有版本差异,如果运行报错,就参考帮助文档修改,其他参数问题不大,主要是要设置HybridFcn结合fmincon来算.
退火
mysaopt = saoptimset('Display','iter','TemperatureFcn ',{@temperatureboltz},'PlotFcn',{@saplotbestf,@saplotf},'HybridFcn',{@fmincon});
[x,fval,exitflag,output] = simulannealbnd(@DeltaSquarePSO,[50;-50],[-100;-100],[100;100],mysaopt);
遗传
mygaopt = gaoptimset('Display','iter','PlotFcn',{@gaplotbestf},'HybridFcn',{@fmincon});
[x,fval,exitflag,output] = ga(@DeltaSquarePSO,2,[],[],[],[],[-100;-100],[100;100],[],[],mygaopt);

More Answers (0)

Categories

Find more on 数学 in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!