matlab 函数拟合 线性插值_利用最小二乘法实现曲线拟合

1. 简介

在工程问题的计算中，我们经常需要处理一些离散数据的拟合问题，而最小二乘法是处理曲线拟合问题的常用方法。目前，许多软件都提供有基于最小二乘法进行曲线拟合的功能，例如在Origin和Excel中均可直接利用离散数据进行曲线拟合。然而，这些软件只能处理一些简单函数的拟合问题，当需要拟合的函数较为复杂时，或者无法用简单的表达式来表述时，则往往无法直接进行拟合。为此，本文将对最小二乘法的基本原理做简单介绍，随后介绍如何利用Matlab的lsqnonlin函数处理复杂函数的拟合问题。

1.1 曲线拟合的最小二乘法原理

利用最小二乘法进行曲线拟合的本质为寻找某个近似函数φ(x)，使得该函数与离散点之间尽可能地逼近。若将偏差定义为近似函数的近似值φ(x_i)与离散点y_i*之间的差值：

显然，偏差的大小可以用来反映近似函数的好坏。若该近似函数可使得偏差的平方和达到最小，即：

则这种求近似函数的方法即为离散数据曲线拟合的最小二乘法，函数φ(x)为最小二乘拟合函数，一般可选为一些简单的函数，例如选为多项式函数，则也称为最小多项式拟合。假定该多项式可表示为如下形式：

其中m为多项式阶次，n为离散数据对的个数。

根据最小二乘法的原理，则该问题实际上为选取参数a_j(j=0,1,…,m)使得如下所示的目标函数达到最小值：

根据多元函数求极值的必要条件，可得如下所示的方程组：

将上式移项并展开可得到如下形式的线性方程组：

求解上述线性方程组即可得到拟合多项式的系数。

2. 利用Matlab处理曲线拟合问题

基于上述计算原理，Matlab提供了polyfit函数用于处理多项式曲线的拟合问题，对于一些较为复杂但仍可通过简单表达式进行表述的函数，也可以利用Matlab的拟合工具箱(Curve fitting)进行拟合。但在某些情况，当拟合函数非常复杂，以致于无法用简单表达式进行表述时(例如分段函数以及涉及到条件语句)，则无法使用拟合工具箱进行拟合。对于此类问题，可以使用Matlab优化工具箱中的lsqnonlin函数进行解决。

2.1 lsqnonlin函数

lsqnonlin函数用于求解以下述形式表示的非线性最小二乘法拟合问题：

在使用该函数进行最小二乘法拟合时，lsqnonlin函数并不需要用户提供min ||f(x)||(平方和)，而是需要用户提供自定义函数fun，用于计算矢量形式表示的f(x)：

lsqnonlin函数常用语法为：

x = lsqnonlin(fun,x0)x = lsqnonlin(fun,x0,lb,ub)

其中fun为用户自定义函数，x0为计算采用的初始值，lsqnonlin函数首先利用x₀通过自定义函数fun计算f_i(x)的取值并计算平方和，随后通过优化算法调整x的取值直至得到平方和的最小值。此外，lb和ub还可以用于定义x的取值范围，使得x满足lb≤x≤ub。

例如，对于节1.1中所述的多项式，根据最小二乘法的定义，则自定义函数f(x)应表示为：

注意此时f(x)中的x为以向量形式表示的多项式P(x)的系数：

在计算时，用户需要指定多项式系数的初始值，则lsqnonlin函数将利用最小二乘法计算多项式系数。

下面，本文将以笔者所在领域常用的NASGRO方程为例，介绍如何利用lsqnonlin函数处理此类复杂函数的曲线拟合问题。

2.2 NASGRO方程简介

在进行基于断裂力学的损伤容限分析时，应力强度因子和裂纹扩展速率模型是最为重要的输入。一般来说，应力强度因子可以通过经验公式或数值方法进行计算，而裂纹扩展速率模型则需要通过裂纹扩展速率试验获得的试验数据拟合得到。例如，大量的试验结果表明，在裂纹扩展的中速率区域，应力强度因子幅值ΔK和裂纹扩展速率da/dN满足良好的对数线性关系，可以通过Paris公式进行描述：

