查看: 1986|回复: 16|关注: 0

[已解决] 求助大神!!对已知数据点进行公式参数拟合!

[复制链接]

新手

12 麦片

财富积分


050


1

主题

8

帖子

0

最佳答案
本帖最后由 何至于 于 2018-6-9 21:57 编辑

已知实验数据点(为了简化就先弄十个数据点,反射率关于温度的数据如下)。然后根据一定的理论模型,要对这些点进行拟合,拟合参数是两个未知数,但是我不会拟合,就先凑了两个数字上去(红色行的蓝色数字)。求教大神,如何编程实现拟合出来那两个参数??

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

论坛优秀回答者

中级

837 麦片

财富积分


5001500


0

主题

2126

帖子

182

最佳答案
  • 关注者: 111
发表于 2018-6-9 22:55:05 | 显示全部楼层
图简单,方便,又不想麻烦学Matlab编程的话,试试1stOpt。

新手

12 麦片

财富积分


050


1

主题

8

帖子

0

最佳答案
 楼主| 发表于 2018-6-9 23:24:22 | 显示全部楼层
shihe 发表于 2018-6-9 22:55
图简单,方便,又不想麻烦学Matlab编程的话,试试1stOpt。

我的拟合公式里有积分,有e指数,所以还是想用matlab编程实现。但是公式太复杂,试过直接套用拟合公式不行,还是需要自己编,对这方面是初学者,想请教大神给我指示和思路~谢谢!!!

论坛优秀回答者

中级

729 麦片

财富积分


5001500


2

主题

572

帖子

126

最佳答案
  • 关注者: 39
发表于 2018-6-10 09:21:37 | 显示全部楼层
你的参数有没有取值范围上的限制?
Eg=6.5-a*(density./1.18-1)-b*(temp./295-1);
a,b ∈ [-Inf,Inf]?
如果没有限制可能不太好处理.

新手

12 麦片

财富积分


050


1

主题

8

帖子

0

最佳答案
 楼主| 发表于 2018-6-10 11:06:06 | 显示全部楼层
TouAkira 发表于 2018-6-10 09:21
你的参数有没有取值范围上的限制?
Eg=6.5-a*(density./1.18-1)-b*(temp./295-1);
a,b ∈ [-Inf,Inf]?

虽然没什么限制,但是可以人为加个-100到100的限制先~

论坛优秀回答者

中级

729 麦片

财富积分


5001500


2

主题

572

帖子

126

最佳答案
  • 关注者: 39
发表于 2018-6-10 13:25:36 | 显示全部楼层
何至于 发表于 2018-6-9 23:06
虽然没什么限制,但是可以人为加个-100到100的限制先~

在你给的取值范围内 a = 7.80555723021025; b = -0.223225262406654; 使残差平方和最小.

  1. clear;
  2. myopt = optimoptions('particleswarm','Display','iter','PlotFcn','pswplotbestf','SwarmSize',100,'HybridFcn',@fmincon);
  3. [x,fval,exitflag,output] = particleswarm(@DeltaSquarePSO,2,[-100;-100],[100;100],myopt);

  4. a = x(1);
  5. b = x(2);
  6. n0=1.39;
  7. me=9*1E-31;
  8. gamma=1.05;
  9. Mol=18;
  10. temp=(5000:25:25000)';
  11. density=(0.03674*1E-3).*temp+2.57706;
  12. pressure=density.*357.34-853;
  13. Ni=3*density./Mol*6*1E29;
  14. Eg=6.5-a*(density./1.18-1)-b*(temp./295-1);
  15. Eg2kT = Eg.*1.6*(1E-19)/2/1.38/(1E-23)./temp;
  16. it = integral(@(x) 2/sqrt(pi)*x.^(1/2)./(1+exp(x+Eg2kT)),0,inf,'ArrayValued',true);
  17. ne = 2*(me*gamma*1.38*(1E-23).*temp./2/pi/(1.055*1E-34)^2).^(3/2).*it;
  18. wp2 = ne.*1.6*1.6*(1E-38)/1.05/me/8.854*1E12;
  19. ve=((1.055*1E-34)*(3.*ne*pi*pi).^(1/3))./(me);
  20. l=2*(3/4/pi./Ni).^(1/3);
  21. relaxtime=gamma.*l./ve;
  22. 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)));
  23. R = (abs((ns-n0)./(ns+n0))).^2;
  24. temperature=[9689.7973;9746.4442;9766.11479;9636.48478;9489.68108;9732.24874;9490.63608;9332.22171;9399.91128;9173.54988];
  25. reflectivity=[0.06542;0.06639;0.06534;0.06887;0.06933;0.06867;0.06816;0.06632;0.06874;0.06649];
  26. figure;

  27. plot(temperature,reflectivity,'*',temp,R,'k-','MarkerSize',12);
  28. axis([5000 25000 0 0.5]);
  29. xlabel('Temperature K'),ylabel('Reflectivity');
