FORCAL数值计算扩展动态库FcIMSL V1.0 测试版
目 录
2 IMSL函数 |
|
错误输出 | |
IMSL::ClearImslErr() | 清空ISML错误输出文件“ImslErr.txt”的内容 |
IMSL::ERSET | 重置错误处理设置 |
IMSL::IERCD() | 返回信息错误代码,无错时返回0 |
IMSL::N1RTY(i) | 返回信息错误类型 |
指数积分和相关函数 | |
IMSL::EI(x) |
![]() |
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 | 用线性方法求解形如ut= f(x, t, u, ux, uxx)的偏微分方程,计算的结果是一个三次的Hermite多项式 |
IMSL::UMAG_DE | 获取或设置偏微分方程求解时的选项 |
公用函数 | |
IMSL::MACH(i) | 返回双精度实型常数 |
IMSL::IFNAN(x) | 判断x是否是NAN,NAN表示不是一个数 |
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公司的产品,在使用前须获得该公司的使用许可。
[返回页首][实数函数] 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):用线性方法求解形如ut= f(x, t, u, ux, 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,
t0),返回时获得ux(x,
tend)。
运行错误: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日