其中C和m为材料常数。

尽管Paris公式已经得到广泛的应用，但是Paris公式仅仅描述了裂纹在中速率区域的扩展行为，没有描述近门槛区域和接近断裂的高速率区域的扩展行为，也没有考虑应力比R和裂纹闭合效应对裂纹扩展速率的影响，因此给出的计算结果将过于保守。另一个常用的裂纹扩展速率模型为Newman提出的NASGRO模型，该模型基于Forman模型改进了裂纹扩展速率模型，同时比Paris和Walker模型更加全面，不仅考虑了应力强度因子门槛值和断裂韧度，还体现了应力比以及裂纹闭合效应对裂纹扩展速率da/dN的影响，如图2.1所示，其表达式如下：

其中R为应力比，ΔK为应力强度因子幅值，ΔK_th为应力强度因子幅值门槛值，K_max为最大应力强度因子，可表示为：

K_c为材料断裂韧度，C, m, p, q为材料常数，f为Newman裂纹张开函数，可表示为：

其中，

其中，α为塑性约束因子，σ_max为最大应力水平，σ_F为流动应力，计算时可以取为单轴屈服应力和单轴强度极限的平均值。

图2.1 NASGRO方程

NASGRO方程中的应力强度因子门槛值ΔKth可采用下面的经验公式进行估算：

其中A₀为裂纹张开函数中的多项式系数，ΔK₁是R=1时的应力强度因子门槛值，C_th是对于正应力(上标为p=positive)和负应力比(上标为m=minus, negative)取不同值的材料常数，a₀是内在小裂纹尺寸(典型值为0.0381mm)。在基于NASGRO方程开发的疲劳裂纹扩展分析软件NASGRO中，正应力比下C_th^p和ΔK₁是保存在数据库里的值，负应力比下C_th^m的默认值为0.1。

2.3 NASGRO方程拟合

图2.2为疲劳裂纹扩展分析软件NASGRO材料库中某铝合金材料的裂纹扩展速率数据，已知试验时采用的试样为中心平板试样(M(T))，σ_max和σ_F的比值为0.3，塑性约束因子α为2.0，材料断裂韧度Kc为65.7，应力比R为1时的门槛值ΔK1为1.23，C_th^p为1.06，C_th^m为0.1，下面需要通过拟合试验数据获得NASGRO方程的参数C, m, p, q。

图2.2 疲劳裂纹扩展数据

拟合NASGRO方程的难点主要有以下几点：

(1)裂纹扩展速率da/dN不仅与应力强度因子幅值ΔK有关，还与使用的应力比R有关，因此实际上为多变量的拟合问题；

(2)裂纹张开函数f为分段函数，并且使用了计算最大值的max函数，该函数在拟合时无法用简单函数进行表述。

针对以上问题，NASGRO软件给出的拟合方法为首先给参数p和q确定一个初始值，并利用最小二乘法确定参数C和m，随后根据工程经验来获得可接受的结果，如果对拟合效果不满意，可以调整任意参数，直至获得满意的结果。

显然，这样的拟合策略具有很大的随意性，如果参数p和q选取不当，很可能对拟合效果有很大的影响。下面，本文将介绍如何利用lsqnonlin函数在不提前定义参数p和q的情况下对NASGRO方程进行拟合。

根据lsqnonlin函数的介绍，首先需要构造自定义函数f(x)使其满足最小二乘法计算的基本原理，由于Paris公式具有对数线性的关系，因此尝试将NASGRO方程两边取对数，可得：

上式可以用如下所示的通式表示：

系数b_j为与C, n, p和q有关(b₀=log(C), b₁=n, b₂=p, b₃=-q)的系数，而g_j为与ΔK、R和NASGRO中所有剩余参数有关的函数。

根据最小二乘法的定义，应选取参数 b_j 使得：

参考lsqnonlin函数对目标函数的定义，则自定义函数f(x)应表示为：

需要注意的是，y_i表示试验数据，应为裂纹扩展速率(da/dN)_i的对数值。