复制代码
把下面代码存储为DeltaSquarePS.m
  1. function [Out] = DeltaSquarePSO (In)
  2. a = In(1);
  3. b = In(2);

  4. syms x;
  5. n0=1.39;
  6. me=9*1E-31;
  7. gamma=1.05;
  8. Mol=18;
  9. temp=[9689.7973;9746.4442;9766.11479;9636.48478;9489.68108;9732.24874;9490.63608;9332.22171;9399.91128;9173.54988];
  10. reflectivity=[0.06542;0.06639;0.06534;0.06887;0.06933;0.06867;0.06816;0.06632;0.06874;0.06649];

  11. density=(0.03674*1E-3).*temp+2.57706;
  12. pressure=density.*357.34-853;
  13. Ni=3*density./Mol*6*1E29;
  14. Eg=6.5-a*(density./1.18-1)-b*(temp./295-1);
  15. Eg2kT = Eg.*1.6*(1E-19)/2/1.38/(1E-23)./temp;
  16. it = integral(@(x) 2/sqrt(pi)*x.^(1/2)./(1+exp(x+Eg2kT)),0,inf,'ArrayValued',true);
  17. ne = 2*(me*gamma*1.38*(1E-23).*temp./2/pi/(1.055*1E-34)^2).^(3/2).*it;
  18. wp2 = ne.*1.6*1.6*(1E-38)/1.05/me/8.854*1E12;
  19. ve=((1.055*1E-34)*(3.*ne*pi*pi).^(1/3))./(me);
  20. l=2*(3/4/pi./Ni).^(1/3);
  21. relaxtime=gamma.*l./ve;
  22. 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)));
  23. R = (abs((ns-n0)./(ns+n0))).^2 ;

  24. Out = sum( 10^6*(R - reflectivity).^2 )';
  25. end
复制代码

新手

12 麦片

财富积分


050


1

主题

8

帖子

0

最佳答案
 楼主| 发表于 2018-6-11 14:52:29 | 显示全部楼层
TouAkira 发表于 2018-6-10 13:25
在你给的取值范围内 a = 7.80555723021025; b = -0.223225262406654; 使残差平方和最小.

把下面代码存储 ...

万分感谢大神的回复!但在运行时一直有以下的报错,不知是什么原因?

错误使用 optimoptions (line 114)
Invalid solver specified. Provide a solver name or handle (such as 'fmincon' or @fminunc).
Type DOC OPTIMOPTIONS for a list of solvers.

论坛优秀回答者

中级

729 麦片

财富积分


5001500


2

主题

572

帖子

126

最佳答案
  • 关注者: 39
发表于 2018-6-11 19:24:20 | 显示全部楼层
何至于 发表于 2018-6-11 02:52
万分感谢大神的回复!但在运行时一直有以下的报错,不知是什么原因?

错误使用 optimoptions (line 114) ...

你用的哪个版本?比较早的版本可能没有particleswarm这个函数
先输入doc particleswarm看看帮助文档里面有没有这个函数

新手

12 麦片

财富积分


050


1

主题

8

帖子

0

最佳答案
 楼主| 发表于 2018-6-12 09:19:22 | 显示全部楼层
TouAkira 发表于 2018-6-11 19:24
你用的哪个版本?比较早的版本可能没有particleswarm这个函数
先输入doc particleswarm看看帮助文档里面有 ...

我用的2014版的,查了下没有这个函数~需要自己写一个吗

论坛优秀回答者

中级

729 麦片

财富积分


5001500


2

主题

572

帖子

126

最佳答案
  • 关注者: 39
发表于 2018-6-12 12:32:09 | 显示全部楼层 |此回复为最佳答案
本帖最后由 TouAkira 于 2018-6-12 00:35 编辑
何至于 发表于 2018-6-11 21:19
我用的2014版的,查了下没有这个函数~需要自己写一个吗

查一下有没有退火算法simulannealbnd或者遗传算法ga的函数,有其中之一就可以代替粒子群算法
具体的options设置可能有版本差异,如果运行报错,就参考帮助文档修改,其他参数问题不大,主要是要设置HybridFcn结合fmincon来算.

退火
  1. mysaopt = saoptimset('Display','iter','TemperatureFcn ',{@temperatureboltz},'PlotFcn',{@saplotbestf,@saplotf},'HybridFcn',{@fmincon});
  2. [x,fval,exitflag,output] = simulannealbnd(@DeltaSquarePSO,[50;-50],[-100;-100],[100;100],mysaopt);
复制代码

遗传
  1. mygaopt = gaoptimset('Display','iter','PlotFcn',{@gaplotbestf},'HybridFcn',{@fmincon});
  2. [x,fval,exitflag,output] = ga(@DeltaSquarePSO,2,[],[],[],[],[-100;-100],[100;100],[],[],mygaopt);
复制代码
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /3 下一条

快速回复 返回顶部 返回列表