求助：曲线拟合的问题

ilxy

采用1stopt拟合曲线，编了一段代码，可能有语法错误，不能运行。我自己检查了半天也没检查出来，作过的朋友有时间的话帮着看一下。
下边是代码：

Parameter b1[1e-12,2e-8];
parameter b2[1e-3,1e3];
parameter b3[1e-11,1e-7];
parameter b4[1e-6,1e-3];
variable x,y;
Constant r0=6.18*0.0001;
Constant wt=28;
constant ce=100;
constant c=0.19;
constant sc=0.19;
constant hesum=109.6;
constant delta=0.5;
startprogram;
var
    i:integer;
    alafa:Double;
    rhand1:Double;
    kd:Double;
    de:Double;
    fenmu1:Double;
    fenmu2:Double;
    fenmu3:Double;
    rhand:Double;
    fenzi:Double;
    fenmu:Double;
Begin
     alafa:=0.0001;
    for i:=0:datalength-1 do begin
    x[i]:=(i+1)*delta;
    rhand1:=3/(0.25*3.2*r0);
    kd:=(b1/alafa^1.5)+b2*alafa^3;
    de:=b3*(ln(1/alafa))^2;
    fenmu1:=1/kd-r0/de;
    fenmu2:=(r0/de)*(1-alafa)^(-0.3333);
    fenmu3:=(1/b4)*(1-alafa)^(-0.6666);
    rhand:=rhand1/(fenmu1+fenmu2+fenmu3);
    fenzi:=rhand*ce*hesum;
    fenmu:=(1-alafa)*(wt+sc*ce)+alafa*c*(wt+ce);
    y[i]:=fenzi/fenmu;
    alafa:=rhand*delta+alafa;
    end;
End;
endprogram;
Data;
0.5  1.54
1   0.41
1.5  0.22
2   0.178
2.5  0.168
3   0.178
3.5  0.21
4   0.252
4.5  0.315
5   0.399
5.5  0.504
6   0.62
6.5  0.757
7  0.904
7.5  1.062
8   1.219
8.5  1.387
9   1.566
9.5  1.745
10  1.955
10.5 2.145
11  2.292
11.5  2.387
12  2.439
12.5 2.449
13  2.418
13.5   2.365
14  2.292
14.5  2.208
15  2.102
15.5  2.008
16  1.903
16.5  1.798
17  1.703
17.5 1.608
18  1.514
18.5  1.4299
19  1.345
19.5 1.272
20  1.209
20.5 1.146
21  1.093
21.5 1.041
22  0.9988
22.5  0.946
23  0.9147
23.5  0.873
24  0.841
24.5  0.8096
25  0.788
25.5  0.767
26  0.736
26.5 0.7255
27  0.704
27.5   0.683
28  0.672
28.5 0.662
29  0.641
29.5 0.63
30  0.62
30.5 0.609
31  0.609
31.5 0.599
32  0.588
32.5 0.578
33  0.568
33.5 0.558
34  0.548
34.5 0.536
35  0.535
35.5 0.526
36   0.515
36.5  0.504
37   0.494
37.5 0.484
38   0.473
38.5  0.4626
39   0.4416
39.5  0.4311
40    0.4310
40.5    0.41
41    0.399
41.5  0.389
42    0.3785
42.5  0.368
43    0.357
43.5  0.3469
44    0.336
44.5  0.325
45    0.315
45.5  0.305
46    0.294
46.5   0.284
47    0.283
47.5   0.273
48    0.262

[[i] 本帖最后由 bainhome 于 2007-1-26 01:20 编辑 [/i]]

shamohu

如果x为自变量，y为因变量，代码当中怎么又在计算x的值？
x:=(i+1)*delta;
拟合一般是由已知x值去计算y值，再与已知y值对比，再调整参数，代码中计算y值时，却好像并没用到x值，是否是以上原因？

ilxy

经过验证之后的正确代码如下：顺便提一句，在1stopt中，有关于时间序列的回归，可以用来确定偏微分方程中的待定系数。对于偏微分方程来说，将其写成差分形式，自变量统一视为时间t，在1stopt中通过差分迭代，即可得到相应的参数。这也是我自己琢磨出来回归的办法，以前一直在寻找这样的办法。1stopt的这个功能是其他的回归软件 matlab，origin等所不具备的，可以算是特色把。
Parameter b1[1e-12,2e-8];
parameter b2[1e-3,1e3];
parameter b3[1e-11,1e-7];
parameter b4[1e-6,1e-3];
variable x,y;
Constant r0=6.18*0.0001;
Constant wt=28;
constant ce=100;
constant c=0.19;
constant sc=0.19;
constant hesum=109.6;
startprogram;
var
    i:integer;
    alafa:Double;
    rhand1:Double;
    kd:Double;
    de:Double;
    fenmu1:Double;
    fenmu2:Double;
    fenmu3:Double;
    rhand:Double;
    fenzi:Double;
    fenmu:Double;
Begin
     alafa:=0.0001;
    for i:=0 to datalength-2 do begin
    rhand1:=3/(0.25*3.2*r0);
    kd:=(b1/alafa^1.5)+b2*alafa^3;
    de:=b3*(ln(1/alafa))^2;
    fenmu1:=1/kd-r0/de;
    fenmu2:=(r0/de)*(1-alafa)^(-0.3333);
    fenmu3:=(1/b4)*(1-alafa)^(-0.6666);
    rhand:=rhand1/(fenmu1+fenmu2+fenmu3);
    fenzi:=rhand*ce*hesum;
    fenmu:=(1-alafa)*(wt+sc*ce)+alafa*c*(wt+ce);
    y:=fenzi/fenmu;
    alafa:=rhand*(x[i+1]-x)+alafa;
    end;
End;
endprogram;
Data;

[ 本帖最后由 ilxy 于 2007-4-2 23:02 编辑 ]

shamohu

代码中有“for i:=0 to datalength-2 do begin”
如果DataLength=100，则最后第100个数据将没有计算值，而1stOpt是自动按DataLength来画图的，因此，ilxy 版主拟合时的图形上最后一点是否有偏差？严格来讲，按上述写法，最后一点是不该要的。

ilxy

Parameter b1 = 1e-12[1e-12,2e-8];
parameter b2 = 1e-3[1e-3,1e3];
parameter b3 = 1e-11[1e-11,1e-7];
parameter b4 = 1e-6[1e-6,1e-3];
Variable x, y;
Constant r0 = 6.18*0.0001, hesum=109.6;
StartProgram [Pascal];
Procedure MainModel;
var
    i: integer;
    alafa, rhand1, kd, de, fenmu1, fenmu2, fenmu3, rhand: Double;
Begin
    alafa := 0.0001;
    for i := 0 to datalength - 1 do begin
        rhand1 := 3/(0.25*3.2*r0);
        kd := (b1/alafa^1.5) + b2*alafa^3;
        de := b3*(ln(1/alafa))^2;
        fenmu1 := 1/kd-r0/de;
        fenmu2 := (r0/de)*(1 - alafa)^(-0.3333);
        fenmu3 := (1/b4)*(1 - alafa)^(-0.6666);
        rhand := rhand1/(fenmu1 + fenmu2 + fenmu3);
        y := hesum*rhand;
        alafa := rhand*0.1 + alafa;
    end;
End;
endprogram;
DataFile "CodeSheet1[A1:B485]";

[ 本帖最后由 ilxy 于 2007-4-3 12:20 编辑 ]