在定义函数 f( x) 时，lsqnonlin 函数明确指出该函数为向量形式表示的函数(vector-valued function )，意味着在计算时必须采用Matlab 特有的点乘计算，对于包含判断条件的向量化计算，可采用如下形式进行处理。例如，假设某一分段函数为如下形式：

则在Matlab中可采用如下语句进行向量化计算：

y= R .^2 * (R >0) + R * (R <= 0)

而裂纹张开函数f中涉及到求取最大值的计算以及分段函数的处理，也可以通过上述语法实现，具体的计算过程可参见程序代码(参见附录)。

此外，由于自定义函数f(x)为关于系数b_j(j=1,2,3,4)的函数，为了将试验数据(不同应力比R下的应力强度因子幅值ΔK和裂纹扩展速率da/dN)传递到函数f(x)中进行计算，可以将试验数据定义为全局变量，以便被f(x)调用。

通过编写程序，可以计算得到NASGRO方程的系数如表1所示。

拟合曲线与试验数据如图2.3所示。

图2.3 试验数据及拟合曲线对比

附录1 NASGRO方程曲线拟合程序

NASGRO_LSQ.m

NASGRO_LSQ用于定义采用最小二乘法拟合NASGRO方程时的自定义函数f(x)，输入参数Coeff为NASGRO方程系数b_j，输出参数为拟合函数与试验数据误差的平方和。

function F=NASGRO_LSQ(Coeff)%程序用于计算最小二乘法拟合NASGRO方程的目标函数%程序返回一个N×1的数值，其中N为数据对的个数%Coeff为拟合时待求的系数(共4个系数)%4个系数分别为log(C)、n、p和-q%DataX(:,1)为应力强度因子幅值%DataX(:,2)为应力比%DataY为裂纹扩展速率%*********全局变量传递**************global S_max_flow alpha DKth1 a0 Cth_p Cth_m a Kcglobal DataX DataY%S_max_flow为施加的最大应力与流动应力的比值%alpha为塑性约束因子%DKth1为应力比为1时对应的门槛值%a0为与门槛值有关的常数%Cth_p为正应力比下与门槛值有关的常数%Cth_m为负应力比下与门槛值有关的常数%a为计算门槛值时使用的裂纹长度，建议取为远大于a0的值%Kc为材料断裂韧度%DataX为应力强度因子幅值%DataY为裂纹扩展速率%***********************************%********Newman裂纹张开函数计算**********R=DataX(:,2);    %应力比DK=DataX(:,1);   %应力强度因子幅值%计算系数A0(与应力比和应力强度因子幅值无关)A0=(0.825-0.34*alpha+0.05*alpha^2)*...    (cos(pi/2*S_max_flow))^(1/alpha);A1=(0.415-0.071*alpha)*S_max_flow;A3=2*A0+A1-1;A2=1-A0-A1-A3;%计算向量形式的裂纹张开函数f1=max(A0+A1*R+A2*R.^2+A3*R.^3,R);f2=A0-2*A1;f3=A0+A1*R;f=f1.*(R>=0)+f2.*(R<-2)+...    f3.*(R>=-2&R<0);  %裂纹张开函数%****************************************%********应力强度因子门槛值计算**********DKth_p1=DKth1*sqrt(a/(a+a0))*((1-R)./(1-f)).^(1+R*Cth_p)./...        (1-A0).^((1-R)*Cth_p);        %正应力比下的门槛值DKth_p2=DKth1*sqrt(a/(a+a0))*((1-R)./(1-f)).^(1+R*Cth_m)./...        (1-A0).^(Cth_p-R*Cth_m);      %负应力比下的门槛值DKth=DKth_p1.*(R>=0)+...    DKth_p2.*(R<0);  %应力强度因子门槛值%****************************************%******根据NASGRO方程计算函数F***********F1=log10((1-f)./(1-R).*DK);   %DataX(1,:)为应力强度因子幅值F2=log10(1-DKth./DK);F3=log10(1-1./(1-R).*(DK./Kc));%****************************************%******根据NASGRO方程计算裂纹扩展速率***********y=Coeff(1)+Coeff(2)*F1+Coeff(3)*F2+Coeff(4)*F3;%***********************************************%*****构造基于最小二乘法的目标函数F**************%最小二乘法应保证目标函数F中所有原始之和达到最小F=(y-log10(DataY)).^2;%***********************************************end

