欢迎访问 Forcal程序设计

FORCAL数值计算扩展动态库Fc IMSL V1.0 测试版

目 录

1 什么是FcIMSL
2 IMSL函数
错误输出 |
IMSL::ClearImslErr() | 清空ISML错误输出文件“ImslErr.txt”的内容
IMSL::ERSET | 重置错误处理设置
IMSL::IERCD() | 返回信息错误代码，无错时返回0
IMSL::N1RTY(i) | 返回信息错误类型
指数积分和相关函数 |
IMSL::EI(x) | x!=0，否则将返回运行错误
IMSL::SI(x) |
错误函数 |
IMSL::ERF(x) |
Bessel函数 |
IMSL::BSJ0(x) | 第一类0阶Bessel函数。运行错误：1：找不到相应的IMSL函数。
IMSL::BSJ1(x) | 第一类1阶Bessel函数。运行错误：1：找不到相应的IMSL函数。
IMSL::BSY0(x) | 第二类0阶Bessel函数。运行错误：1：找不到相应的IMSL函数。
IMSL::BSY1(x) | 第二类1阶Bessel函数。运行错误：1：找不到相应的IMSL函数。
IMSL::BSI0(x) | 修正的第一类0阶Bessel函数。运行错误：1：找不到相应的IMSL函数。
IMSL::BSI1(x) | 修正的第一类1阶Bessel函数。运行错误：1：找不到相应的IMSL函数。
IMSL::BSK0(x) | 修正的第三类0阶Bessel函数。运行错误：1：找不到相应的IMSL函数。
IMSL::BSK1(x) | 修正的第三类1阶Bessel函数。运行错误：1：找不到相应的IMSL函数。
IMSL::BSJNS(x,n,ns) | 计算一列第一类Bessel函数，阶数为整数，参数为实数。其中x的绝对值要小于100000；ns为长度不小于n的数组，用于返回计算值，第i个参数将存放第i阶Bessel函数值。 运行错误：1：找不到相应的IMSL函数；2：非法的数组长度；3：实数数组不存在；4：数组太小。
IMSL::CBJNS(x,y,n,ns) | 计算一列第一类Bessel函数，阶数为整数，参数为复数(x,y)。其中复数(x,y)的绝对值要小于10000；ns为长度不小于n×2的数组，用于返回计算值，第i×2和i×2+1个参数将存放第i阶Bessel函数值，因为是复数，占两个存储单元。 运行错误：1：找不到相应的IMSL函数；2：非法的数组长度；3：实数数组不存在；4：数组太小。
IMSL::BSINS(x,n,ns) | 计算一列修正的第一类Bessel函数，阶数为整数，参数为实数。其中ns为长度不小于n的数组，用于返回计算值，第i个参数将存在第i阶Bessel函数值。 运行错误：1：找不到相应的IMSL函数；2：非法的数组长度；3：实数数组不存在；4：数组太小。
IMSL::CBINS(x,y,n,ns) | 计算一列修正的第一类Bessel函数，阶数为整数，参数为复数(x,y)。其中复数(x,y)的绝对值要小于10000；ns为长度不小于n×2的数组，用于返回计算值，第i×2和i×2+1个参数将存放第i阶Bessel函数值，因为是复数，占两个存储单元。 运行错误：1：找不到相应的IMSL函数；2：非法的数组长度；3：实数数组不存在；4：数组太小。
IMSL::BSJS(xun,x,n,ns) | 计算阶和参数均为实数的Bessel函数数列。其中0<=xun<1，x为非负数，ns为长度不小于n的数组，用于返回计算值，第i个参数将存放第i+xun阶Bessel函数值。 运行错误：1：找不到相应的IMSL函数；2：非法的数组长度；3：实数数组不存在；4：数组太小。
IMSL::BSYS(xun,x,n,ns) | 计算阶和参数均为实数的Bessel函数数列。其中0<=xun<1，x为正数，ns为长度不小于n的数组，用于返回计算值，第i个参数将存放第i+xun阶Bessel函数值。 运行错误：1：找不到相应的IMSL函数；2：非法的数组长度；3：实数数组不存在；4：数组太小。
IMSL::BSIS(xun,x,n,ns) | 计算阶和参数均为实数的Bessel函数数列。其中0<=xun<1，x为实数，ns为长度不小于n的数组，用于返回计算值，第i个参数将存放第i+xun阶Bessel函数值。 运行错误：1：找不到相应的IMSL函数；2：非法的数组长度；3：实数数组不存在；4：数组太小。
IMSL::BSKS(xun,x,n,ns) | 计算阶和参数均为实数的Bessel函数数列。其中0<xun<1，x为实数，ns为长度不小于n的数组，用于返回计算值，第i个参数将存放第i+xun阶Bessel函数值。 运行错误：1：找不到相应的IMSL函数；2：非法的数组长度；3：实数数组不存在；4：数组太小。
IMSL::CBJS(xun,x,y,n,ns) | 计算阶为实数，参数为复数的Bessel函数数列。其中-0.5<xun，x为复数(x,y)，ns为长度不小于n×2的数组，用于返回计算值，第i个参数将存放第i+xun阶Bessel函数值。 运行错误：1：找不到相应的IMSL函数；2：非法的数组长度；3：实数数组不存在；4：数组太小。
IMSL::CBYS(xun,x,y,n,ns) | 计算阶为实数，参数为复数的Bessel函数数列。其中-0.5<xun，x为复数(x,y)，ns为长度不小于n×2的数组，用于返回计算值，第i个参数将存放第i+xun阶Bessel函数值。 运行错误：1：找不到相应的IMSL函数；2：非法的数组长度；3：实数数组不存在；4：数组太小。
IMSL::CBIS(xun,x,y,n,ns) | 计算阶为实数，参数为复数的Bessel函数数列。其中-0.5<xun，x为复数(x,y)，ns为长度不小于n×2的数组，用于返回计算值，第i个参数将存放第i+xun阶Bessel函数值。 运行错误：1：找不到相应的IMSL函数；2：非法的数组长度；3：实数数组不存在；4：数组太小。
IMSL::CBKS(xun,x,y,n,ns) | 计算阶为实数，参数为复数的Bessel函数数列。其中-0.5<xun，x为复数(x,y)，ns为长度不小于n×2的数组，用于返回计算值，第i个参数将存放第i+xun阶Bessel函数值。 运行错误：1：找不到相应的IMSL函数；2：非法的数组长度；3：实数数组不存在；4：数组太小。
插值 |
IMSL::CSIEZ | cubic spline 三次样条插值
积分 |
IMSL::QDAGS | 计算一个函数的积分（可能在端点存在奇异）
IMSL::QDAG | 基于Gauss-Kronrod法则使用自适应算法计算一个函数的积分
IMSL::QDAGP | 计算一个给定奇异点的函数的积分
IMSL::QDAGI | 计算一个函数在无穷区间或半无穷区间的积分
IMSL::TWODQ | 计算二重积分
IMSL::QAND | 计算函数在超矩形下的积分
微分 |
IMSL::DERIV | 计算一元函数的一阶、二阶或三阶导数
常微分方程 |
IMSL::IVPRK | 求解一次常微分方程[y'=f(t,y) ,y(t0)=y0]，Runge-Kutta方法
IMSL::IVMRK | 求解一次常微分方程[y'=f(t,y) ,y(t0)=y0]，Runge-Kutta方法
IMSL::IVPAG | 求解一次常微分方程[Ay'=f(t,y) ,y(t0)=y0]或[y'=f(t,y) ,y(t0)=y0]，Adams 或者 Gear 方法
偏微分方程 |
IMSL::FPS2H | 在一个二维的矩形区域求解Poisson或Helmholtz方程
IMSL::FPS3H | 在一个三维的矩形区域求解Poisson或Helmholtz方程
IMSL::MOLCH | 用线性方法求解形如 u t= f(x, t, u, u x, uxx)的偏微分方程，计算的结果是一个三次的Hermite多项式
IMSL::UMAG_DE | 获取或设置偏微分方程求解时的选项
公用函数 |
IMSL::MACH(i) | 返回双精度实型常数
IMSL::IFNAN(x) | 判断x是否是NAN，NAN表示不是一个数
3 如何封装IMSL库函数
3.1 用IMSL32.dll输出IMSL库函数 |
3.2 FcIMSL32W.dll中封装IMSL库函数的例子 |

1 什么是FcSyIMSL [返回页首]

FcIMSL32W.dll是一个Forcal数值计算扩展动态库，该库对IMSL库中的双精度函数进行了封装。

在FcIMSL中的函数是通过二级函数命名空间“IMSL”输出的，所有函数均具有类似“IMSL::FPS3H(...)”的格式 ，都是实数函数。使用!using("IMSL");可简化FcIMSL中的函数访问。

FcIMSL32W.dll需要FcData32W.dll和IMSL32.dll（或者是其他包含IMSL库函数的支持库）的支持。FcData32W.dll要先于FcIMSL32W.dll加载；IMSL32.dll要在搜索目录中，或者需要指定IMSL32.dll的位置。

能够被FcIMSL32W.dll使用的IMSL32.dll（或者是其他包含IMSL库函数的支持库）应满足以下要求：用Fortran创建该动态库，输出函数名以“F_”开头，且是双精度函数。例如FcIMSL支持函数FPS3H，在IMSL32.dll中输出的函数名应是F_DFPS3H，且使用双精度函数DFPS3H进行计算。

显式加载IMSL32.dll的方法：在加载Forcal扩展动态库列表中使用如下方法，注意冒号后即IMSL32.dll的路径及文件名，也可以使用其他包含IMSL库函数的支持库。

"dll\FcIMSL32W.dll:dll\IMSL32.dll"

一旦FcIMSL32W.dll加载成功，将在当前目录创建或打开文件“ImslErr.txt”，调用ISML函数的所有错误都自动输出到该文件。FcIMSL32W.dll在工作时，不能从外部清空或删除文件“ImslErr.txt”。

如果FcIMSL32W.dll在IMSL32.dll中找不到要调用的函数，所有函数都将返回运行错误代码1。

注意1：IMSL库函数要求传递Fortran格式的数组。可以使用函数 sys::rearray 在Forcal数组（C数组格式）和Fortran数组之间进行转换。

注意2：有些IMSL库函数不允许嵌套调用，在需要嵌套调用这些库函数时，请同时加载 FcIMSL32W1.dll，该库的输出函数与 FcIMSL32W.dll 完全相同，但使用命名空间“IMSL1”：

"dll\FcIMSL32W1.dll:dll\IMSL321.dll"

警告：IMSL32.dll是由Fortran生成的，而Fortran中没有异常捕获机制，无法回卷堆栈，一旦发生了计算错误，都是直接终止程序。

FcIMSL32W.dll封装的IMSL库函数仅作演示学习之用 ，参考IMSL库的说明文档能更好地理解这些函数的用法。IMSL库是Visual Numerics公司的产品，在使用前须获得该公司的使用许可。

2 IMSL函数 [返回页首]

[返回页首][实数函数] IMSL::ERSET(IERSVR, IPACT, ISACT)：重置错误处理设置

IERSVR：IMSL错误级别为1～5，数字越大，错误越严重。若IERSVR=0，包括了所有错误；否则仅包含自IERSVR以上更严重的错误。
IPACT：IPACT=-1：不改变当前设置；IPACT=0：不输出错误；IPACT=1：输出错误；IPACT=2：重置缺省设置。
IPACT：IPACT=-1：不改变当前设置；IPACT=0：不停止程序；IPACT=1：忽略该参数；IPACT=2：重置缺省设置。

[返回页首][实数函数] IMSL::CSIEZ(NDATA, XDATA, FDATA, XVEC1,XVEC2,...,XVECm, &VALUE1,&VALUE2,...,&VALUEm)：cubic spline 三次样条插值

NDATA：输入的数据点数目。
XDATA： 长度不小于 NDATA 的数组。存放数据点的X轴值。
FDATA：长度不小于 NDATA 的数组。存放数据点的Y轴值。
XVEC1,XVEC2,...,XVECm：所要插值的X轴位置，共m个位置。
&VALUE1,&VALUE2,...,&VALUEm：返回插值得到的m个结果，须使用取地址运算符“&”。

运行错误：1：找不到相应的IMSL函数；2：数据点太少或参数不匹配；3：数组不存在；4：数组太小。

[例子]

> main(:x,y,v1,v2)=
> {
> x=new[rtoi(real_s),rtoi(5),rtoi(EndType),1.615, 1.634, 1.702, 1.828, 1.921],
> y=new[rtoi(real_s),rtoi(5),rtoi(EndType),2.4145, 2.46459, 2.65271, 3.03035, 3.34066],
> IMSL::CSIEZ(5,x,y,1.682,1.813,&v1,&v2),
> printff{"函数f(1.682)={1,r}\r\n",v1},
> printff{"函数f(1.813)={1,r}\r\n",v2},
> delete[x],delete[y]
> };

[返回页首][实数函数] IMSL::QDAGS(F,A,B,ERRABS,ERRREL,&ERREST)：计算一个函数的积分（可能在端点存在奇异）

F：Forcal一元函数句柄。该函数由用户定义。
A,B：积分下限及积分上限。
ERRABS： 期望的绝对精度。ERRABS>=0。
ERRREL：期望的相对精度。ERRREL>=0。
ERREST： 返回计算的绝对误差。
返回值：积分值。

运行错误：1：找不到相应的IMSL函数；2：指定的表达式不存在；3：不是一元函数；4：参数不符合要求 ；5：内存错误；6：不允许嵌套调用。

[例子] 求函数f(x)=ln(x)/sqrt(x)在0～1上的 积分值

> f(x)=ln(x)/sqrt(x);
> IMSL::QDAGS[HFor("f"),0,1,0,0.001,0]; //得到积分值
> (:err)=IMSL::QDAGS[HFor("f"),0,1,0,0.001,&err],err; //得到积分值的绝对误差

[返回页首][实数函数] IMSL::QDAG(F,A,B,ERRABS,ERRREL,IRULE,&ERREST)：基于Gauss-Kronrod法则使用自适应算法计算一个函数的积分

F：Forcal一元函数句柄。该函数由用户定义。
A,B：积分下限及积分上限。
ERRABS：期望的绝对精度。ERRABS>=0。
ERRREL：期望的相对精度。ERRREL>=0。
IRULE：积分法则。Gauss-Kronrod法则用于下面的点数（IRULE值：点数）：（1：7～15）、（2：11～21）、（3：15～31）、（4：20～41）、（5：25～51）、（6：30～61）。大多数函数推荐取IRULE=2；如果函数有一个奇异的峰值，取IRULE=1；如果函数是摆动的，取IRULE=6。
ERREST：返回计算的绝对误差。
返回值：积分值。

运行错误：1：找不到相应的IMSL函数；2：指定的表达式不存在；3：不是一元函数；4：参数不符合要求 ；5：内存错误；6：不允许嵌套调用。

[返回页首][实数函数] IMSL::QDAGP(F,A,B,ERRABS,ERRREL,&ERREST,p1,p2,...)：计算一个给定奇异点的函数的积分

F：Forcal一元函数句柄。该函数由用户定义。
A,B：积分下限及积分上限。
ERRABS：期望的绝对精度。ERRABS>=0。
ERRREL：期望的相对精度。ERRREL>=0。
ERREST：返回计算的绝对误差。
p1,p2,...：积分区间的间断点，至少要有一个间断点。通常这些间断点是被积函数出现奇异值的地方。
返回值：积分值。

运行错误：1：找不到相应的IMSL函数；2：参数太少（至少7个参数）；3：指定的表达式不存在；4：不是一元函数；5：参数不符合要求 ；6：内存错误；7：不允许嵌套调用。

[返回页首][实数函数] IMSL::QDAGI(F,BOUND,INTERV,ERRABS,ERRREL,&ERREST)：计算一个函数在无穷区间或半无穷区间的积分

F：Forcal一元函数句柄。该函数由用户定义。
BOUND： 积分部分有限区间的边界，如果 INTERV=2则被忽略。
INTERV：指示积分区间的标志。如果 INTERV=-1：（-∞，BOUND）；如果 INTERV=1：（BOUND，+∞）；如果 INTERV=2：（-∞，+∞）。
ERRABS：期望的绝对精度。ERRABS>=0。
ERRREL：期望的相对精度。ERRREL>=0。
ERREST：返回计算的绝对误差。
返回值：积分值。

运行错误：1：找不到相应的IMSL函数；2：指定的表达式不存在；3：不是一元函数；4：参数不符合要求 ；5：内存错误；6：不允许嵌套调用。

[返回页首][实数函数] IMSL::TWODQ(F, A, B, G, H, ERRABS, ERRREL, IRULE, ERREST)：计算二重积分

F：Forcal二元 被积函数句柄。该函数由用户定义。
A,B：积分下限与积分上限。
G,H：Forcal一元函数句柄，由用户定义 ，用来计算内层积分的下限与上限。
ERRABS：期望的绝对精度。
ERRREL：期望的相对精度。
IRULE：选择积分的法则。Gauss-Kronrod法则用于下面的点数（IRULE值：点数）：（1：7～15）、（2：11～21）、（3：15～31）、（4：20～41）、（5：25～51）、（6：30～61）。大多数函数推荐取IRULE=2；如果函数有一个奇异的峰值，取IRULE=1；如果函数是摆动的，取IRULE=6。
ERREST：返回计算的绝对误差。
返回值：积分值。

运行错误：1：找不到相应的IMSL函数；2：指定的表达式不存在；3： 自定义函数参数不符合要求；4：内存错误；5：不允许嵌套调用。

[例子] 计算函数f(x,y)=y*cos(x+y*y)的近似积分值。外层区间取[0,1]，内层区间取[-2*x,5*x]。

> f(x,y)=y*cos(x+y*y);
> g(x)=-2*x;
> h(x)=5*x;
> IMSL::TWODQ[HFor("f"),0,1,HFor("g"),HFor("h"),0.001,0,6,0];

[返回页首][实数函数] IMSL::QAND(F, A1, A2, ..., An, B1, B2, ..., Bn, ERRABS, ERRREL, MAXFCN, ERREST)：计算函数在超矩形下的积分

F：Forcal多元函数句柄 ，设有n个自变量。该函数由用户定义。
A1, A2, ..., An：积分下限，有n个下限值。
B1, B2, ..., Bn：积分上限，有n个上限值。
ERRABS：期望的绝对精度。
ERRREL：期望的相对精度。
MAXFCN： 可允许的函数估计的最大数目。必须 MAXFCN<=256^n。
ERREST：返回计算的绝对误差。
返回值：积分值。

运行错误：1：找不到相应的IMSL函数；2：指定的表达式不存在；3：不是 多元函数；4：参数不匹配；5：内存错误；6：不允许嵌套调用。

[例子] 计算函数exp[-(x1^2+x2^2+x3^2)]在扩展的立方体（整个三维空间）上的积分的近似值。积分区间为[-2,-2,-2～2,2,2]。

> f(x1,x2,x3)=exp[-(x1^2+x2^2+x3^2)];
> IMSL::QAND[HFor("f"),-2,-2,-2,2,2,2,0.0001,0.001,100000,0]; //得到积分值
> (:err)=IMSL::QAND[HFor("f"),-2,-2,-2,2,2,2,0.0001,0.001,100000,&err],err; //得到积分值的绝对误差

[返回页首][实数函数] IMSL::DERIV(FCN,KORDER,X,BGSTEP,TOL)：计算一元函数的一阶、二阶或三阶导数

FCN：Forcal一元函数句柄，要计算该函数在X的导数。该函数由用户定义。
KORDER：导数的阶（1、2或3）。
X：被求导的点。
BGSTEP：用来计算步长的起始值，这个步长被用于计算导数。BGSTEP>0。
TOL：允许的相对误差（0～1）。
返回值：函数的一阶、二阶或三阶导数。

运行错误：1：找不到相应的IMSL函数；2：指定的表达式不存在；3：不是一元函数；4： 参数不符合要求；5：内存错误；6：不允许嵌套调用。

[例子] 求函数f(x)=2*x^4+3*x在0.75处的三阶导数

> f(x)=2*x^4+3*x;
> IMSL::DERIV[HFor("f"),3,0.75,0.1,0.01];

[返回页首][实数函数] IMSL::IVPRK(&IDO, FCN, &T, TEND, TOL, PARAM, &Y1,&Y2,...,&Yn)：求解一次常微分方程[y'=f(t,y) ,y(t0)=y0]，Runge-Kutta方法

IDO：指示计算的状态。IDO=1：初次进入；IDO=2：再次正常进入；IDO=3：最后调用，释放工作空间。一般情况下，最初调用时将IDO置为1，然后函数自动将它置为2，这个值将一直使用，直到最后一次调用置为3。为了使函数自动将IDO置为2，需要使用取地址运算符“&”。
FCN：用户定义的函数，必须由二级函数HFor获得该函数的句柄，用于计算微分方程组中各方程右端函数值，由用户自编。该表达式有2*n+1个参数，第一个参数为自变量，随后n个参数为函数值，最后n个参数为右端函数值（即微分方程的值）。n为微分方程组中方程的个数，也是未知函数的个数。形式如下：

> f(t,y1,y2,y3,d1,d2,d3)=
> { d1=y2,
> d2=-y1,
> d3=-y3
> };

T：独立变量t。输入初始时间t0，返回积分被更新的时间，一般，这个新值是TEND。为了使函数自动更新T，需要使用取地址运算符“&”。
TEND：时间变量终值，用户希望得到这一点的解。
TOL：允许误差。
PARAM： 长度不小于50的数组。用来控制函数 IVPRK的行为，数组中每个元素值都有它代表的设置值，当元素值为0时，会使用默认值。该参数可以为0，则全部使用默认值。
&Y1,&Y2,...,&Yn：存放n个未知函数在起始点处的函数值Y(n)=Yn(T)，返回时存放微分方程的解Y(n)=Yn(TEND)。为了使函数 能返回值，需要使用取地址运算符“&”。

运行错误：1：找不到相应的IMSL函数；2：指定的表达式不存在；3： 方程参数不符合要求；4：数组不存在；5：不是一维数组或数组长度太小；6：参数不符合要求；7：不能嵌套调用该函数。

[例子] 设一阶微分方程组及初值为：
r'=2r-2rf, r(0)=1
f'=-f+rf, f(0)=3
计算t=1,2,...,10时的r、f的值。
程序如下：

> f(t,r,f,dr,df)={dr=2*r-2*r*f, df=-f+r*f}; //函数定义
> main(:IDO,T,TEND,r,f)=
> {
> IMSL::ClearImslErr(),
> printff["\r\n时间 r f\r\n"],
> IDO=1, T=0, TEND=1, r=1, f=3,
> while{TEND<=10,
> IMSL::IVPRK[&IDO,HFor("f"),&T,TEND,1e-5,0,&r,&f],
> printff["{1,r,-14.6}{2,r,-14.6}{3,r,-14.6}\r\n",T,r,f],
> TEND++
> },
> IMSL::IVPRK[3,HFor("f"),&T,TEND,1e-5,0,&r,&f] //释放空间
> };

[返回页首][实数函数] IMSL::IVMRK(&IDO, FCN, &T, TEND, &Y1,&Y2,...,&Yn, &YY1,&YY2,...,&YYn)：求解一次常微分方程[y'=f(t,y) ,y(t0)=y0]，Runge-Kutta方法

IDO：指示计算的状态。IDO=1：初次进入；IDO=2：再次正常进入；IDO=3：最后调用，释放工作空间。一般情况下，最初调用时将IDO置为1，然后函数自动将它置为2，这个值将一直使用，直到最后一次调用置为3。为了使函数自动将IDO置为2，需要使用取地址运算符“&”。
FCN：用户定义的函数，必须由二级函数HFor获得该函数的句柄，用于计算微分方程组中各方程右端函数值，由用户自编。该表达式有2*n+1个参数，第一个参数为自变量，随后n个参数为函数值，最后n个参数为右端函数值（即微分方程的值）。n为微分方程组中方程的个数，也是未知函数的个数。形式如下：

> f(t,y1,y2,y3,d1,d2,d3)=
> { d1=y2,
> d2=-y1,
> d3=-y3
> };

T：独立变量t。输入初始时间t0，返回积分被更新的时间，一般，这个新值是TEND。为了使函数自动更新T，需要使用取地址运算符“&”。
TEND：时间变量终值，用户希望得到这一点的解。
&Y1,&Y2,...,&Yn：存放n个未知函数在起始点处的函数值Y(n)=Yn(T)，返回时存放微分方程的解Y(n)=Yn(TEND)。为了使函数 能返回值，需要使用取地址运算符“&”。
&YY1,&YY2,...,&YYn：返回时存放微分方程的解Yn'(TEND)。为了使函数 能返回值，需要使用取地址运算符“&”。

运行错误：1：找不到相应的IMSL函数；2：指定的表达式不存在；3： 方程参数不符合要求；4：参数不符合要求；5：不能嵌套调用该函数。

[例子] 设一阶微分方程组及初值为：
r'=2r-2rf, r(0)=1
f'=-f+rf, f(0)=3
计算t=1,2,...,10时的r、f的值。
程序如下：

> f(t,r,f,dr,df)={dr=2*r-2*r*f, df=-f+r*f}; //函数定义
> main(:IDO,T,TEND,r,f,rr,ff)=
> {
> IMSL::ClearImslErr(),
> printff["\r\n时间 r f r' f'\r\n"],
> IDO=1, T=0, TEND=1, r=1, f=3,
> while{TEND<=10,
> IMSL::IVMRK[&IDO,HFor("f"),&T,TEND,&r,&f,&rr,&ff],
> printff["{1,r,-14.6}{2,r,-14.6}{3,r,-14.6}{4,r,-14.6}{5,r,-14.6}\r\n",T,r,f,rr,ff],
> TEND++
> },
> IMSL::IVMRK[3,HFor("f"),&T,TEND,&r,&f,&rr,&ff] //释放空间
> };

[返回页首][实数函数] IMSL::IVPAG(&IDO, FCN, FCNJ, A, &T, TEND, TOL, PARAM, &Y1,&Y2,...,&Yn)：求解一次常微分方程[Ay'=f(t,y) ,y(t0)=y0]或[y'=f(t,y) ,y(t0)=y0]，Adams 或者 Gear 方法

IDO：指示计算的状态。IDO=1：初次进入；IDO=2：再次正常进入；IDO=3：最后调用，释放工作空间。一般情况下，最初调用时将IDO置为1，然后函数自动将它置为2，这个值将一直使用，直到最后一次调用置为3。为了使函数自动将IDO置为2，需要使用取地址运算符“&”。
FCN：用户定义的函数，必须由二级函数HFor获得该函数的句柄，用于计算微分方程组中各方程右端函数值，由用户自编。该表达式有2*n+1个参数，第一个参数为自变量，随后n个参数为函数值，最后n个参数为右端函数值（即微分方程的值）。n为微分方程组中方程的个数，也是未知函数的个数。形式如下：

> f(t,y1,y2,y3,d1,d2,d3)=
> { d1=y2,
> d2=-y1,
> d3=-y3
> };

FCNJ：用户定义的函数，必须由二级函数HFor获得该函数的句柄，用于计算 Jacobian 矩阵，由用户自编 。该函数仅在数组PARAM的第13个参数（即PARAM[12]）为1时提供。该函数有1+n+m个参数，第一个参数为自变量，随后n个参数为函数值，最后m个参数返回 Jacobian 矩阵的值。n为微分方程组中方程的个数 。m的大小与数组PARAM的第14个参数（即PARAM[13]，记为MTYPE）有关，另记NLC为数组PARAM的第15个参数（即PARAM[14]），NUC为数组PARAM的第16个参数（即PARAM[15]）；若MTYPE=0或MTYPE=2，则m=n*n；若MTYPE=1，则m=(2 * NLC + NUC + 1)*(2 * NLC + NUC + 1)；若MTYPE=3，则m=(NUC + 1)*(NUC + 1)；注意：Jacobian 矩阵的m个元素在 J1,J2,... ...,Jm 中是按列存放的。FCNJ 形式如下：

> J(t,y1,y2,y3,J1,J2,... ...,Jm)=
> { J1=y2,
> J2=-y1,
> ... ...
> Jm=-y3
> };

A：矩阵k*n。该矩阵仅在数组PARAM的第19个参数（即PARAM[18]）非0时提供。n为微分方程组中方程的个数。若数组PARAM的第20个参数（即PARAM[19]）非0，则k=PARAM[19]，否则k的值与MTYPE有关：若MTYPE=0或MTYPE=2，则k=n；若MTYPE=1，则k=(NLC + NUC + 1)；若MTYPE=3，则k=(NUC + 1)。
T：独立变量t。输入初始时间t0，返回积分被更新的时间，一般，这个新值是TEND。为了使函数自动更新T，需要使用取地址运算符“&”。
TEND：时间变量终值，用户希望得到这一点的解。
TOL：允许误差。
PARAM： 长度不小于50的数组。用来控制函数 IVPAG的行为，数组中每个元素值都有它代表的设置值，当元素值为0时，会使用默认值。该参数可以为0，则全部使用默认值。
&Y1,&Y2,...,&Yn：存放n个未知函数在起始点处的函数值Y(n)=Yn(T)，返回时存放微分方程的解Y(n)=Yn(TEND)。为了使函数能返回值，需要使用取地址运算符“&”。

运行错误：1：找不到相应的IMSL函数；2：指定的表达式不存在；3：方程参数不符合要求；4：数组不存在；5：不是一维数组或数组长度太小；6：参数不符合要求 ；7：不能嵌套调用该函数。

[例子1] 设一阶微分方程组及初值为：
y1'=y2*y3, y1(0)=0
y2'=-y1*y3, y2(0)=1
y3'=-0.51*y1*y2, y3(0)=1
计算t=1,2,...,10时的y1、y2、y3的值。
程序如下：

> f(t,y1,y2,y3,dy1,dy2,dy3)= //函数定义
> {
> dy1=y2*y3,
> dy2=-y1*y3,
> dy3=-0.51*y1*y2
> };
> main(:IDO,T,TEND,y1,y2,y3)=
> {
> IMSL::ClearImslErr(),
> printff["\r\n时间 y1 y2 y3\r\n"],
> IDO=1, T=0, TEND=1, y1=0, y2=1, y3=1,
> while{TEND<=10,
> IMSL::IVPAG[&IDO,HFor("f"),0,0,&T,TEND,1e-5,0,&y1,&y2,&y3],
> printff["{1,r,-14.6}{2,r,-14.6}{3,r,-14.6}{4,r,-14.6}\r\n",T,y1,y2,y3],
> TEND++
> },
> IMSL::IVPAG[3,HFor("f"),0,0,&T,TEND,1e-5,0,&y1,&y2,&y3] //释放空间
> };

[例子2] 设一阶刚性微分方程组及初值为：
y1'=-c1*y1+c2*y2*y3, y1(0)=1
y2'=-dy1-dy3, y2(0)=0
y3'=c3*y2*y2, y3(0)=0
其中c1=0.04,c2=10000,c3=3e7。计算t=1,2,...,10时的y1、y2、y3的值。
程序如下：

> f(t,y1,y2,y3,dy1,dy2,dy3:c1,c2,c3)= //函数定义
> {
> c1=0.04, c2=10000, c3=3e7,
> dy1=-c1*y1+c2*y2*y3,
> dy3=c3*y2*y2,
> dy2=-dy1-dy3
> };
> j(t,y1,y2,y3,J11,J21,J31,J12,J22,J32,J13,J23,J33:c1,c2,c3)=
> {
> c1=0.04, c2=10000, c3=3e7,
> J11=-c1,
> J21=-J11,
> J31=0,
> J12=c2*y3,
> J22=-J12-J32,
> J32=2*c3*y2,
> J13=c2*y2,
> J23=-J13,
> J33=0
> };
> main(:IDO,T,TEND,PARAM,i,y1,y2,y3)=
> {
> IMSL::ClearImslErr(),
> PARAM=new[rtoi(real_s),rtoi(50)],
> i=0, while{i<50, PARAM.XSLSF::setrai[i,0] ,i++}, //参数初始化
> PARAM.XSLSF::setrai[9,1], //选择绝对错误控制
> PARAM.XSLSF::setrai[11,2], //选择BDF方法
> PARAM.XSLSF::setrai[12,1], //选择chord方法，用函数j提供Jacobian矩阵
> printff["\r\n时间 y1 y2 y3\r\n"],
> IDO=1, T=0, TEND=1, y1=1, y2=0, y3=0,
> while{TEND<=10,
> IMSL::IVPAG[&IDO,HFor("f"),HFor("j"),0,&T,TEND,1e-5,PARAM,&y1,&y2,&y3],
> printff["{1,r,-14.6}{2,r,-14.6}{3,r,-14.6}{4,r,-14.6}\r\n",T,y1,y2,y3],
> TEND++
> },
> IMSL::IVPAG[3,HFor("f"),HFor("j"),0,&T,TEND,1e-5,PARAM,&y1,&y2,&y3], //释放空间
> delete[PARAM]
> };

[返回页首][实数函数] IMSL::FPS2H(PRHS,BRHS,COEFU,AX,BX,AY,BY,IBCTY1,IBCTY2,IBCTY3,IBCTY4,IORDER,U)：在一个二维的矩形区域求解Poisson或Helmholtz方程

求解如下形式的偏微分方程：

其中c为常数，p为x和y的函数。

PRHS：用户定义的函数，计算偏微分方程的右端项。函数形式为：PRHS(x,y)
BRHS：用户定义的函数，计算边界条件的右端项。函数形式为：BRHS(ISIDE,x,y)，其中ISIDE为边界数目。
COEFU：微分方程中u的系数的值。
AX,BX,AY,BY：区域X的左边、右边的值；区域Y的底边、顶边的值。
IBCTY1,IBCTY2,IBCTY3,IBCTY4：区域右边、底边、左边、顶边的边界条件类型，或者指示解是周期性的。有三类边界条件：
1：给定u的值（第一类边界条件或狄利克雷边界条件）。
2：给定du/dx（边1和/或3）、du/dy（边2和/或4）的值（诺依曼边界条件） 。
3：解是周期性的。
IORDER：精确的有限差分近似的阶，取2或4，一般取4。
U：大小为NX*NY的数组，返回网格点的解。必须NX>=4，NY>=4。 要求(BX - AX)/(NX - 1) = (BY -AY)/(NY - 1)。

运行错误：1：找不到相应的IMSL函数；2：指定的表达式不存在；3：不是二元或三元函数；4： 数组不存在；5：不是二维数组；6：参数不符合要求；7：内存错误；8：不允许嵌套调用。

[例子]

Forcal代码：

> i: OutMatrix(p:n1,n2,k,i,j)= //输出二维数组（矩阵）
> {
> FCDLen(p,&k,&n1,&n2),
> if[k!=2,printff{"\r\n*** 地址为{1,i}的对象不是一个矩阵！ ***\r\n",p.ptoi()},return(0)],
> i=0,(i<n1).while{
> printff{"\r\n"},
> j=0,(j<n2).while{printff{"{1,r,14.6}",get[p,i,j]},j++}, //可修改数据输出宽度和有效数字个数以满足实际要求
> i++
> },
> printff{"\r\n"}
> };
> PRHS(x,y)=-2*sin[x+2*y]+16*exp[2*x+3*y];
> BRHS(ISIDE,x,y)=
> {
> which{
> ISIDE==2 : 2*cos[x+2*y]+3*exp[2*x+3*y],
> sin[x+2*y]+exp[2*x+3*y]
> }
> };
> main(:U)=
> {
> U=new[rtoi(real_s),rtoi(33),rtoi(17)].sys::rearray(),
> IMSL::FPS2H(HFor("PRHS"),HFor("BRHS"),3,0,0.25,0,0.5,1,2,1,1,4,U),
> U.sys::rearray().OutMatrix[],
> delete[U]
> };

[返回页首][实数函数] IMSL::FPS3H(PRHS,BRHS,COEFU,AX,BX,AY,BY,AZ,BZ,IBCTY1,IBCTY2,IBCTY3,IBCTY4,IBCTY5,IBCTY6,IORDER,U)：在一个三维的矩形区域求解Poisson或Helmholtz方程

求解如下形式的偏微分方程：

其中c为常数，p为x、y和z的函数。

PRHS：用户定义的函数，计算偏微分方程的右端项。函数形式为：PRHS(x,y,z)
BRHS：用户定义的函数，计算边界条件的右端项。函数形式为：BRHS(ISIDE,x,y,z)，其中ISIDE为边界数目。
COEFU：微分方程中u的系数的值。
AX,BX,AY,BY,AZ,BZ：区域X的左边、右边的值 ；区域Y的底边、顶边的值；区域Z的前边、后边的值。
IBCTY1,IBCTY2,IBCTY3,IBCTY4,IBCTY5,IBCTY6：区域右边、底边、左边、顶边 、前边、后边的边界条件类型，或者指示解是周期性的。有三类边界条件：
1：给定u的值（第一类边界条件或狄利克雷边界条件）。
2：给定du/dx（边1和/或3）、du/dy（边2和/或4）的值（诺依曼边界条件） 、du/dz（边5和/或6）。
3：解是周期性的。
IORDER：精确的有限差分近似的阶，取2或4，一般取4。
U：大小为NX*NY*NZ的数组，返回网格点的解。必须NX>=4，NY>=4，NZ>=4。要求(BX - AX)/(NX - 1) = (BY -AY)/(NY - 1) = (BZ -AZ)/(NZ - 1)。

运行错误：1：找不到相应的IMSL函数；2：指定的表达式不存在；3：不是二元或三元函数；4： 数组不存在；5：不是三维数组；6：参数不符合要求；7：内存错误；8：不允许嵌套调用。

[例子]

Forcal代码：

> PRHS(x,y,z)=-4*cos[3*x+y-2*z]+12*exp[x-z]+10;
> BRHS(ISIDE,x,y,z)=
> {
> which{
> ISIDE==5 : -2*sin[3*x+y-2*z]-exp[x-z],
> cos[3*x+y-2*z]+exp[x-z]+1
> }
> };
> main(:U,t)=
> {
> U=new[rtoi(real_s),rtoi(17),rtoi(17),rtoi(9)].sys::rearray(),
> IMSL::FPS3H(HFor("PRHS"),HFor("BRHS"),10,0,0.125,0,0.25,0,0.25,1,1,1,1,2,1,4,U),
> U.sys::rearray(), t=U.XSLSF::getrai[16,16,8],
> delete[U], t
> };

[返回页首][实数函数] IMSL::MOLCH(&IDO,FCNUT,FCNBC,NPDES,&T,TEND,XBREAK,TOL,HINIT,Y)：用线性方法求解形如 u t= f(x, t, u, u x, uxx)的偏微分方程，计算的结果是一个三次的Hermite多项式

求解如下形式的偏微分方程：

IDO：指示计算的状态。IDO=1：初次进入；IDO=2：再次正常进入；IDO=3：最后调用，释放工作空间。一般情况下，最初调用时将IDO置为1，然后函数自动将它置为2，这个值将一直使用，直到最后一次调用置为3。为了使函数自动将IDO置为2，需要使用取地址运算符“&”。

FCNUT：用户定义的函数，计算ut。函数形式为（M=NPDES）：

FCNUT(x,t,u1,u2,...,uM,ux1,ux2,...,uxM,uxx1,uxx2,...,uxxM,ut1,ut2,...,utM)=
{
ut1=f1(x,t,u1,u2,...,uM,ux1,ux2,...,uxM,uxx1,uxx2,...,uxxM),
ut2=f2(x,t,u1,u2,...,uM,ux1,ux2,...,uxM,uxx1,uxx2,...,uxxM),
... ...
utM=fM(x,t,u1,u2,...,uM,ux1,ux2,...,uxM,uxx1,uxx2,...,uxxM)
};

FCNBC：用户定义的函数，计算边界条件。函数形式为（M=NPDES）：

FCNBC(x,t,α1,α2,...,αM,β1,β2,...,βM,γ1,γ2,...,γM)=
{
α1=...,
α2=...,
... ...
αM=...,
β1=...,
β2=...,
... ...
βM=...,
γ1=...,
γ2=...,
... ...
γM=...
};

NPDES：微分方程的个数。
T：独立变量t。输入初始时间t0，返回积分被更新的时间，一般，这个新值是TEND。 为了使函数自动更新T，需要使用取地址运算符“&”。
TEND：时间变量终值，用户希望得到这一点的解。
XBREAK：长为NX的数组，包含的是三维Hermite样条的断点，这些断点用来对x离散化。 这些点应是递增的。
TOL：允许误差。
HINIT：在t积分时的初始步长。
Y：大小为NPDES*NX（或NPDES*2*NX）的数组。在输入时，Y包含的是初始值，必须满足边界条件；输出时，返回计算结果。 如果Y的大小为NPDES*2*NX，须在 Y(k,i+NX)提供一阶导数初值 ux（x, t 0），返回时获得 ux（x, t end）。

运行错误：1：找不到相应的IMSL函数；2：指定的表达式不存在；3： 指定的函数参数不匹配；4：数组不存在；5：不是一维或二维数组；6：参数不符合要求；7：内存错误；8：不能嵌套调用该函数。

[例子1] 求解以下偏微分方程：

Forcal代码：

> FCNUT(x,t,u,ux,uxx,ut:D,V)=
> {
> which{
> x<=0.5 : [ D=5, V=1000 ],
> [ D=1, V=1 ]
> },
> ut=D*uxx-V*ux
> };
> FCNBC(x,t,α,β,γ)=
> {
> α=1, β=0, γ=0
> };
> main(:XBREAK,Y,NPDES,IDO,NX,TOL,HINIT,NSTEP,T,TEND,U0,i)=
> {
> IMSL::ClearImslErr(),
> NPDES=1, NX=100,
> XBREAK=new[rtoi(real_s),rtoi(NX)],
> Y=new[rtoi(real_s),rtoi(NX),rtoi(NPDES)],
> U0=1,
> i=0,(i<NX).while{
> XBREAK.XSLSF::setrai[i,i/(NX-1)],
> Y.XSLSF::setrai[i,0,0], //设置函数值
> i++
> },
> Y.XSLSF::setrai[0,0,U0],
> Y.sys::rearray(), //转换为Fortran数组
> IDO=1, TOL=sqrt(IMSL::MACH(4)), HINIT=0.01*TOL, NSTEP=10, T=0, TEND=0,
> i=0,(i<=NSTEP).while{
> TEND=(i+1)/NSTEP,
> if[i==NSTEP, IDO=3],
> IMSL::MOLCH ( &IDO, HFor("FCNUT"), HFor("FCNBC"), NPDES, &T, TEND, XBREAK, TOL,HINIT, Y),
> i++
> },
> Y.sys::rearray(), //转换为Forcal数组
> printff{"\r\n"},
> i=0,(i<NX).while{
> printff{"{1,i,-5}: {2,r,-16.8}",i+1, Y.XSLSF::getrai[i,0]},
> if{floor[(i+1)/10]*10==i+1, printff{"\r\n"}},
> i++
> },
> delete[XBREAK],delete[Y]
> };

[例子2] 线性标准双曲偏微分方程utt=uxx，也就是摆线方程，这是第一类偏微分方程系统。定义一个独立变量ut=v，那么vt=uxx就是这个系统的第二个方程。取初始值u(x,0)=sin(px)和ut(x,0)=v(x,0)=0，而绳子的另一端固定，从而u(0,t)=u(1,t)=v(0,t)=v(1,t)=0。这个方程的精确解为u(x,t)=sin(px)cos(pt)，在这个方程里，取t的范围为0<t<=2，步长为0.01，则在200步后可获得结果。 若取步长为0.001，则在2000步后可获得更优结果。

Forcal代码：

> FCNUT(x,t,u1,u2,ux1,ux2,uxx1,uxx2,ut1,ut2)=
> {
> ut1=u2,
> ut2=uxx1
> };
> FCNBC(x,t,α1,α2,β1,β2,γ1,γ2)=
> {
> α1=1, β1=0, γ1=0,
> α2=1, β2=0, γ2=0
> };
> main(:PARAM,XBREAK,Y,NPDES,IDO,NX,TOL,HINIT,NSTEP,T,TEND,ERRU,PI,i,j,t)=
> {
> IMSL::ClearImslErr(),
> NPDES=2, NX=10, PI=3.141592653589793,
> PARAM=new[rtoi(real_s),rtoi(50)],
> XBREAK=new[rtoi(real_s),rtoi(NX)],
> Y=new[rtoi(real_s),rtoi(2*NX),rtoi(NPDES)],
> IMSL::UMAG_DE[0,0,PARAM], //获得参数设置
> PARAM.XSLSF::setrai[16,1], //表示要传递导数
> IMSL::UMAG_DE[1,0,PARAM], //重新设置参数
> i=0,(i<NX).while{
> t=i/(NX-1),
> XBREAK.XSLSF::setrai[i,t],
> Y.XSLSF::setrai[i,0,sin(PI*t)], Y.XSLSF::setrai[i,1,0], //设置函数值
> Y.XSLSF::setrai[i+NX,0,PI*cos(PI*t)], Y.XSLSF::setrai[i+NX,1,0], //设置第一个导数值
> i++
> },
> IDO=1, TOL=0.1*sqrt(IMSL::MACH(4)), HINIT=0.01*TOL, NSTEP=200, T=0, TEND=0, ERRU=0,
> i=0,(i<=NSTEP).while{
> TEND=TEND+0.01,
> Y.sys::rearray(), //转换为Fortran数组
> if[i==NSTEP, IDO=3],
> IMSL::MOLCH ( &IDO, HFor("FCNUT"), HFor("FCNBC"), NPDES, &T, TEND, XBREAK, TOL,HINIT, Y),
> Y.sys::rearray(), //转换为Forcal数组
> j=0,(j<NX).while{
> t = Y.XSLSF::getrai[j,0] - sin(PI*XBREAK.XSLSF::getrai(j))*cos(PI*TEND),
> ERRU = max(ERRU,abs(t)),
> j++
> },
> i++
> },
> printff{"\r\n最终结果： 计算值 精确值\r\n",ERRU, PARAM.XSLSF::getrai[32]},
> j=0,(j<NX).while{
> printff{"\r\n{1,r,25.16}{2,r,25.16}",Y.XSLSF::getrai[j,0], sin(PI*XBREAK.XSLSF::getrai(j))*cos(PI*TEND)},
> j++
> },
> IMSL::UMAG_DE[0,0,PARAM], //获得参数设置
> printff{"\r\n\r\nMaximum error in u(x,t) : {1,r}\r\nMaximum step size used is : {2,r}",ERRU, PARAM.XSLSF::getrai[32]},
> IMSL::UMAG_DE[1,1,PARAM], //重新设置缺省参数
> delete[PARAM],delete[XBREAK],delete[Y]
> };

[返回页首][实数函数] IMSL::UMAG_DE(GetOrSet,IsDefault,Param)：获取或设置偏微分方程求解时的选项

GetOrSet：GetOrSet=0：获得选项值；GetOrSet!=0：设置选项值。
IsDefault：IsDefault=0：获取或设置选项的缺省值；IsDefault!=0：获取或设置选项的当前值。
Param：长度不小于50的一维数组，包含选项值，其意义参考IMSL说明。

运行错误：1：找不到相应的IMSL函数；2：数组不存在；3：不是一维数组 或数组太小。

3 如何封装IMSL库函数 [返回页首]

源代码下载：http://www.forcal.net/xiazai/forcal9/forcal9code.rar

3.1 用IMSL32.dll输出IMSL库函数 [返回页首]

使用支持Fortran90以上的编译器（如 Compaq Visual Fortran(CVF)或Intel Visual Fortran(IVF)等）创建动态库。用以下格式输出IMSL库函数：

> real(8) FUNCTION F_DSI(x)
> !DEC$ ATTRIBUTES DLLEXPORT::F_DSI
> USE Numerical_Libraries
> real(8) x
> F_DSI=DSI(x)
> END FUNCTION
> !!!!!!!!!!!!
> real(8) FUNCTION F_DDERIV(FCN,KORDER,X,BGSTEP,TOL)
> !DEC$ ATTRIBUTES DLLEXPORT::F_DDERIV
> EXTERNAL FCN
> integer(4) KORDER
> real(8) X,BGSTEP,TOL
> F_DDERIV=DDERIV(FCN,KORDER,X,BGSTEP,TOL)
> END FUNCTION

3.2 FcIMSL32W.dll中封装IMSL库函数的例子 [返回页首]

在FORCAL表达式中，函数的一般形式为：fun(x1,x2,...,xn)
与此相对应，FORCAL外部函数的通用形式为：double _stdcall fc_fun(fcINT mm,double *xx,void vFor);
参数说明：

mm：整数。mm指出自变量的个数。
xx：自变量数组的指针。若mm=-1，表示该函数没有自变量，若mm=0，表示有一个自变量，若mm=1，表示有2个自变量，... ...，依次类推。xx[i]即该函数的第i+1个自变量。
vFor：表达式句柄。通过该句柄可以获得表达式的详细信息。

根据函数所实现的功能，每一个自变量可以代表不同的意义，可以是一个双精度实数，也可以对它取整（(long)xx[i]）或使用部分字节表示其他意义，例如表示一个整数，或者表示一个指向任意对象的指针：数组、字符串、表达式或者其他任何一种对象类型。

由于IMSL函数是为编译型语言设计的，在设计为脚本调用时，需要检查传入的参数是否符合要求，即：检查参数是否有意义、验证函数句柄是否有效、验证数组指针是否有效等等。如果IMSL函数接受一个函数指针参数，还要为该函数设计一个专门的函数，并在该函数中转换为Forcal函数的调用，见下面的例子。

为了简便，以下均假定预先获得了IMSL库函数的句柄。

3.2.1 封装DSI函数：最简单的例子

> double _stdcall rfc_DSI(fcINT mm,double *xx,void *vFor)
> {
> return pDSI(xx); //pDSI是指向IMSL函数DSI的指针，下同
> }

3.2.2 封装DEI函数：需检查参数是否符合要求

> double _stdcall rfc_DEI(fcINT mm,double *xx,void *vFor)
> {
> static wchar_t ErrName[]=L"IMSL::EI";
> if(*xx==0.0)
> {
> if(!*FCERROR) SetRunErr(2,ErrName,1,0,vFor); return 0.0; //参数不能为0
> }
> return pDEI(xx);
> }

3.2.3 封装DDERIV函数：需检查参数是否符合要求、检查函数句柄是否符合要求 ；考虑到该函数可嵌套调用，故每次进入该函数，表达式的信息都要进行保存

> struct S_DERIV{ //定义保存表达式信息的结构
> void *rDERIV; //Forcal表达式句柄
> double *xDERIV; //Forcal表达式自变量数组
> S_DERIV *pNext; //指向下一个结构
> };
> S_DERIV *pMainS_DERIV=NULL; //指向保存表达式信息的结构的全局指针
> double _stdcall funDERIV(double *x) //为DDERIV专门定制的函数
> {
> *(pMainS_DERIV->xDERIV)=*x;
> return RealCal(pMainS_DERIV->rDERIV,pMainS_DERIV->xDERIV); //调用Forcal表达式（Forcal函数）
> }
> double _stdcall rfc_DDERIV(fcINT mm,double *xx,void *vFor)
> {
> static wchar_t ErrName[]=L"IMSL::DERIV";
> void *vvFor,*rPara;
> fcVOID nModule;
> fcINT nPara;
> long KORDER;
> S_DERIV *pS_DERIV;
> void *rDERIV;
> double *xDERIV,val;
> if(*FCERROR) return 0.0; //*FCERROR函数获取Forcal的运行错误，若出现过任何运行错误，不再进行计算；
> vvFor=(void *)*(fcVOID *)&xx[0];
> if(!GetFor(NULL,Key_RealFor,vvFor,nModule,rDERIV,rPara,nPara))
> {
> if(!*FCERROR) SetRunErr(2,ErrName,1,0,vFor); return 0.0; //指定的表达式不存在；
> }
> if(nPara!=0) {if(!*FCERROR) SetRunErr(2,ErrName,2,0,vFor); return 0.0;} //不是一元方程；
> xDERIV=(double *)rPara;
> KORDER=(long)xx[1];
> if(KORDER<1 || KORDER>3 || xx[3]<=0.0 || xx[4]<=0.0 || xx[4]>=1.0)
> {
> if(!*FCERROR) SetRunErr(2,ErrName,3,0,vFor); return 0.0; //参数不符合要求
> }
> pS_DERIV=new S_DERIV;
> if(!pS_DERIV)
> {
> if(!*FCERROR) SetRunErr(2,ErrName,4,0,vFor); return 0.0; //内存错误
> }
> pS_DERIV->rDERIV=rDERIV; pS_DERIV->xDERIV=xDERIV; pS_DERIV->pNext=pMainS_DERIV; pMainS_DERIV=pS_DERIV;
> val=pDDERIV(funDERIV,&KORDER,xx+2,xx+3,xx+4);
> pMainS_DERIV=pMainS_DERIV->pNext; delete pS_DERIV;
> return val;
> }

版权所有© Forcal程序设计 2002-2011，保留所有权利
E-mail: forcal@sina.com QQ：630715621
最近更新: 2011年05月03日