参数方程二阶求导,其中,先求出dy/dx 这个dy/dx 是代表关于t的一阶导吗?...参数方程二阶求导,其中,先求出dy/dx 这个dy/dx 是代表关于t的一阶导吗?如果是,那为什么求二阶导需要求关于x的二阶导呢?
pbNX64XN19
dy/dx 表示 y对x的一阶导数,此处它是 t 的函数,dy/dx = y '(t) / x '(t) = g(t) （记作 g(t) ） d²y/dx² 表示 y对x的二阶导数,也就是 dy/dx 对x 的导数,于是 d²y/dx² = g '(t) / x '(t) = [ y ''(t) x '(t) - y '(t) x ''(t) ] / [ x '(t) ] ³ (公式)
为您推荐：
其他类似问题
扫描下载二维码参数方程二阶导数的求法_百度文库
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
参数方程二阶导数的求法
上传于||文档简介
&&参​数​方​程​二​阶​导​数​的​求​法
阅读已结束，如果下载本文需要使用0下载券
想免费下载更多文档？
定制HR最喜欢的简历
你可能喜欢Forcal微分方程参数优化（拟合）
Forcal 微分方程参数优化（拟合）
例子1 数学模型：
&&&&& dy/dt=k*y*z+0.095*b*z
&&&&& dz/dt=-b*z-0.222*z
&&& 实验数据（ yi）：
&ti&&&&&&&
0.68&&& 3.
1.1&&&& 4.
1.63&&& 6.
2.07&&& 9.
2.67&&& 9.
3.09&&& 11.
3.64&&& 10.
4.65&&& 16.
5.1&&&& 18.
5.58&&& 21.
6.11&&& 25.
6.63&&& 26.
7.06&&& 26.
7.62&&& 27.
8.66&&& 35.
9.04&&& 35.5562035
9.63&&& 36.
&&& 已知z在t=0.01时的初值为5000，求k和b的最优值？
&&& Forcal代码1：
!using[&IMSL&,&fcopt&,&math&];
//使用命名空间
实验数据(::tyArray,max)=
&&& max=18,&&&&&&&&&&&&&&&&&&&
//实验数据组数
&&& tyArray=arrayinitns{max,2,
//存放实验数据ti,yi
&&&&&&& &0.01 2.
&&&&&&& 0.68 3.
&&&&&&& 1.1& 4.
&&&&&&& 1.63 6.
&&&&&&& 2.07 9.
&&&&&&& 2.67 9.
&&&&&&& 3.09 11.
&&&&&&& 3.64 10.
&&&&&&& 4.65 16.
&&&&&&& 5.1& 18.
&&&&&&& 5.58 21.
&&&&&&& 6.11 25.
&&&&&&& 6.63 26.
&&&&&&& 7.06 26.
&&&&&&& 7.62 27.
&&&&&&& 8.66 35.
&&&&&&& 9.04 35.5562035
&&&&&&& 9.63 36.&
&&& }.free()
f(t,y,z,dy,dz::k,b)=
&&& dy=k*y*z+0.095*b*z,
&&& dz=-b*z-0.222*z
目标函数(_k,_b : i,s,IDO,y,z,yy,t,tt : tyArray,max,k,b)=
&&& k=_k,b=_b,&&&&&&&&&&&&&&&&
//传递优化变量，函数f中要用到k,b
&&& i=0, s=0, IDO=1, t=tyArray[0,0], y=tyArray[0,1], z=5000,
&&& (++i&max).while{
&&&&&&& tt=tyArray[i,0], yy=tyArray[i,1],
&&&&&&& IVPRK[&IDO,HFor(&f&),&t,tt,1e-6,0,&y,&z],
&&&&&&& s=s+[y-yy]^2
&&& IVPRK[3,HFor(&f&),&t,tt+0.1,1e-6,0,&y,&z],
SimOpt[HFor(&目标函数&)];&&&&&
//求n维极值的单形调优法
&&& 结果（前面的数为最优参数，最后一个数是极小值。下同）：
1.297e-004
-1.575e-004 24.1
&&& Forcal代码2：优化（拟合）并绘图
!using[&IMSL&,&fcopt&,&math&,&fc2d&];
实验数据(::tyArray,max,tA,tyA)=
&&& max=18,
&&& tyArray=arrayinitns{max,2,
&&&&&&& &0.01 2.
&&&&&&& 0.68 3.
&&&&&&& 1.1& 4.
&&&&&&& 1.63 6.
&&&&&&& 2.07 9.
&&&&&&& 2.67 9.
&&&&&&& 3.09 11.
&&&&&&& 3.64 10.
&&&&&&& 4.65 16.
&&&&&&& 5.1& 18.
&&&&&&& 5.58 21.
&&&&&&& 6.11 25.
&&&&&&& 6.63 26.
&&&&&&& 7.06 26.
&&&&&&& 7.62 27.
&&&&&&& 8.66 35.
&&&&&&& 9.04 35.5562035
&&&&&&& 9.63 36.&
&&& }.free(),
&&& tA=tyArray(neg:0).free(),
tyA=tyArray(neg:1).free()&
//获得列向量
f(t,y,z,dy,dz::k,b)=
&&& dy=k*y*z+0.095*b*z,
&&& dz=-b*z-0.222*z
目标函数(_k,_b : i,s,IDO,y,z,yy,t,tt : tyArray,max,k,b)=
&&& k=_k,b=_b,
&&& i=0, s=0, IDO=1, t=tyArray[0,0], y=tyArray[0,1], z=5000,
&&& (++i&max).while{
&&&&&&& tt=tyArray[i,0], yy=tyArray[i,1],
&&&&&&& IVPRK[&IDO,HFor(&f&),&t,tt,1e-6,0,&y,&z],
&&&&&&& s=s+[y-yy]^2
&&& IVPRK[3,HFor(&f&),&t,tt+0.1,1e-6,0,&y,&z],
yi(x,y,n::tA,tyA,max)=
x=tA, y=tyA, n=max;&
//绘制实验数据点(ti,yi)的函数
yyi(tt : y,z,t : tyArray)=&
&&&&&&&&&&&&&&&
//绘制拟合曲线的函数
&&& t=tyArray[0,0], y=tyArray[0,1], z=5000,
&&& IVPRK[1,HFor(&f&),&t,tt,1e-6,0,&y,&z],
&&& IVPRK[3,HFor(&f&),&t,tt+0.1,1e-6,0,y,z],
main(::k,b)=
SimOpt[HFor(&目标函数&), optstarter,&k,&b,0],
Plot{Ix : 0,10,
&&&& Igarray : HFor(&yi&), Acolor,Vred,
Adot, Athick,8, //绘制函数yi的图形，将绘制散点图
&&&& Iufun : HFor(&yyi&), Acolor,Vblue, Adots,100
&& && //绘制函数yyi的图形
&&& 结果：
例子2 数学模型：三房室给药后，药物在房室1和2，2和3间转运，并经房室3排出。假定转运均符合一级速率模型，则各房室药物量（X1，X2，X3）与时间t的关系如下：
&&&&&&& dX1/dt=-K12*X1+K21*X2
&&&&&&& dX2/dt=(-K21-K23)*X2+K12*X1+K32*X3
&&&&&&& dX3/dt=(-K32-K30)*X3+K23*X2
&&& 已知初值t=0时：X1=100, X2=X3=0 。
&&& 实验值：
&&&&&&& t=0: X1=100, X2=0， X3=0
&&&&&&& t=1: X1=90,& X2=8， X3=2
&&&&&&& t=3: X1=73,& X2=19，X3=7
&&&&&&& t=5: X1=60,& X2=14，X3=23
&&& 求参数K12，K21，K23，K32，和K30。
&&& 限制条件：K12，K21，K23，K32，和K30均大于零；X1＋X2＋X3小于等于100（总量不能超过给药量）。
&&& 拟合可以用非线形最小二乘法，即求取minΣ（y实验值-y真实值）2时的K值。
&&& 当然这个问题的结构相对简单，可以求出各方程式，然后进行拟合。但是对于更复杂的，如十房室以上，积分式将非常复杂，手工求方程积分式几近不可能。
Forcal代码1：
!using[&IMSL&,&fcopt&,&math&,&sys&];
实验数据(::tArray,max)={&&& max=4,&&& tArray=arrayinit{2,max,4
&&&&&&& 0, 100, 0,& 0,&&&&&&& 1, 90,& 8,& 2,&&&&&&& 3, 73,& 19, 7,&&&&&&& 5, 60,& 14, 23&&& }.free()};
f(t,X1,X2,X3,dX1,dX2,dX3::k12,k21,k23,k32,k30)={&&& dX1=-k12*X1+k21*X2,&&& dX2=(-k21-k23)*X2+k12*X1+k32*X3,
&&& dX3=(-k32-k30)*X3+k23*X2};J(_k12,_k21,_k23,_k32,_k30 :
IDO,t1,t2,X1,X2,X3,X11,X22,X33,i,s : k12,k21,k23,k32,k30,tArray,max)=&&&
//目标函数{&&&
k12=_k12, k21=_k21, k23=_k23, k32=_k32, k30=_k30,&&& tArray.GA(0,&t1,&X1,&X2,&X3),&&&
IDO=1, s=0, i=0, (++i&max).while{&&&&&&&
tArray.GA(i*4,&t2,&X11,&X22,&X33),&&&&&&&
IVPRK[&IDO,HFor(&f&),&t1,t2,1e-6,0,&X1,&X2,&X3],&&&&&&& s=s+(X1-X11)^2+(X2-X22)^2+(X3-X33)^2&&& },&&&
IVPRK[3,HFor(&f&),&t1,t2+0.1,1e-6,0,&X1,&X2,&X3],
&&& s};S(_k12,_k21,_k23,_k32,_k30 : IDO,t1,t2,X1,X2,X3,i,s :
k12,k21,k23,k32,k30,tArray,max)=&&&
//约束函数{&&&
k12=_k12, k21=_k21, k23=_k23, k32=_k32, k30=_k30,&&& tArray.GA(0,&t1,&X1,&X2,&X3),&&&
IDO=1, s=1, i=0, (++i&max).while{&&&&&&&
IVPRK[&IDO,HFor(&f&),&t1,t2,1e-6,0,&X1,&X2,&X3],&&&&&&& s=s & (X1+X2+X3&=100),
&&&& //计算约束条件&&&&&&& i++&&& },&&&
IVPRK[3,HFor(&f&),&t1,t2+0.1,1e-6,0,&X1,&X2,&X3],
&&& s};main(:k12,k21,k23,k32,k30)={&&& k12=0.1, k21=0.1, k23=0.1, k32=0.1, k30=0.1,&&&//参数初值
&&& //获得一组符合约束的初值
&&& ROptIni[HFor(&J&),HFor(&S&),optmax,1000,
optstarter : &k12, &k21, &k23, &k32, &k30, 0, optrange
: 0,1e10, 0,1e10, 0,1e10, 0,1e10, 0,1e10],
&&& //优化
&&& RGOpt[HFor(&J&),HFor(&S&),optmax,2000, optmode,1,
optstarter : &k12, &k21, &k23, &k32, &k30, 0, optrange
: 0,1e10, 0,1e10, 0,1e10, 0,1e10, 0,1e10]};
&&& 结果：
8.57e-002 0..54e-003
&&& Forcal代码2：优化（拟合）并绘图
!using[&IMSL&,&fcopt&,&math&,&sys&,&fc2d&];
实验数据(::tArray,max,tX,tX1,tX2,tX3)={&&& max=4,&&& tArray=arrayinit{2,max,4
&&&&&&& 0, 100, 0,& 0,&&&&&&& 1, 90,& 8,& 2,&&&&&&& 3, 73,& 19, 7,&&&&&&& 5, 60,& 14, 23&&& }.free(),&&&
tX=tArray(neg:0).free(),
tX1=tArray(neg:1).free(),
tX2=tArray(neg:2).free(),
tX3=tArray(neg:3).free()
f(t,X1,X2,X3,dX1,dX2,dX3::k12,k21,k23,k32,k30)={&&& dX1=-k12*X1+k21*X2,&&& dX2=(-k21-k23)*X2+k12*X1+k32*X3,
&&& dX3=(-k32-k30)*X3+k23*X2};J(_k12,_k21,_k23,_k32,_k30 :
IDO,t1,t2,X1,X2,X3,X11,X22,X33,i,s : k12,k21,k23,k32,k30,tArray,max)={&&&
k12=_k12, k21=_k21, k23=_k23, k32=_k32, k30=_k30,&&& tArray.GA(0,&t1,&X1,&X2,&X3),&&&
IDO=1, s=0, i=0, (++i&max).while{&&&&&&&
tArray.GA(i*4,&t2,&X11,&X22,&X33),&&&&&&&
IVPRK[&IDO,HFor(&f&),&t1,t2,1e-6,0,&X1,&X2,&X3],&&&&&&& s=s+(X1-X11)^2+(X2-X22)^2+(X3-X33)^2&&& },&&&
IVPRK[3,HFor(&f&),&t1,t2+0.1,1e-6,0,&X1,&X2,&X3],
&&& s};S(_k12,_k21,_k23,_k32,_k30 : IDO,t1,t2,X1,X2,X3,i,s :
k12,k21,k23,k32,k30,tArray,max)={&&&
k12=_k12, k21=_k21, k23=_k23, k32=_k32, k30=_k30,&&& tArray.GA(0,&t1,&X1,&X2,&X3),&&&
IDO=1, s=1, i=0, (++i&max).while{&&&&&&&
IVPRK[&IDO,HFor(&f&),&t1,t2,1e-6,0,&X1,&X2,&X3],&&&&&&& s=s & (X1+X2+X3&=100),&&&&&&& i++&&& },&&&
IVPRK[3,HFor(&f&),&t1,t2+0.1,1e-6,0,&X1,&X2,&X3],
&&& s};X1(x,y,n::tX,tX1,max)=
//绘制实验数据点(t,X1)的函数
X2(x,y,n::tX,tX2,max)=
//绘制实验数据点(t,X2)的函数
X3(x,y,n::tX,tX3,max)=
//绘制实验数据点(t,X3)的函数
XX1(tt : X1,X2,X3,t : tArray)=&
&&&&&&&&&&&
//绘制拟合曲线X1的函数
&&& tArray.GA(0,&t,&X1,&X2,&X3),
&&& IVPRK[1,HFor(&f&),&t,tt,1e-6,0,&X1,&X2,&X3],
&&& IVPRK[3,HFor(&f&),&t,tt+0.1,1e-6,0,X1,X2,X3],
XX2(tt : X1,X2,X3,t : tArray)=&
&&&&&&&&&&&
//绘制拟合曲线X2的函数
&&& tArray.GA(0,&t,&X1,&X2,&X3),
&&& IVPRK[1,HFor(&f&),&t,tt,1e-6,0,&X1,&X2,&X3],
&&& IVPRK[3,HFor(&f&),&t,tt+0.1,1e-6,0,X1,X2,X3],
XX3(tt : X1,X2,X3,t : tArray)=&
&&&&&&&&&&&
//绘制拟合曲线X3的函数
&&& tArray.GA(0,&t,&X1,&X2,&X3),
&&& IVPRK[1,HFor(&f&),&t,tt,1e-6,0,&X1,&X2,&X3],
&&& IVPRK[3,HFor(&f&),&t,tt+0.1,1e-6,0,X1,X2,X3],
XXX(tt : X1,X2,X3,t : tArray)=&
&&&&&&&&&&&
//绘制拟合曲线X1+X2+X3的函数
&&& tArray.GA(0,&t,&X1,&X2,&X3),
&&& IVPRK[1,HFor(&f&),&t,tt,1e-6,0,&X1,&X2,&X3],
&&& IVPRK[3,HFor(&f&),&t,tt+0.1,1e-6,0,X1,X2,X3],
&&& X1+X2+X3
main(::k12,k21,k23,k32,k30)={&&& k12=0.1, k21=0.1, k23=0.1, k32=0.1, k30=0.1,
&&& ROptIni[HFor(&J&),HFor(&S&),optmax,1000,
optstarter : &k12, &k21, &k23, &k32, &k30, 0, optrange
: 0,1e10, 0,1e10, 0,1e10, 0,1e10, 0,1e10],&&& RGOpt[HFor(&J&),HFor(&S&),optmax,2000, optmode,1,
optstarter : &k12, &k21, &k23, &k32, &k30, 0, optrange
: 0,1e10, 0,1e10, 0,1e10, 0,1e10, 0,1e10],&&&
Plot{Ix : 0,5, Iy : 0,110, Iydynamic :
&&&& Igarray : HFor(&X1&), Acolor,Vred,
Adot, Athick,8, Alegend,&X1实验数据&,&&
//绘制函数X1的图形，将绘制散点图
&&&& Igarray : HFor(&X2&), Acolor,Vblue,
Adot, Athick,8, Alegend,&X2实验数据&,&
//绘制函数X2的图形，将绘制散点图
&&&& Igarray : HFor(&X3&), Acolor,Vgreen,
Adot, Athick,8, Alegend,&X3实验数据&,
//绘制函数X3的图形，将绘制散点图
&&&& Iufun : HFor(&XX1&), Acolor,Vred, Adots,100, Alegend,&X1模拟数据&,
//绘制函数XX1的图形
&&&& Iufun : HFor(&XX2&), Acolor,Vblue, Adots,100, Alegend,&X2模拟数据&,
&& &&& //绘制函数XX2的图形
&&&& Iufun : HFor(&XX3&), Acolor,Vgreen, Adots,100, Alegend,&X3模拟数据&,
&& && //绘制函数XX3的图形
&&&& Iufun : HFor(&XXX&), Acolor,RGB(255,128,0), Adots,100,
Athick,8, Alegend,&X1+X2+X3模拟数据&&
//绘制函数XXX的图形
&&& 结果：
例子3 数学模型：需要拟合的微分方程组为:
df1/dt=a1*(x-f2)^a2,
&&&&&&& df2/dt=a3*df1/dt-a4*f2^a5
其中：x=0.3。拟合参数为a1,a2,a3,a4,a5。试验获取的数据为（t,f1）；t=0时，f2=0.15。
f10&& , 0&&&&&&& ,10& , 0.002174 ,20& , 0.002663 ,30& , 0.002934 ,40& , 0.003113 ,50& , 0.003248 ,60& , 0.003354 ,70& , 0.003442 ,80& , 0.003514 ,90& , 0.003578 ,100 , 0.003635 ,110 , 0.003686 ,120 , 0.00373 ,130 , 0.003774 ,140 , 0.003813 ,150 , 0.003847 ,160 , 0.003882 ,170 , 0.003913 ,180 , 0.003942 ,190 , 0.00397& ,200 , 0.003996 ,210 , 0.00402& ,220 , 0.004044 ,230 , 0.004067 ,240 , 0.004087 ,250 , 0.004107 ,260 , 0.004126 ,270 , 0.004143 ,280 , 0.004159 ,290 , 0.004174 ,300 , 0.00419& ,310 , 0.004207 ,320 , 0.00422& ,330 , 0.004235 ,340 , 0.004248 ,350 , 0.004263 ,360 , 0.004276 ,370 , 0.004287 ,380 , 0.004301 ,390 , 0.00431& ,400 , 0.004323 ,410 , 0.004333 ,420 , 0.004344 ,430 , 0.004352 ,440 , 0.004362 ,450 , 0.004375 ,460 , 0.004383 ,470 , 0.00439& ,480 , 0.004399 ,490 , 0.004408 ,500 , 0.004416 ,510 , 0.004424 ,520 , 0.00443& ,530 , 0.00444& ,540 , 0.004449 ,550 , 0.004453 ,560 , 0.004458 ,570 , 0.004468 ,580 , 0.004474 ,590 , 0.004483 ,600 , 0.004488 ,610 , 0.004494 ,620 , 0.004499 ,630 , 0.004505 ,640 , 0.00451& ,650 , 0.004516 ,660 , 0.004522 ,670 , 0.004527 ,680 , 0.004532 ,690 , 0.004537 ,700 , 0.004541 ,710 , 0.004548 ,720 , 0.004552 ,730 , 0.004557 ,740 , 0.004561 ,750 , 0.004566 ,760 , 0.004569 ,770 , 0.004575 ,780 , 0.004579 ,790 , 0.004584 ,800 , 0.004589 ,810 , 0.004594 ,820 , 0.004596 ,830 , 0.004601 ,840 , 0.004604 ,850 , 0.004609 ,860 , 0.004611 ,870 , 0.004614 ,880 , 0.00462& ,890 , 0.004622 ,900 , 0.004626 ,910 , 0.00463& ,920 , 0.004632 ,930 , 0.004636 ,940 , 0.004638 ,950 , 0.004641 ,960 , 0.004647 ,970 , 0.004649 ,980 , 0.004653 ,990 , 0.004655 ,658
&&& Forcal代码1：
!using[&IMSL&,&fcopt&,&math&];
init(::tArray,max)={&&& max=101,&&& tArray=arrayinit{2,max,2
&&&&&&& 0 ,& 0 ,&&&&&&& 10 , 0.002174 ,&&&&&&& 20 , 0.002663 ,&&&&&&& 30 , 0.002934 ,&&&&&&& 40 , 0.003113 ,&&&&&&& 50 , 0.003248 ,&&&&&&& 60 , 0.003354 ,&&&&&&& 70 , 0.003442 ,&&&&&&& 80 , 0.003514 ,&&&&&&& 90 , 0.003578 ,&&&&&&& 100 , 0.003635 ,&&&&&&& 110 , 0.003686 ,&&&&&&& 120 , 0.00373 ,&&&&&&& 130 , 0.003774 ,&&&&&&& 140 , 0.003813 ,&&&&&&& 150 , 0.003847 ,&&&&&&& 160 , 0.003882 ,&&&&&&& 170 , 0.003913 ,&&&&&&& 180 , 0.003942 ,&&&&&&& 190 , 0.00397 ,&&&&&&& 200 , 0.003996 ,&&&&&&& 210 , 0.00402 ,&&&&&&& 220 , 0.004044 ,&&&&&&& 230 , 0.004067 ,&&&&&&& 240 , 0.004087 ,&&&&&&& 250 , 0.004107 ,&&&&&&& 260 , 0.004126 ,&&&&&&& 270 , 0.004143 ,&&&&&&& 280 , 0.004159 ,&&&&&&& 290 , 0.004174 ,&&&&&&& 300 , 0.00419 ,&&&&&&& 310 , 0.004207 ,&&&&&&& 320 , 0.00422 ,&&&&&&& 330 , 0.004235 ,&&&&&&& 340 , 0.004248 ,&&&&&&& 350 , 0.004263 ,&&&&&&& 360 , 0.004276 ,&&&&&&& 370 , 0.004287 ,&&&&&&& 380 , 0.004301 ,&&&&&&& 390 , 0.00431 ,&&&&&&& 400 , 0.004323 ,&&&&&&& 410 , 0.004333 ,&&&&&&& 420 , 0.004344 ,&&&&&&& 430 , 0.004352 ,&&&&&&& 440 , 0.004362 ,&&&&&&& 450 , 0.004375 ,&&&&&&& 460 , 0.004383 ,&&&&&&& 470 , 0.00439 ,&&&&&&& 480 , 0.004399 ,&&&&&&& 490 , 0.004408 ,&&&&&&& 500 , 0.004416 ,&&&&&&& 510 , 0.004424 ,&&&&&&& 520 , 0.00443 ,&&&&&&& 530 , 0.00444 ,&&&&&&& 540 , 0.004449 ,&&&&&&& 550 , 0.004453 ,&&&&&&& 560 , 0.004458 ,&&&&&&& 570 , 0.004468 ,&&&&&&& 580 , 0.004474 ,&&&&&&& 590 , 0.004483 ,&&&&&&& 600 , 0.004488 ,&&&&&&& 610 , 0.004494 ,&&&&&&& 620 , 0.004499 ,&&&&&&& 630 , 0.004505 ,&&&&&&& 640 , 0.00451 ,&&&&&&& 650 , 0.004516 ,&&&&&&& 660 , 0.004522 ,&&&&&&& 670 , 0.004527 ,&&&&&&& 680 , 0.004532 ,&&&&&&& 690 , 0.004537 ,&&&&&&& 700 , 0.004541 ,&&&&&&& 710 , 0.004548 ,&&&&&&& 720 , 0.004552 ,&&&&&&& 730 , 0.004557 ,&&&&&&& 740 , 0.004561 ,&&&&&&& 750 , 0.004566 ,&&&&&&& 760 , 0.004569 ,&&&&&&& 770 , 0.004575 ,&&&&&&& 780 , 0.004579 ,&&&&&&& 790 , 0.004584 ,&&&&&&& 800 , 0.004589 ,&&&&&&& 810 , 0.004594 ,&&&&&&& 820 , 0.004596 ,&&&&&&& 830 , 0.004601 ,&&&&&&& 840 , 0.004604 ,&&&&&&& 850 , 0.004609 ,&&&&&&& 860 , 0.004611 ,&&&&&&& 870 , 0.004614 ,&&&&&&& 880 , 0.00462 ,&&&&&&& 890 , 0.004622 ,&&&&&&& 900 , 0.004626 ,&&&&&&& 910 , 0.00463 ,&&&&&&& 920 , 0.004632 ,&&&&&&& 930 , 0.004636 ,&&&&&&& 940 , 0.004638 ,&&&&&&& 950 , 0.004641 ,&&&&&&& 960 , 0.004647 ,&&&&&&& 970 , 0.004649 ,&&&&&&& 980 , 0.004653 ,&&&&&&& 990 , 0.004655 ,&&&&&&& 658
&&& }.free()};
f(t,f1,f2,df1,df2 : x : a1,a2,a3,a4,a5)={&&& x=0.3,&&& df1=a1*(x-f2)^a2,&&& df2=a3*df1-a4*f2^a5};J(_a1,_a2,_a3,_a4,_a5 :
IDO,i,s,f1,f2,f11,t1,t2 : tArray,max,a1,a2,a3,a4,a5)={&&& a1=_a1, a2=_a2, a3=_a3, a4=_a4, a5=_a5,&&&
t1=tArray(0,0), f1=tArray(0,1), f2=0.15,&&& IDO=1, i=0,
s=0,&&& (++i&max).while{&&&&&&&
t2=tArray(i,0), f11=tArray(i,1),&&&&&&&
IVPRK[&IDO,HFor(&f&),&t1,t2,1e-6,0,&f1,&f2],&&&&&&& s=s+[f1-f11]^2&&& },&&&
IVPRK[3,HFor(&f&),&t1,t2+0.1,1e-6,0,&f1,&f2],
&&& s};ClearImslErr(),&&&&&&&&&&&&
//清空IMSL错误输出
ERSET(0,0,0),&&&&&&&&&&&&&&
//关闭IMSL所有警告
Opt[HFor(&J&)],&&&&&
//Opt函数全局优化
ERSET(0,2,2), ERSET(0,1,0);
//恢复IMSL警告
&&& 结果（没有多次运行，或许还有更优解）：
8.728 18.71 -4.431e-006
4.509e-015 1.018e-009
&&& Forcal代码2：优化（拟合）并绘图
!using[&IMSL&,&fcopt&,&math&,&fc2d&];
init(::tArray,max,tA,f1A)={&&& max=101,&&& tArray=arrayinit{2,max,2
&&&&&&& 0 ,& 0 ,&&&&&&& 10 , 0.002174 ,&&&&&&& 20 , 0.002663 ,&&&&&&& 30 , 0.002934 ,&&&&&&& 40 , 0.003113 ,&&&&&&& 50 , 0.003248 ,&&&&&&& 60 , 0.003354 ,&&&&&&& 70 , 0.003442 ,&&&&&&& 80 , 0.003514 ,&&&&&&& 90 , 0.003578 ,&&&&&&& 100 , 0.003635 ,&&&&&&& 110 , 0.003686 ,&&&&&&& 120 , 0.00373 ,&&&&&&& 130 , 0.003774 ,&&&&&&& 140 , 0.003813 ,&&&&&&& 150 , 0.003847 ,&&&&&&& 160 , 0.003882 ,&&&&&&& 170 , 0.003913 ,&&&&&&& 180 , 0.003942 ,&&&&&&& 190 , 0.00397 ,&&&&&&& 200 , 0.003996 ,&&&&&&& 210 , 0.00402 ,&&&&&&& 220 , 0.004044 ,&&&&&&& 230 , 0.004067 ,&&&&&&& 240 , 0.004087 ,&&&&&&& 250 , 0.004107 ,&&&&&&& 260 , 0.004126 ,&&&&&&& 270 , 0.004143 ,&&&&&&& 280 , 0.004159 ,&&&&&&& 290 , 0.004174 ,&&&&&&& 300 , 0.00419 ,&&&&&&& 310 , 0.004207 ,&&&&&&& 320 , 0.00422 ,&&&&&&& 330 , 0.004235 ,&&&&&&& 340 , 0.004248 ,&&&&&&& 350 , 0.004263 ,&&&&&&& 360 , 0.004276 ,&&&&&&& 370 , 0.004287 ,&&&&&&& 380 , 0.004301 ,&&&&&&& 390 , 0.00431 ,&&&&&&& 400 , 0.004323 ,&&&&&&& 410 , 0.004333 ,&&&&&&& 420 , 0.004344 ,&&&&&&& 430 , 0.004352 ,&&&&&&& 440 , 0.004362 ,&&&&&&& 450 , 0.004375 ,&&&&&&& 460 , 0.004383 ,&&&&&&& 470 , 0.00439 ,&&&&&&& 480 , 0.004399 ,&&&&&&& 490 , 0.004408 ,&&&&&&& 500 , 0.004416 ,&&&&&&& 510 , 0.004424 ,&&&&&&& 520 , 0.00443 ,&&&&&&& 530 , 0.00444 ,&&&&&&& 540 , 0.004449 ,&&&&&&& 550 , 0.004453 ,&&&&&&& 560 , 0.004458 ,&&&&&&& 570 , 0.004468 ,&&&&&&& 580 , 0.004474 ,&&&&&&& 590 , 0.004483 ,&&&&&&& 600 , 0.004488 ,&&&&&&& 610 , 0.004494 ,&&&&&&& 620 , 0.004499 ,&&&&&&& 630 , 0.004505 ,&&&&&&& 640 , 0.00451 ,&&&&&&& 650 , 0.004516 ,&&&&&&& 660 , 0.004522 ,&&&&&&& 670 , 0.004527 ,&&&&&&& 680 , 0.004532 ,&&&&&&& 690 , 0.004537 ,&&&&&&& 700 , 0.004541 ,&&&&&&& 710 , 0.004548 ,&&&&&&& 720 , 0.004552 ,&&&&&&& 730 , 0.004557 ,&&&&&&& 740 , 0.004561 ,&&&&&&& 750 , 0.004566 ,&&&&&&& 760 , 0.004569 ,&&&&&&& 770 , 0.004575 ,&&&&&&& 780 , 0.004579 ,&&&&&&& 790 , 0.004584 ,&&&&&&& 800 , 0.004589 ,&&&&&&& 810 , 0.004594 ,&&&&&&& 820 , 0.004596 ,&&&&&&& 830 , 0.004601 ,&&&&&&& 840 , 0.004604 ,&&&&&&& 850 , 0.004609 ,&&&&&&& 860 , 0.004611 ,&&&&&&& 870 , 0.004614 ,&&&&&&& 880 , 0.00462 ,&&&&&&& 890 , 0.004622 ,&&&&&&& 900 , 0.004626 ,&&&&&&& 910 , 0.00463 ,&&&&&&& 920 , 0.004632 ,&&&&&&& 930 , 0.004636 ,&&&&&&& 940 , 0.004638 ,&&&&&&& 950 , 0.004641 ,&&&&&&& 960 , 0.004647 ,&&&&&&& 970 , 0.004649 ,&&&&&&& 980 , 0.004653 ,&&&&&&& 990 , 0.004655 ,&&&&&&& 658
&&& }.free(),&&& tA=tArray(neg:0).free(),
f1A=tArray(neg:1).free()
f(t,f1,f2,df1,df2 : x : a1,a2,a3,a4,a5)={&&& x=0.3,&&& df1=a1*(x-f2)^a2,&&& df2=a3*df1-a4*f2^a5};J(_a1,_a2,_a3,_a4,_a5 :
IDO,i,s,f1,f2,f11,t1,t2 : tArray,max,a1,a2,a3,a4,a5)={&&& a1=_a1, a2=_a2, a3=_a3, a4=_a4, a5=_a5,&&&
t1=tArray(0,0), f1=tArray(0,1), f2=0.15,&&& IDO=1, i=0,
s=0,&&& (++i&max).while{&&&&&&&
t2=tArray(i,0), f11=tArray(i,1),&&&&&&&
IVPRK[&IDO,HFor(&f&),&t1,t2,1e-6,0,&f1,&f2],&&&&&&& s=s+[f1-f11]^2&&& },&&&
IVPRK[3,HFor(&f&),&t1,t2+0.1,1e-6,0,&f1,&f2],
&&& s};f1i(x,y,n::tA,f1A,max)=
x=tA, y=f1A, n=max;&
//绘制实验数据点(t,f1)的函数
ff1i(tt : f1,f2,t : tArray)=&
&&&&&&&&&&&&&&
//绘制拟合曲线的函数
&&& t=tArray[0,0], f1=tArray[0,1], f2=0.15,
&&& IVPRK[1,HFor(&f&),&t,tt,1e-6,0,&f1,&f2],
&&& IVPRK[3,HFor(&f&),&t,tt+0.1,1e-6,0,f1,f2],
main(::a1,a2,a3,a4,a5)=
&&& ClearImslErr(),&&&&&&&&&&&&
//清空IMSL错误输出
ERSET(0,0,0),&&&&&&&&&&&&&&
//关闭IMSL所有警告
Opt[HFor(&J&), optstarter,&a1,&a2,&a3,&a4,&a5,0],&
//Opt函数全局优化
Plot{Iclear, Ix : 0,1000,
&&&& Igarray : HFor(&f1i&), Acolor,Vred,
Adot,&&&&&&&&
//绘制函数f1i的图形，将绘制散点图
&&&& Iufun : HFor(&ff1i&), Acolor,Vblue, Adots,100
&& &//绘制函数ff1i的图形
&&& ERSET(0,2,2), ERSET(0,1,0)
&//恢复IMSL警告
&&& 结果：
例子4 数学模型：
有关于y1,y2,y3,y4的微分方程：
k1 = 3.5*10^4;
k2 = 1.0*10^3;
k3 = a*ka+b*kb+c*
r1 = k1*y1*y2;
r2 = k2*y1*y3;
r3 = k3*y1*y4;
dy1/dt = -r1-r2-r3;
dy2/dt = -r1;
dy3/dt = -r2+r1;
dy4/dt = -r3+r2;
其中ka，kb，kc为拟合参数；常数a、b、c为已知数据；k1,k2,k3,r1,r2,r3为中间变量；t=0时，y1=30.0*10^-6;
y2=10.2*10^-6; y3=4.1*10^-6; y4=6.7。
有4组数据：
当[a b c]=[100 5 20]时
t=[0, 10, 20, 30, 40]
y1=[30.0*10^-6, 17.5*10^-6, 15*10^-6, 14*10^-6, 13.5*10^-6 ]
当[a b c]=[120 15 25]时
t=[0 10 20 30 40]
y1=[30.0*10^-6, 12.5*10^-6, 10*10^-6, 8*10^-6, 7*10^-6 ]
当[a b c]=[130 25 30]时
t=[0 10 20 30 40]
y(1)=[30.0*10^-6, 11.5*10^-6, 9*10^-6, 7.5*10^-6, 6*10^-6 ]
当[a b c]=[140 30 35]时
t=[0 10 20 30 40]
y(1)=[30.0*10^-6, 10.5*10^-6, 8.5*10^-6, 6.5*10^-6, 5*10^-6 ]
求ka，kb，kc？
分析：观察式子 k3 = a*ka+b*kb+c*kc
，当[a b c]实验数据少于3组时，将有无穷组最优的ka，kb，kc。现在有4组数据，故只有一组最优解。
Forcal代码：
!using[&IMSL&,&fcopt&,&math&];
实验数据(::tArray,y11Array,y12Array,y13Array,y14Array,max)=
&&& max=5,
&&& tArray=arrayinit[1,max : 0, 10, 20, 30, 40].free(),
&&& y11Array=arrayinit[1,max : 30.0e-6, 17.5e-6, 15e-6, 14e-6,
13.5e-6].free(),
&&& y12Array=arrayinit[1,max : 30.0e-6, 12.5e-6, 10e-6, 8e-6,
7e-6].free(),
&&& y13Array=arrayinit[1,max : 30.0e-6, 11.5e-6, 9e-6, 7.5e-6,
6e-6].free(),
&&& y14Array=arrayinit[1,max : 30.0e-6, 10.5e-6, 8.5e-6, 6.5e-6,
5e-6].free()
f(t,y1,y2,y3,y4,z1,z2,z3,z4 : k1,k2,k3,r1,r2,r3 : a,b,c,ka,kb,kc)=
&&& k1 = 3.5e4,
&&& k2 = 1.0e3,
&&& k3 = a*ka+b*kb+c*kc,
&&& r1 = k1*y1*y2,
&&& r2 = k2*y1*y3,
&&& r3 = k3*y1*y4,
&&& z1 = -r1-r2-r3,
&&& z2 = -r1,
&&& z3 = -r2+r1,
&&& z4 = -r3+r2
目标函数(kka,kkb,kkc : i,s,IDO,t,tt,y1,y2,y3,y4,z1 :
tArray,y11Array,y12Array,y13Array,y14Array,max,a,b,c,ka,kb,kc)=
&&& ka=kka, kb=kkb, kc=kkc,
&&& a=100, b=5, c=20, //[a
b c]=[100 5 20]
&&& i=0, IDO=1, t=tArray[0], y1=30.0e-6, y2=10.2e-6, y3=4.1e-6,
&&& (++i&max).while{
&&&&&&& tt=tArray[i], z1=y11Array[i],
&&&&&&& IVPRK[&IDO,HFor(&f&),&t,tt,1e-6,0,&y1,&y2,&y3,&y4],
&&&&&&& s=s+(y1-z1)^2
&&& IVPRK[3,HFor(&f&),&t,tt+0.1,1e-6,0,&y1,&y2,&y3,&y4],
&&& a=120, b=15, c=25, //[a
b c]=[120, 15, 25]
&&& i=0, IDO=1, t=tArray[0], y1=30.0e-6, y2=10.2e-6, y3=4.1e-6,
&&& (++i&max).while{
&&&&&&& tt=tArray[i], z1=y12Array[i],
&&&&&&& IVPRK[&IDO,HFor(&f&),&t,tt,1e-6,0,&y1,&y2,&y3,&y4],
&&&&&&& s=s+(y1-z1)^2
&&& IVPRK[3,HFor(&f&),&t,tt+0.1,1e-6,0,&y1,&y2,&y3,&y4],
&&& a=130, b=25, c=30, //[a
b c]=[130, 25, 30]
&&& i=0, IDO=1, t=tArray[0], y1=30.0e-6, y2=10.2e-6, y3=4.1e-6,
&&& (++i&max).while{
&&&&&&& tt=tArray[i], z1=y13Array[i],
&&&&&&& IVPRK[&IDO,HFor(&f&),&t,tt,1e-6,0,&y1,&y2,&y3,&y4],
&&&&&&& s=s+(y1-z1)^2
&&& IVPRK[3,HFor(&f&),&t,tt+0.1,1e-6,0,&y1,&y2,&y3,&y4],
&&& a=140, b=30, c=35, //[a
b c]=[140, 30, 35]
&&& i=0, IDO=1, t=tArray[0], y1=30.0e-6, y2=10.2e-6, y3=4.1e-6,
&&& (++i&max).while{
&&&&&&& tt=tArray[i], z1=y14Array[i],
&&&&&&& IVPRK[&IDO,HFor(&f&),&t,tt,1e-6,0,&y1,&y2,&y3,&y4],
&&&&&&& s=s+(y1-z1)^2
&&& IVPRK[3,HFor(&f&),&t,tt+0.1,1e-6,0,&y1,&y2,&y3,&y4],
&&& sqrt[s/(max*4-4)]
ClearImslErr(),&&&&&&&&&&&&
//清空IMSL错误输出
ERSET(0,0,0),&&&&&&&&&&&&&&
//关闭IMSL所有警告
Opt[HFor(&目标函数&), optrange,-1e5,1e5,-1e5,1e5,-1e5,1e5],
//Opt函数优化
ERSET(0,2,2), ERSET(0,1,0); //恢复IMSL警告
结果（ka,kb,kc，目标函数终值）：
3.708e-004 -5.95e-004 1.296e-006
验证（t，实验值，拟合值）：
当[a b c]=[100 5 20]时
10. 1.75e-005 1.7
20. 1.5e-005& 1.5
30. 1.4e-005& 1.4
40. 1.35e-005 1.2
当[a b c]=[120 15 25]时
10. 1.25e-005 1.4
20. 1.e-005&& 1.0
30. 8.e-006&& 8.2
40. 7.e-006&& 6.3
当[a b c]=[130 25 30]时
10. 1.15e-005 1.2
20. 9.e-006&& 8.7
30. 7.5e-006& 5.9
40. 6.e-006&& 4.0
当[a b c]=[140 30 35]时
10. 1.05e-005 1.3
20. 8.5e-006& 8.8
30. 6.5e-006& 6.0
40. 5.e-006&& 4.1
例子5 数学模型：
微分方程组为：
dn1/dt = a*n1*n2-b*n1;
dn2/dt = c*n1*n2;
其中a，b，c为拟合参数；t=0时，n1=n0; n2=n0。
2&&&& 3&&&& 4&&&& 5&&&&
3.5&& 2.8&& 2.4&& 1.5&& 0.9
分析：由于n0未知，故共有4个拟合参数：n0,a,b,c。
Forcal代码：
!using[&IMSL&,&fcopt&,&math&];
实验数据(::tArray,n1Array,max)=
&&& max=6,
&&& tArray=arrayinitns{max : &1 2 3 4 5 6&}.free(),
&&& n1Array=arrayinitns{max : &1 3.5 2.8 2.4 1.5 0.9&}.free()
f(t,n1,n2,dn1,dn2::a,b,c)=
&&& dn1=a*n1*n2-b*n1,
&&& dn2=c*n1*n2
目标函数(n0,aa,bb,cc : i,s,IDO,n1,n2,nn1,t,tt : tArray,n1Array,max,a,b,c)=
&&& a=aa, b=bb, c=cc,
&&& i=0, s=0, IDO=1, t=0, n1=n0, n2=n0,
&&& (i&max).while{
&&&&&&& tt=tArray[i], nn1=n1Array[i],
&&&&&&& IVPRK[&IDO,HFor(&f&),&t,tt,1e-6,0,&n1,&n2],
&&&&&&& s=s+[n1-nn1]^2,
&&&&&&& i++
&&& IVPRK[3,HFor(&f&),&t,tt+0.1,1e-6,0,&n1,&n2],
&&& sqrt[s/max]
ClearImslErr(),&&&&&&&&&&&&
//清空IMSL错误输出
ERSET(0,0,0),&&&&&&&&&&&&&&
//关闭IMSL所有警告
Opt[HFor(&目标函数&), optwaysimdeep, optwayconfra, optrange,-1e10,1e10,-1e10,1e10,-1e10,1e10,-1e10,1e10],
//Opt函数优化
ERSET(0,2,2), ERSET(0,1,0); //恢复IMSL警告
结果（n0,a,b,c，误差）：
4.411e-002
74.58 0.5 0.2884
验证（t，实验值，拟合值）：
1. 1.&& 1.00662
2. 3.5& 3.47867
3. 2.8& 2.96124
4. 2.4& 2.12753
5. 1.5& 1.50219
6. 0.9& 1.05743
版权所有& Forcal程序设计
，保留所有权利