NASGRO_nonfit.m

NASGRO_nonfit用于读取处理试验数据、定义材料参数以及调用lsqnonlin函数进行曲线拟合，并通过拟合得到的参数绘制NASGRO方程曲线。

%程序用于非线性拟合NASGRO方程的系数clear,clc%***************计算参数输入*****************global S_max_flow alpha DKth1 a0 Cth_p Cth_m a Kc  %定义全局变量global DataX DataY S_max_flow=0.3;         %最大应力与流动应力的比值alpha=2;                %塑性约束因子Kc=65.7;                %材料断裂韧度DKth1=1.23;            %应力比R=1对应的门槛值Cth_p=1.06;             %应力比为正时的门槛值计算参数Cth_m=0.1;              %应力比为负时的门槛值计算参数a0=0.0381;              %门槛值计算参数a=a0*1000;              %计算门槛值时采用的裂纹长度%************************************************%**************试验数据导入*****************%数据格式:应力强度因子,裂纹扩展速率,应力比FCG_data=importdata('FCG_data.txt');%数据处理data_num=size(FCG_data,1);         %数据对总个数FCG_data=sortrows(FCG_data,3);     %将数据按照应力比顺排  tbl=tabulate(FCG_data(:,3));       %数据统计R_eval=tbl(:,1);                        %试验出现的应力比(统计数据)R_num=size(R_eval,1);                   %试验使用的应力比数量DataX(:,1)=FCG_data(:,1);               %应力强度因子幅值DataX(:,2)=FCG_data(:,3);               %应力比DataY=FCG_data(:,2);                    %裂纹扩展速率clear tbl%********************************************%*************最小二乘法拟合*****************Coeff0=[-9 3 1 -1];      %拟合方程系数的迭代初始值[Coeff,resnorm,residual,exitflag,output]=lsqnonlin(@NASGRO_LSQ,Coeff0);%NASGRO方程系数输出fprintf('NASGRO方程系数C为%s\n',10^Coeff(1));fprintf('NASGRO方程系数m为%s\n',Coeff(2));fprintf('NASGRO方程系数p为%s\n',Coeff(3));fprintf('NASGRO方程系数q为%s\n',-1*Coeff(4));%********************************************%*************绘制试验数据及拟合曲线***************Output_DK=cell(R_num,1);       %用于输出数据的元胞数组初始化Output_DaDn=cell(R_num,1);%绘制试验数据figurefor i=1:R_num     %试样采用的应力比数量    data_index=find(FCG_data(:,3)==R_eval(i));    loglog(FCG_data(data_index,1),FCG_data(data_index,2),'o');    hold onend%绘制拟合曲线for i=1:R_num     %计算裂纹张开函数R=R_eval(i);     %应力比[f,A0,~]=Crack_opening_f(R,S_max_flow,alpha);%计算应力强度因子门槛值[DKth]=...    DK1th_calc(DKth1,R,f,A0,a0,Cth_p,Cth_m,a);%计算用于绘制NASGRO方程的应力强度因子幅值%应力强度因子幅值应大于门槛值并且最大应力强度因子应小于断裂韧度DK=linspace(DKth,(1-R)*Kc,500);   %应力强度因子幅值%计算裂纹扩展速率C=10^Coeff(1);m=Coeff(2);p=Coeff(3);q=-1*Coeff(4);Crack_rate=C*(((1-f)/(1-R))*DK).^m.*...    (1-(DKth./DK)).^p./...    (1-(1/(1-R)*(DK/Kc))).^q;loglog(DK,Crack_rate)%拟合曲线数据输出Output_DK{i}=DK';Output_DaDn{i}=Crack_rate';endlegend_str=cell(2*R_num,1);   %定义图例的元胞数组预分配for i=1:R_num         legend_str{i}=['R=' num2str(R_eval(i))];    %获取应力比用于图例显示    legend_str{i+R_num}=legend_str{i};endlegend(legend_str);title('NASGRO方程拟合');xlabel('应力强度因子幅值ΔK');ylabel('裂纹扩展速率da/dN');hold off%**************************************************

