首页：http://www.7d-soft.com/index.htm
4.0新功能 (预定2010年8月6日)：
1：支持复数拟合、复数方程组计算；
2：支持微分方程拟合求解；
3：通用全局优化求解器变异功能，优化能力提高20%以上；
4：新的编程模式计算引擎；
5：强大易用的数据批处理拟合功能
6：公式自动搜索：增加更多的二维、三维函数库；
7：改进的积分计算，拟合，解方程可含有积分函数，支持高斯积分和辛普森积分算法
8：三维图形旋转、缩放、移动等功能
9：？号输入，可动态输入常数。
10：参数定义更加方便*：Parameter 0<=a<=10, b=[1,3];
11：更多的数学函数支持：Wrap、Wrap0…
12：支持更多功能的关键字：FileWeight，OutWeight…
13：重复计算时自动记录每次结果
14：Exp函数计算修正，与Matlab等保持一致：Exp(-3^0.23)-> Exp(-(3^0.23))
15：….
3.0新功能 (2009年5月1日)：
New in 3.0
1：重新设计的与其它高级语言的接口，更加方便与C++, Fortran, Basic, Pascal等语言的浑编联动。
2：增加新的算法：稳健全局优化算法。
3：改进了离子群和最大继承算法，优化能力更强。
4：增加了常微分方程（ODE）的求解功能，算法包括：龙格-库塔-费尔博格法(Runge-Kutta-Fehlberg Method)、欧拉算法（Euler
Method）、2-5阶龙格-库塔算法（Runge-Kutta Method），不仅能求解一般的初值ODE方程，还能解特殊形式的ODE方程，对边值问题的ODE方程也能轻松求解。
5：对线性规划问题自动判断识别，速度更快。
6：更加灵活的LoopConstant定义：LoopConstant d=[2,(max(x,1))];
7：与Vista兼容
8：编程模式增加对特殊函数的支持（Erf, Erfc, Gamma, Bessel…）
9：Parameter对For的支持。
10：拟合计算结束进行预测时，可计算每一点的导数
11：SubDivision、RunNext与Inherit功能
12：LogFile自动保存功能
13：RowData、RowDataSet与EndRowDataSet关键字
14：更加方便的Sum()，Prod()和For()语句
15：3D图形新格式：点状三维图
16：“恢复刚关闭的文件“功能
….
2.5新功能 (2006年10月7日)：
1：更加出色、稳健的通用全局优化能力
1：对Basic的全面支持
2：参数型变量问题的拟合（未知中间变量）：ParVariable
3：带积分的拟合和函数优化
4：隐函数优化算法的改进，速度增加10倍
5：隐函数拟合算法的订正：TradImplicit， ImplicitRange
6：BatchFile： 文件批处理功能
7：StepReg：逐步拟合功能
8：CodeSheet：代码本表格，支持直接从表格中读取数据
9：代码本显示形式：单业、多业和下拉
10：LoopConstant、FullLoopModel：自动循环计算功能
11：Constant a(1:3)=[1,2,3] -> Constant a = [1,2,3]
12：WeightReg：灵活多变、任意形式的带权重拟合
13：PassParameter：编程模式下支持返回计算变量
14：参数初值自动选择更加智能、健壮，适应范围更广
15：RegType：最小二乘法、最小一乘法等不同形式拟合
16：MDataSet，EndMDataset：网络节点数据自动转至矩阵数据
17：HotRun：设定自动热计算及计算次数
18：Sum，Prod，For更简洁写法
19：编程模式下可以直接定义二维参变量
20：Plot、PlotLoopData：迭代计算过程中更加丰富、强大的动态图形表示方式
21：众多改进及Bug修正
2.0新功能 (2006年10月7日)：
1：求解非线性方程组功能大幅改进，【麦夸特法+通用全局优化算法】已成为解非线性方程组的首选算法，其改进后的求解能力总体上强于任何其它算法。
2：最大最小优化问题求解 (MinMax)：一种多目标优化求解功能。
3：智能拟合功能：该功能特别适合于数据量很大时的拟合，可数倍甚至数十倍缩短计算时间，数据量越大，效果越明显。
4：改进的非常容易实现的带等式或不等式约束的拟合
5：算法自动选择功能：对于刚接触1stOpt的用户而言，由于不清楚各算法的特点及适用范围，常无法确定如何选择算法，该功能可根据问题的类型自动选择算法。
6：函数表达式以脚本语言描述表达功能：对于复杂、繁琐、冗长的问题，可通过脚本语言来描述
7：常字符串数组表达功能：定义字符串数组功能
例：ConstStr S(1:3) = [x1^2+x2, x1*x2-x2^2, sin(x1)+x2];
等同于：ConstStr S1 = x1^2+x2, S2 = x1*x2-x2^2, S3 = sin(x1)+x2;
例：ConstStr S(1:3) = x2*[x1^2+x2, x1*x2-x2^2, sin(x1)+x2];
等同于：ConstStr S1 = x2*(x1^2+x2), S2 = x2*(x1*x2-x2^2), S3 = x2*(sin(x1)+x2);
8：公式拟合自动搜索时稳健模式搜索功能
9：0-1规划，修正数值范围溢出问题
10：公式自动拟合库中，增加众多峰函数
11：约束函数连续表达功能：
例：10.3>=x1+sin(x2)*x3>=0
等同于：
x1+sin(x2)*x3>=0;
x1+sin(x2)*x3<=10.3;
例：Parameter x1[0.5,0.66], x4[0.04,0.2], x7[,0.035];
MinFunction 0.44*x1+0.94*x2+0.88*x3+0.48*x4+4*x5+3.4*x6+2.3*x7+0.12*x8+1.6*x9+19*x10+25*x11;
3230*x1+2640*x2+2500*x3+1730*x4+2900*x5+2230*x6+2500*x7>2750;
8.27*x1+43*x2+40*x3+15.4*x4+62*x5+50*x6+45*x7>15;
8.27*x1+43*x2+40*x3+15.4*x4+62*x5+50*x6+45*x7<16;
0.038*x1+0.32*x2+0.32*x3+0.14*x4+3.91*x5+4.6*x6+33.4*x8+21*x9>2.85;
0.038*x1+0.32*x2+0.32*x3+0.14*x4+3.91*x5+4.6*x6+33.4*x8+21*x9<3;
0.058*x1+0.15*x2+0.14*x3+0.32*x4+2.9*x5+2.15*x6+0.14*x8+18.5*x9>0.5;
0.058*x1+0.15*x2+0.14*x3+0.32*x4+2.9*x5+2.15*x6+0.14*x8+18.5*x9<0.55;
0.26*x1+2.45*x2+2.41*x3+0.54*x4+4.35*x5+3.28*x6+2.6*x7+99*x11>0.8;
0.125*x1+0.48*x2+0.51*x3+0.18*x4+1.65*x5+1.31*x6+0.65*x7+99*x10>0.31;
0.298*x1+1.08*x2+1.4*x3+0.58*x4+2.21*x5+1.74*x6+0.83*x7+99*x10>0.58;
0.298*x1+1.08*x2+1.4*x3+0.58*x4+2.21*x5+1.74*x6+0.83*x7+99*x10<0.63;
0.077*x1+0.6*x2+0.6*x3+0.27*x4+0.8*x5+0.64*x6>0.19;
x2+x3>0.1;
x2+x3<0.22;
x5+x6>0.03;
x5+x6<0.07;
x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11=1;
可写为：
Parameter x1[0.5,0.66], x4[0.04,0.2],x7[,0.035];
MinFunction 0.44*x1+0.94*x2+0.88*x3+0.48*x4+4*x5+3.4*x6+2.3*x7+0.12*x8+1.6*x9+19*x10+25*x11;
3230*x1+2640*x2+2500*x3+1730*x4+2900*x5+2230*x6+2500*x7>2750;
16>8.27*x1+43*x2+40*x3+15.4*x4+62*x5+50*x6+45*x7>15;
3>0.038*x1+0.32*x2+0.32*x3+0.14*x4+3.91*x5+4.6*x6+33.4*x8+21*x9>2.85;
0.55>0.058*x1+0.15*x2+0.14*x3+0.32*x4+2.9*x5+2.15*x6+0.14*x8+18.5*x9>0.5;
0.26*x1+2.45*x2+2.41*x3+0.54*x4+4.35*x5+3.28*x6+2.6*x7+99*x11>0.8;
0.125*x1+0.48*x2+0.51*x3+0.18*x4+1.65*x5+1.31*x6+0.65*x7+99*x10>0.31;
0.63>0.298*x1+1.08*x2+1.4*x3+0.58*x4+2.21*x5+1.74*x6+0.83*x7+99*x10>0.58;
0.077*x1+0.6*x2+0.6*x3+0.27*x4+0.8*x5+0.64*x6>0.19;
0.22>x2+x3>0.1;
0.07>x5+x6>0.03;
x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11=1;
12：矩阵计算，基本函数求导计算
13：带权重的拟合功能
14：带约束的超越方程求解
15：For语句，支持循环表达式
16：支持自动重复计算
17：改进的预测/验证功能
18：DataSet，AutoData定义数据时，可指定起始基数：
缺省时，起始基数为1
例：
DataSet;
EndDataSet:
例：AutoData x = 1:1:10;
例：定义起始基数为0
DataSet [Base = 0];
EndDataSet:
例：AutoData[Base = 0] x = 1:1:10;
19：增加IFF关键字
20：代码中直接从Excel表单和1stOpt电子表格中读取数据：必须指定文件名、表单名及数据范围
例：从Excel文件“C:\Data1.xls”中的“Sheet1”中读取数据进行拟合计算，数据范围从A1到B20
Function y = a + b*x + Exp(c*x);
DataFile C:\Data1.xls[Sheet1[A1:B20]];
21：常数连续定义：
例：Constant A(1:3) = 2;
等同于 Constant A1 = 2, A2 = 2, A3 = 2;
例：Constant A(1:3) = 10*[1,2,3];
等同于 Constant A(1:3) = [10,20,30];
22：增强的编程模式，可完善自动处理任意多的等式及不等式约束，对于复杂的带约束的工程问题，可轻易求解。
23：修正定义多维常数、参数时出现的问题
1.5新功能 (2006年4月18日)：
1：单纯形线性规划算法中，可进行整数规划、混合整数规划计算。
2：编程模式中，对约束条件的自动处理功能。
3: 权重拟合功能
4：结果数据自动保存功能。
5：同一代码本中，所有问题同时求解功能。
6：函数优化预测检验功能
7：数据自动产生功能: 关键字： AutoData
例：AutoData X = 1:1:10, Y = X^2+X；
等同于：Constant X(1:10) = [1,2,3,4,5,6,7,8,9,10];
Constant Y(1:10) = [2,6,12,20,30,42,56,72,90,110];
8: 循环语句关键字：For，支持无穷镶套
例：For(i=1:3)(x[i]>=A[i]*i);
等同于： x[1] >=A[1]*1;
x[2] >=A[2]*2;
x[3] >=A[3]*3;
9: 新增特殊密度分布函数：BetaCDF, BetaPDF, BinoCDF, BinoPDF, Chi2CDF, Chi2PDF, ExpCDF,
ExpPDF, PoissCDF, PoissPDF, TCDF, TPDF
10：增加函数求导计算功能
例：(x*exp(x+sin(x)))’ ==>
diff(x*exp(x+sin(x)),x) = exp(x+sin(x))+exp(x+sin(x))*(1+cos(x))*x
diff(x*exp(x+sin(x)),x=3) = exp(x+sin(x))+exp(x+sin(x))*(1+cos(x))*x = 23.82417126
11：新增：
BinParameter: 定义0-1变量；
IntParameter: 定义正整数变量；
ParameterDomain：定义变量范围；
PlotXYZData：画三维数据图；
PlotMeshData：画三维网格数据图；
PlotPoint3D：画三维点图；
12:众多改进，运行更快捷、稳定。
错误纠正：
1：函数表达式中出现空格显错的问题。
2：拟合时，用“DataFile”调用外部数据文件出错。
3：用超过两次“DataSet- EndDataSet”定义数据时出错
4：拟合时，用“SkipStep“出错。
5：函数中出现诸如“2E+10“时显错的问题。
6: 其它诸多Bugs
//There are many test examples located in the folder: examples\ //Nonlinear regression: any No. of variables and parameters, any user-defined function
Title "NLReg Demo - 1";
Parameter b1, b2, b3, b4;
Variable x, y;
Function y=exp(b1*x)+exp(b2*x)+b3*sin(b4*x);
Data;
//x y 2.1 4.1
2.2 4.3
2.3 4.5
2.4 4.6
2.5 4.8
2.6 5.0
2.7 5.1
2.8 5.3
2.9 5.5 NewDivision "NLReg Demo - 2";
//Data may be in codesheet, the simplified formation is like below:
Function y=exp(b1*x)+exp(b2*x)+b3*sin(b4*x);
DataFile "CodeSheet1[A2:B11]"; //Press F2 to view data NewDivision "NLReg Demo - 3";
//The data may be in row type:
Function y=exp(b1*x)+exp(b2*x)+b3*sin(b4*x);
RowData;
,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9;
,4.1,4.3,4.5,4.6,4.8,5.0,5.1,5.3,5.5; NewDivision "NLReg Demo - 4";
//In program model (Pascal):
Variable x, y;
Parameter b();
StartProgram [Pascal];
Procedure MainModel;
var i: integer;
Begin
for i := to DataLength - do
y[i] := exp(b1*x[i])+exp(b2*x[i])+b3*sin(b4*x[i]);
End;
EndProgram;
RowData;
,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9;
,4.1,4.3,4.5,4.6,4.8,5.0,5.1,5.3,5.5; NewDivision "NLReg Demo - 5";
//In program model (Basic):
Variable x, y;
Parameter b();
StartProgram [Basic];
Sub MainModel
Dim i as integer
for i = to DataLength -
y(i) = exp(b1*x(i))+exp(b2*x(i))+b3*sin(b4*x(i))
Next
ConstrainedResult = b1+b2+b3=
End Sub
EndProgram;
DataFile "CodeSheet1[A2:B11]"; //Press F2 to view data NewDivision "NLReg Demo - 6";
//Constrained regression: b1, b2, b3 >=0; b1+b2+b3=1;
Parameter b() =[,];
Function y=exp(b1*x)+exp(b2*x)+b3*sin(b4*x);
b1+b2+b3=;
DataFile "CodeSheet1[A2:B11]"; //Press F2 to view data NewDivision "NLReg Demo - 7";
//parameters are shared by various function
SharedModel;
Variable x,y1,y2;
Function y1=a+g*t^*x^/(+t^*x^)+c*x^d;
y2=b+g*t*x/(+t^*x^)+c*x^d1;
DataFile "CodeSheet1[C2:E23]"; //Press F2 to view data NewDivision "NLReg Demo 8";
//curve fitting with integration
Variables x, y;
Parameters a, b, c, d;
Function y=a-b*exp(-c*x^d)*int((t+x)*sin(x),t=0.1:c);
data;
0.05 0.13
0.15 0.13
0.25 0.19
0.35 0.34
0.45 0.53
0.55 0.71
0.65 1.06
0.75 1.6
0.85 1.64
0.95 1.83
1.05 2.09
1.15 2.05
1.25 2.13
1.35 2.12
1.45 2.09 NewDivision "NLReg Demo 9";
{Two times of x-y data with different T, k1, k2, k3, k4 are parameters,
regression equation: y=k1*k3*exp((k2+k4)/T)*x/(+k3*exp(k4/T)*x) Data：
T= 288K T=303K
x y x y 3.41586 0.04887 8.31345 0.0901
9.47655 0.13519 28.43464 0.26128
24.02448 0.29356 59.3231 0.4764
56.42759 0.51347 94.48793 0.66564
93.07483 0.7103 135.639 0.836
132.8176 0.8863 179.1531 0.97904
172.8579 1.03949 218.6459 1.0981
210.2345 1.16252 254.779 1.18808
238.6986 1.27173 283.5697 1.25296
272.3641 1.35695 308.7593 1.30915
297.1748 1.42801
313.8014 1.48757 How to estimate paremeters of k？ }
Variable x,y;
Parameter k();
VarConstant t=[,];
StartProgram [Pascal];
Procedure MainModel;
var j: integer;
Begin
for j := to DataLength - do
y[j] := k1*k3*exp((k2+k4)/T)*x[j]/(+k3*exp(k4/T)*x[j]); End;
EndProgram;
Data;
//x y 3.41586 0.04887
9.47655 0.13519
24.02448 0.29356
56.42759 0.51347
93.07483 0.7103
132.8176 0.8863
172.8579 1.03949
210.2345 1.16252
238.6986 1.27173
272.3641 1.35695
297.1748 1.42801
313.8014 1.48757
Data;
//x y 8.31345 0.0901
28.43464 0.26128
59.3231 0.4764
94.48793 0.66564
135.639 0.836
179.1531 0.97904
218.6459 1.0981
254.779 1.18808
283.5697 1.25296
308.7593 1.30915 NewDivision "NLReg Demo 10";
Function y=c*(+Erf(((x-b)*exp(-a*(x/b)))/(sqrt()*a*b)));
Data; 0.04
0.08
0.14104
0.18635
0.24
0.3
0.38272
0.44474
0.5573
0.64532
0.69053
0.76975
0.82394
0.86369
0.90259
0.93255
0.96256
0.98994
0.99023
0.99621
0.99821 NewDivision "NLReg Demo 11";
Parameter a,b,c,d;
Function y = (A-D)/(+(x/C)^B) + D;
Data;
0.01953125 0.624
0.078125 0.728
0.15625 0.827
0.3125 0.976
0.390625 1.138
0.625 1.243
2.5 1.292
1.414
1.496 NewDivision "Equation Solve - 1";
Function x*0.090949*0.274666*( - y*Exp(-*z)) - a = 0.00097;
x*0.206949*0.274666*( - y*Exp(-*z)) - a = 0.00681;
x*0.230149*0.274666*( - y*Exp(-*z)) - a = 0.00785;
x*0.241749*0.274666*( - y*Exp(-*z)) - a = 0.01131; NewDivision "Equation Solve - 2";
Constant R=,D=,bta=0.75,p=10e6;
Parameter x,y,m,n;
Function
p+*bta*R^*x/(*R+*D-m)^+*bta*R^*x/(*R+*D-m)^+*bta*R^*y/(R+D-n)^=R^*y/(R-n)^;
p+*bta*R^*x/(*R+*D-m)^+*bta*R^*x/(R+D-m)^+*bta*R^*y/(*R+D-m)^=R^*y/(R+n)^;
p+*bta*R^*x/(*R+*D-m)^+*bta*R^*y/(*R+*D-n)^+*bta*R^*y/(R+D+n)^=R^*x/(R-m)^;
p+*bta*R^*x/(*R+*D-m)^+*bta*R^*y/(*R+*D-n)^+*bta*R^*y/(*R+D+n)^=R^*x/(R+m)^; NewDivision "Equation Solve - 3";
{System equations:
-*cos(a1)+*cos(a2)-*cos(a3)+*cos(a4)-*cos(a5)+*cos(a6)=1.60;
-*cos(*a1)+*cos(*a2)-*cos(*a3)+*cos(*a4)-*cos(*a5)+*cos(*a6)=;
-*cos(*a1)+*cos(*a2)-*cos(*a3)+*cos(*a4)-*cos(*a5)+*cos(*a6)=;
-*cos(*a1)+*cos(*a2)-*cos(*a3)+*cos(*a4)-*cos(*a5)+*cos(*a6)=;
-*cos(*a1)+*cos(*a2)-*cos(*a3)+*cos(*a4)-*cos(*a5)+*cos(*a6)=;
-*cos(*a1)+*cos(*a2)-*cos(*a3)+*cos(*a4)-*cos(*a5)+*cos(*a6)=;}
ParameterDomain = [,pi/];
Constant p=[,,,,,], v=[1.6,,,,,];
Function For(i=:,p,v)(-*cos(p*a1)+*cos(p*a2)-*cos(p*a3)+*cos(p*a4)-*cos(p*a5)+*cos(p*a6)=v); NewDivision "ODE Function";
Variable t = [,], y=, z=;
Plot y[x], z;
ODEOptions = [SN=,A=,P=];
ODEFunction y' = cos(z'-y')*t*y/(1+z);
z' = (t*z+y*exp(y)*sin(y'))/; NewDivision "Function Optimization - 1";
MinFunction Ba;
Ba*SO4=;
BaOH/Ba*OH=4.8;
HSO4/SO4*H=0.98;
H*OH=;
Ba+^(-)*BaOH=SO4+^(-)*HSO4;
*Ba+1E-*BaOH+1E-2H=*SO4+1E-*HSO4+1E-*OH; NewDivision "Rosenbrok function";
Constant n=;
Parameters x(:n+)[-1.5,1.5];
Minimum = true;
Function sum(i=:n)(*(x[i+]-x[i]^)^+(-x[i])^); NewDivision "Ackley function 2 - Optimization";
Constant n = ;
Parameters x(:n)[-,];
Minimum = true;
Function sum(i=:n-)(*(cos(*x[i])+sin(*x[i+]))+sqrt(x[i+]^+x[i]^)); NewDivision "NLReg Demo 12";
//nonlinear regression with multi-files
Variable t,h;
Parameter a, a1, k, k1;
VarConstant b=[,,,];
StartProgram [Pascal];
Procedure MainModel;
var i: integer;
temb, h0: double;
Begin
temb := b;
h0 := ;
for i := to DataLength - do begin
h0 := h0+ a*(temb-a1*tanh(k1*h0))*exp(k*(temb-a1*tanh(k1*h0)))*;
h[i] := h0;
end;
End;
EndProgram;
Data;
//0 0
0.226066238
0.281510409
0.314326458
0.351110333
0.375454391
0.398296935
0.417922218
0.435617147
0.451381719
Data;
//0 0
0.5653775
0.8946725
1.1468305
1.3735105
1.5691755
1.7503255
1.9100265
2.0534495
2.1794185
Data;
//0 0
1.760498
2.784063
3.657308
4.438827
5.171771
5.860787
6.513054
7.140139
7.746008
Data;
//0 0
4.03011
6.880921
9.285056
11.364176
13.288876
14.892276
16.389606
17.784386
19.036266 NewDivision "ODE Regression";
//Eurle method
//x'=dx/dt=p1*(20-x)+p2*(p3-x);
//y'=dy/dt=p1*(x-y)+p2*(p3-y)
//z'=dz/dt=p1*(y-z)+p2*(p3-z)
Constant V = [,,]; //Initial Values
Constant dt=;
Parameter p();//=[0.5,100];
Variable t, x[OutPut], y[OutPut], z[OutPut];
StartProgram [Pascal];
Procedure MainModel;
var i: integer;
x0,y0,z0: double;
Begin
x0 := V[];
y0 := V[];
z0 := V[];
for i := to DataLength - do begin
x0 := (p1*(-x0)+p2*(p3-x0))*dt+x0;
y0 := (p1*(x0-y0)+p2*(p3-y0))*dt+y0;
z0 := (p1*(y0-z0)+p2*(p3-z0))*dt+z0;
x[i] := x0;
y[i] := y0;
z[i] := z0;
end;
End;
EndProgram;
Data;
//0 10 15 20
11.80 15.37 20.77
14.04 17.04 22.23
15.44 16.85 21.16
17.80 18.07 22.49
18.60 19.43 24.46
19.22 20.12 23.58
21.77 21.83 24.51
22.17 22.11 25.38
23.41 23.37 25.53
23.17 24.92 26.69
24.56 25.55 29.21
25.85 26.06 28.00
24.64 28.94 30.32
25.15 28.94 30.41
26.92 30.06 31.87
26.37 29.29 31.87
26.71 31.48 34.60
26.61 31.86 33.57
27.13 33.84 36.75
29.32 32.95 36.36
28.24 33.03 36.23
28.42 32.50 36.01
28.11 33.12 39.19
28.98 35.32 37.29
30.23 35.56 37.79
30.21 34.86 42.05
29.11 35.40 42.67
30.42 36.20 40.74
28.84 35.91 40.53
29.44 36.50 43.33 NewDivision "ODE Regression";
//4th-order RK
Constant y1_1=, y2_1=, y3_1=; //Initial Values
Constant m=, dx=; //m: function number
Parameter p();//=[0.5,100];
Variable t, x[OutPut], y[OutPut], z[OutPut];
Minimum;
StartProgram [Pascal];
Procedure MainModel;
var i, j: integer;
y1_2, y2_2, y3_2: double;
y11, y21, y31: double;
rk1, rk2, rk3, rk4, y_y, yy, dydx: array[..m-] of double;
x, tx: double;
Procedure Drive(x_d: double; y_d: array of double; var dydx_d: array of double);
begin
dydx_d[] := p1*(-y_d[])+p2*(p3-y_d[]);
dydx_d[] := p1*(y_d[]-y_d[])+p2*(p3-y_d[]);
dydx_d[] := p1*(y_d[]-y_d[])+p2*(p3-y_d[]);
end;
begin
yy[]:=y1_1;
yy[]:=y2_1;
yy[]:=y3_1;
y_y[]:=y1_1;
y_y[]:=y2_1;
y_y[]:=y3_1;
for i := to DataLength- do begin
for j := to m- do y_y[j] := yy[j];
tx := dx*i;
x := tx;
Drive(x, yy, dydx);
for j := to m- do rk1[j] := dydx[j];
x := tx+0.5*dx;
for j := to m- do
yy[j] := y_y[j] + 0.5*rk1[j]*dx;
Drive(x, yy, dydx);
for j := to m- do
rk2[j] := dydx[j];
x := tx+0.5*dx;
for j := to m- do
yy[j] := y_y[j] + 0.5*rk2[j]*dx;
Drive(x, yy, dydx);
for j := to m- do
rk3[j] := dydx[j];
x := tx+dx;
for j := to m- do
yy[j] := y_y[j] + rk3[j]*dx;
Drive(x, yy, dydx);
for j := to m- do
rk4[j] := dydx[j];
for j := to m- do
yy[j] := y_y[j] + (/)*(rk1[j]+*rk2[j]+*rk3[j]+rk4[j])*dx;
x[i] := yy[];
y[i] := yy[];
z[i] := yy[];
end;
end;
EndProgram;
Data;
//0 10 15 20
11.80 15.37 20.77
14.04 17.04 22.23
15.44 16.85 21.16
17.80 18.07 22.49
18.60 19.43 24.46
19.22 20.12 23.58
21.77 21.83 24.51
22.17 22.11 25.38
23.41 23.37 25.53
23.17 24.92 26.69
24.56 25.55 29.21
25.85 26.06 28.00
24.64 28.94 30.32
25.15 28.94 30.41
26.92 30.06 31.87
26.37 29.29 31.87
26.71 31.48 34.60
26.61 31.86 33.57
27.13 33.84 36.75
29.32 32.95 36.36
28.24 33.03 36.23
28.42 32.50 36.01
28.11 33.12 39.19
28.98 35.32 37.29
30.23 35.56 37.79
30.21 34.86 42.05
29.11 35.40 42.67
30.42 36.20 40.74
28.84 35.91 40.53
29.44 36.50 43.33
5.0版本例子
使用1.5的小例子：
比如我知道一个二次函数的一系列离散点，相求这个二次函数的系数：
Title "Type your title here";
//Parameters ;
//Variable ;
//Function ;
Variable x,y;
Parameters a,b,c;
Function y=a*x*x+b*x+c;
Data;
9.8
23.3
42.8
68.3
99.8
137.3
180.8
230.3
285.8
347.3
414.8
488.3
567.8
653.3
744.8
842.3
945.8
1055.3
1170.8
1292.3
1419.8
1553.3
1692.8
1838.3
1989.8
2147.3
2310.8
2480.3
2655.8
2837.3
3024.8
是不是挺完美的，哈哈
实战：
Title "Type your title here";
//Parameters ;
//Variable ;
//Function ;
Parameters a,b,c;
Variable x,y;
Function y=-(a+b*c)*LN(c-x)-b*x+(a+b*c)*LN(c);
Data;
//x y
38.3612 38.09276591
71.3481 96.65369746
101.926 154.780949
133.801 156.4442223
164.935 334.699552
195.698 594.6346182
232.391 1529.185089
论文取点
蓝线是离散点，红线是计算点