附录2 NASGRO方程绘制程序

程序仅用于计算绘制NASGRO方程时需要的参数，未参与自定义函数的计算。

Crack_opening_f.m

function [f,A0,A1,A2,A3] =...    Crack_opening_f(R,S_max_flow,alpha)%程序用于计算Newman裂纹张开函数%R为应力比%S_max_flow为最大应力与流动应力的比值S_max/S_flow%alpha为塑性约束因子A0=(0.825-0.34*alpha+0.05*alpha^2)*...    (cos(pi/2*S_max_flow))^(1/alpha);A1=(0.415-0.071*alpha)*S_max_flow;A3=2*A0+A1-1;A2=1-A0-A1-A3;if R>=0    f=A0+A1*R+A2*R^2+A3*R^3;    f=max([f;R]);elseif R    f=A0-2*A1;else    f=A0+A1*R;endend

DK1th_calc.m

function [DKth]=...    DK1th_calc(DKth1,R,f,A0,a0,Cth_p,Cth_m,a)%程序通过应力比为1的门槛值计算不同应力比下的门槛值%DKth1为应力比为1时对应的门槛值%R为应力比%f为裂纹张开函数,A0为裂纹张开函数中的多项式系数%a0为材料常数，Cth_p和Cth_m分别为正负应力比下的材料常数%a为裂纹长度if R>=0    DKth=DKth1*sqrt(a/(a+a0))*((1-R)/(1-f))^(1+R*Cth_p)/...        (1-A0)^((1-R)*Cth_p);else    DKth=DKth1*sqrt(a/(a+a0))*((1-R)/(1-f))^(1+R*Cth_m)/...        (1-A0)^(Cth_p-R*Cth_m);endend

试验数据

FCG_data.txt

应力强度因子幅值      裂纹扩展速率       应力比

11.9896  2.43E-05  0.711.6632  2.51E-05  0.710.5603  1.71E-05  0.710.916  1.43E-05  0.710.5604  1.43E-05  0.710.1043  1.10E-05  0.79.25035  8.18E-06  0.78.89987  5.91E-06  0.78.61001  5.10E-06  0.78.46855  6.09E-06  0.78.05821  4.95E-06  0.77.83889  4.40E-06  0.77.41798  3.91E-06  0.77.54187  3.48E-06  0.76.98103  3.38E-06  0.76.86652  1.99E-06  0.76.25152  1.99E-06  0.76.01454  2.44E-06  0.76.11514  1.17E-06  0.75.66032  1.62E-06  0.75.35634  1.82E-06  0.75.47599  1.24E-06  0.75.18191  1.39E-06  0.75.32693  1.24E-06  0.74.93077  1.44E-06  0.75.09686  1.07E-06  0.74.48922  9.23E-07  0.74.61486  8.20E-07  0.74.43997  7.73E-07  0.74.36708  6.48E-07  0.74.06469  6.87E-07  0.74.64051  4.41E-07  0.74.82333  4.04E-07  0.73.59998  4.04E-07  0.73.78335  3.49E-07  0.73.8678  3.70E-07  0.73.91075  3.39E-07  0.73.86786  2.38E-07  0.73.62002  1.72E-07  0.73.31402  1.93E-07  0.73.2238  2.05E-07  0.73.20604  2.38E-07  0.73.10161  2.24E-07  0.73.05068  2.05E-07  0.72.83944  2.31E-07  0.72.6722  1.82E-07  0.72.59947  1.57E-07  0.72.61388  1.24E-07  0.72.68703  1.17E-07  0.72.7929  9.53E-08  0.72.31499  1.28E-07  0.72.315  1.10E-07  0.72.34071  9.81E-08  0.72.26445  1.10E-07  0.72.20281  1.10E-07  0.72.00554  8.22E-08  0.71.96175  7.75E-08  0.71.92954  7.10E-08  0.71.88741  6.50E-08  0.71.86669  6.31E-08  0.71.82592  6.89E-08  0.71.74707  4.98E-08  0.71.70894  3.93E-08  0.71.68088  3.60E-08  0.71.64418  3.50E-08  0.71.63523  6.32E-09  0.71.58198  5.00E-09  0.743.3791  6.61E-04  039.9322  6.61E-04  035.9575  3.36E-04  033.2835  3.36E-04  031.4968  1.97E-04  031.3239  1.43E-04  028.8351  1.23E-04  028.6768  8.15E-05  026.6911  9.17E-05  026.5444  6.83E-05  024.9811  4.94E-05  023.7706  4.02E-05  023.2517  3.37E-05  022.6188  3.37E-05  021.7614  3.90E-05  020.8215  3.57E-05  020.8217  2.51E-05  021.4044  2.36E-05  021.7617  2.66E-05  019.9226  1.81E-05  019.062  1.71E-05  018.6457  1.98E-05  019.2735  2.16E-05  019.5951  2.51E-05  017.7421  1.76E-05  017.9392  1.52E-05  017.0698  1.61E-05  016.6051  1.61E-05  016.3325  1.47E-05  016.5139  1.27E-05  015.888  1.24E-05  015.7137  1.10E-05  015.6275  6.65E-06  014.465  1.16E-05  014.4652  9.20E-06  014.3065  6.65E-06  013.8403  9.76E-06  013.8402  1.13E-05  013.2425  8.18E-06  013.3159  7.48E-06  013.2427  6.09E-06  012.5315  7.27E-06  012.3939  7.71E-06  012.1233  6.85E-06  011.9902  6.27E-06  012.5316  5.10E-06  011.7285  4.95E-06  010.8563  5.91E-06  010.7371  5.41E-06  010.7966  4.95E-06  010.6193  4.67E-06  010.3875  4.15E-06  09.99386  3.28E-06  09.40513  3.28E-06  09.30192  2.91E-06  09.25075  2.67E-06  08.99896  2.51E-06  08.90013  2.75E-06  08.19305  1.93E-06  07.97006  1.71E-06  07.79604  1.62E-06  07.796  1.87E-06  07.54204  1.87E-06  06.9813  1.24E-06  07.05883  1.07E-06  06.53398  7.96E-07  06.829  7.08E-07  06.67993  5.76E-07  06.49812  5.27E-07  06.32121  5.93E-07  06.08171  4.29E-07  06.32129  4.16E-07  06.32131  3.81E-07  06.35634  3.29E-07  05.85128  3.10E-07  05.66068  3.01E-07  05.75522  2.67E-07  05.26876  2.67E-07  05.26878  2.38E-07  05.41621  2.38E-07  05.35678  2.05E-07  05.35682  1.72E-07  05.29803  1.57E-07  05.21103  1.53E-07  05.15379  1.72E-07  05.21097  1.99E-07  05.09716  2.24E-07  05.06913  1.99E-07  04.95844  1.99E-07  04.82348  1.82E-07  04.71814  1.93E-07  04.74422  2.31E-07  04.58968  2.38E-07  04.61513  1.77E-07  04.9312  1.48E-07  04.6407  1.57E-07  04.66641  1.36E-07  04.48956  1.28E-07  04.39151  1.40E-07  04.29561  1.40E-07  04.31936  1.72E-07  04.29556  1.99E-07  04.20179  1.62E-07  04.20181  1.40E-07  04.24848  1.21E-07  04.15571  1.21E-07  04.11007  1.36E-07  04.06492  1.57E-07  04.04256  1.40E-07  03.97618  1.40E-07  03.9325  1.67E-07  03.84666  1.32E-07  03.88939  1.04E-07  03.88942  8.47E-08  04.04264  8.47E-08  03.70087  1.32E-07  03.7009  1.10E-07  03.60016  1.10E-07  03.62006  1.32E-07  03.52153  1.24E-07  03.52155  1.04E-07  03.38806  1.21E-07  03.50218  8.98E-08  03.6002  8.22E-08  03.5022  7.75E-08  03.42573  7.53E-08  03.36947  7.53E-08  03.31411  9.25E-08  03.2239  9.25E-08  03.18854  7.31E-08  03.25972  6.89E-08  03.31417  5.77E-08  03.36953  4.84E-08  03.08471  5.13E-08  03.08473  4.43E-08  03.1536  4.17E-08  03.31424  3.40E-08  03.20635  1.83E-08  03.01752  1.14E-08  03.01755  9.28E-09  03.17125  8.49E-09  02.96804  6.14E-09  03.08497  5.79E-09  03.10206  4.85E-09  03.05116  3.40E-09  02.87152  1.37E-09  02.79336  1.37E-09  021.6401  2.50E-04  0.521.8805  1.97E-04  0.520.8201  1.97E-04  0.519.9211  1.23E-04  0.520.7061  9.17E-05  0.517.9382  6.83E-05  0.517.9383  5.56E-05  0.517.6438  4.94E-05  0.516.6044  5.09E-05  0.516.6965  3.90E-05  0.515.8873  4.02E-05  0.515.6265  3.37E-05  0.515.0343  3.27E-05  0.514.5447  2.51E-05  0.514.3059  2.29E-05  0.513.0967  2.10E-05  0.512.0561  1.71E-05  0.512.3253  1.47E-05  0.512.4623  1.13E-05  0.511.7281  1.13E-05  0.511.1597  1.24E-05  0.511.0983  1.13E-05  0.510.5023  1.24E-05  0.510.5024  9.76E-06  0.510.33  9.20E-06  0.510.4447  8.42E-06  0.510.2732  7.71E-06  0.59.93858  6.65E-06  0.59.50927  6.65E-06  0.58.89981  7.06E-06  0.58.99867  5.74E-06  0.58.8022  5.25E-06  0.58.28374  4.27E-06  0.58.32966  3.38E-06  0.57.92594  3.91E-06  0.57.71025  3.28E-06  0.56.90432  4.40E-06  0.57.21618  2.51E-06  0.56.75371  2.37E-06  0.56.67961  1.99E-06  0.56.21708  2.23E-06  0.56.04791  1.77E-06  0.56.04792  1.66E-06  0.55.38609  1.10E-06  0.55.29767  8.96E-07  0.55.23953  8.20E-07  0.55.12513  7.73E-07  0.55.01323  7.08E-07  0.55.26863  4.97E-07  0.54.9038  5.93E-07  0.54.66611  7.29E-07  0.54.41561  4.68E-07  0.54.31924  3.59E-07  0.54.20166  3.59E-07  0.53.99806  3.19E-07  0.53.66015  2.45E-07  0.53.50207  2.18E-07  0.53.48282  1.77E-07  0.53.18843  1.77E-07  0.53.03396  1.17E-07  0.53.01723  1.53E-07  0.52.83949  1.44E-07  0.52.80833  1.32E-07  0.52.65752  1.36E-07  0.52.45989  1.48E-07  0.52.54274  1.14E-07  0.52.59952  9.53E-08  0.52.44638  1.04E-07  0.52.39295  1.14E-07  0.52.3798  9.25E-08  0.52.2896  8.98E-08  0.52.22728  9.53E-08  0.52.21503  8.47E-08  0.52.16666  8.47E-08  0.52.13107  9.25E-08  0.52.11936  7.53E-08  0.52.09609  7.10E-08  0.52.05032  7.10E-08  0.52.07309  6.31E-08  0.52.0731  5.77E-08  0.52.05034  5.77E-08  0.52.01667  5.77E-08  0.51.99453  5.28E-08  0.51.99454  4.84E-08  0.51.99454  4.56E-08  0.51.98357  4.17E-08  0.51.92958  3.82E-08  0.51.98358  3.71E-08  0.51.98361  2.68E-08  0.51.9727  2.38E-08  0.52.0279  2.18E-08  0.52.02791  2.06E-08  0.51.98367  1.14E-08  0.51.86688  4.71E-09  0.51.84641  3.21E-09  0.51.83625  3.12E-09  0.5