常微分方程边值问题的数值解法_百度文库
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
常微分方程边值问题的数值解法
上传于||暂无简介
阅读已结束，如果下载本文需要使用2下载券
想免费下载本文？
定制HR最喜欢的简历
下载文档到电脑，查找使用更方便
还剩8页未读，继续阅读
定制HR最喜欢的简历
你可能喜欢查看: 8078|回复: 19|关注: 0
Matlab 微分方程组参数拟合求助！
<h1 style="color:# 麦片财富积分
新手, 积分 5, 距离下一级还需 45 积分
我的问题很类似前不久的一个帖子：“求助:Matlab常微分方程组参数的拟合，有点难度！” 由于是新新新手，（刚装了几天）真是看不明白。所以还是寄望于高手相助。我的问题如下：三房室模型给药后，药物在在房室1和2，2和3间转运， 并经房室3排出。假定转运均符合一级速率模型，则各房室药物量（X1，X2，X3）与时间t的关系如下：
dX1/dt=-K1,2*X1+K2,1*X2
dX2/dt=(-K2,1-K2,3)*X2+K1,2*X1+K3,2*X3
dX3/dt=(-K3,2-K3,0)*X3+K2,3*X2
已知初值(t=0)：X1=100, X2=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
求拟合的K1,2， K2,1， K2,3， K3,2，和K3,0。
限制条件：K1,2， K2,1， K2,3， K3,2，和K3,0 均大于零。X1＋X2＋X3 小于等于100（总量不能超过给药量）。
拟合可以用非线形最小二乘法，即求取minΣ(y实验值-y真实值)2时的k值。
先行谢过！
<h1 style="color:# 麦片财富积分
当然这个问题的结构相对简单，所以也可以求出各方程积分式，然后用EXCEL规划求解拟合。这是药物动力学中的方法。但是对于更复杂的系统，如十房室以上，积分式将非常复杂，手工求方程积分式几近不可能。
[ 本帖最后由 sunnoom 于
03:16 编辑 ]
<h1 style="color:#2 麦片财富积分
关注者： 4
用Forcal求解，代码：
//这里是代码窗口，请将Forcal代码写在下面
//若长时间没有结果，按Ctrl+Alt+q退出Forcal运行
i: OutVector(p:k,i)= k=FCDLen(p),printff{&\r\n&},i=0,(i&k).while{printff{&{1,r,14.6}&,get[p,i]},i++},printff
{&\r\n&};& & //输出一维数组
!using[&XSLSF&]; //使用命名空间XSLSF
f(t,X1,X2,X3,dX1,dX2,dX3::k12,k21,k23,k32,k30)={ //函数定义，连分式法对微分方程组积分一步函数pbs1中要用到
&&dX1=-k12*X1+k21*X2,
&&dX2=(-k21-k23)*X2+k12*X1+k32*X3,
&&dX3=(-k32-k30)*X3+k23*X2
};
t_i(hf,a,step,eps,t1,t2,x1,x2,x3:h,i)=&&//用于验证初值
{
& & h=(t2-t1)/step,
& & {& &pbs1[hf,t1,a,h,eps],&&//连分式法对微分方程组积分一步函数pbs1，hf为函数f的句柄
& && &&&t1=t1+h
& & }.until[abs(t1-t2)&h/2],
& & a.getra(0,&x1,&x2,&x3)
};
t_i_2(hf,a,step,eps,t1,t2,x_1,x_2,x_3:x1,x2,x3,h,i)=& & //用于计算目标函数
{
& & h=(t2-t1)/step,
& & {& &pbs1[hf,t1,a,h,eps],&&//连分式法对微分方程组积分一步函数pbs1，hf为函数f的句柄
& && &&&t1=t1+h
& & }.until[abs(t1-t2)&h/2],
& & a.getra(0,&x1,&x2,&x3),
& & (x1-x_1)^2+(x2-x_2)^2+(x3-x_3)^2
};
J(_k12,_k21,_k23,_k32,_k30:t1:hf,Array,step,eps,k12,k21,k23,k32,k30)={& & //目标函数定义
&&k12=_k12,k21=_k21,k23=_k23,k32=_k32,k30=_k30,
&&t1=0, Array.setra(0,100,0,0),
&&t_i_2(hf,Array,step,eps: &t1, 1 : 90,&&8,&&2)+
&&t_i_2(hf,Array,step,eps: &t1, 3 : 73, 19,&&7)+
&&t_i_2(hf,Array,step,eps: &t1, 5 : 60, 14, 23)
};
S(_k12,_k21,_k23,_k32,_k30,c0,c1,c2,d0,d1,d2,w0,w1,w2:t1,x1,x2,x3:hf,Array,step,eps,k12,k21,k23,k32,k30)=& & //约束条件定义
{
&&c0=0,& && && &c1=0,& && && &c2=0,
&&d0=100,& && & d1=100,& && & d2=100,
&&k12=_k12,k21=_k21,k23=_k23,k32=_k32,k30=_k30,
&&t1=0, Array.setra(0,100,0,0),
&&t_i(hf,Array,step,eps: &t1, 1 :&&&x1, &x2, &x3), w0=x1+x2+x3,
&&t_i(hf,Array,step,eps: &t1, 3 :&&&x1, &x2, &x3), w1=x1+x2+x3,
&&t_i(hf,Array,step,eps: &t1, 5 :&&&x1, &x2, &x3), w2=x1+x2+x3
};
main(:a,b,alpha,_eps,k,x,xx,dx,i:hf,Array,step,eps)=
{
& & hf=HFor(&f&),
& & Array=new[rtoi(real_s),rtoi(45)],
& & step=20,eps=1e-6,& & //step越大，eps越小越精确，用于积分一步的变步长欧拉方法
& & a=new[rtoi(real_s),rtoi(5),rtoi(EndType),1e-50,1e-50,1e-50,1e-50,1e-50],&&//存放常量约束条件中的变量的下界
& & b=new[rtoi(real_s),rtoi(5),rtoi(EndType), 1, 1, 1, 1, 1],&&//存放常量约束条件中的变量的上界
& & x=new[rtoi(real_s),rtoi(6),rtoi(EndType), 0..1..514e-007,9.],&&//其中前n个分量存放初始复形的第一个顶点坐标值（要求满足所有的约束条件），返回时存放极小值点各坐标值。最后一个分量返回极小值。
& & xx=new[rtoi(real_s),rtoi(6),rtoi(10)],
& & _eps=1e-10, alpha=1.01,k=800,dx=1e-4,& &//_eps控制精度要求；反射系数alpha一般取1.3左右；k为允许的最大迭代次数；dx为微小增量。需选择合适的alpha,_eps,k,dx值求解。
& & i=cplx[HFor(&J&),HFor(&S&),a,b,alpha,_eps,x,xx,k,dx],
& & printff{&\r\n实际迭代次数={1,r}\r\n&,i},
& & OutVector[x],
& & delete[Array],delete[a],delete[b],delete[x],delete[xx]
};
验证初值(_k12,_k21,_k23,_k32,_k30:c0,c1,c2,d0,d1,d2,w0,w1,w2:hf,Array,step,eps)=
{
& & hf=HFor(&f&),
& & Array=new[rtoi(real_s),rtoi(45)],
& & step=20,eps=1e-6,& & //step越大，eps越小越精确，用于积分一步的变步长欧拉方法
& & S(_k12,_k21,_k23,_k32,_k30,c0,c1,c2,d0,d1,d2,&w0,&w1,&w2),
& & printff{&\r\n验证总量：d0={1,r}, d1={2,r}, d2={3,r}, \r\n&,w0,w1,w2},
& & printff{&\r\n验证目标函数：J={1,r}\r\n&,J(_k12,_k21,_k23,_k32,_k30)},
& & delete[Array]
};
验证初值(0.1,0.1,0.1,0.1,0.1);&&//第一次选取的初值
验证初值(0..1..514e-007,9.); //第二次选取的初值
复制代码
实际迭代次数=321.
& && &0.112973& &7.& && &0..3& &9.& && & 33.4976
0.
验证总量：d0=99.226, d1=99.794, d2=98.664,
验证目标函数：J=531.58
0.
验证总量：d0=99.795, d1=99.537, d2=99.708,
验证目标函数：J=33.386
0.
复制代码
故最优值可能是：0.112973& &7.& && &0..3& &9.
目标函数值：33.4976
选择合适的alpha,_eps,k,dx值求解时，精度还可能提高。
Forcal下载：
<h1 style="color:#2 麦片财富积分
关注者： 4
以上用的是“求约束条件下n维极值的复形调优法”。
若不考虑约束条件，用“求n维极值的单形调优法”求解为：
//这里是代码窗口，请将Forcal代码写在下面
//若长时间没有结果，按Ctrl+Alt+q退出Forcal运行
i: OutVector(p:k,i)= k=FCDLen(p),printff{&\r\n&},i=0,(i&k).while{printff{&{1,r,14.6}&,get[p,i]},i++},printff
{&\r\n&};& & //输出一维数组
!using[&XSLSF&]; //使用命名空间XSLSF
f(t,X1,X2,X3,dX1,dX2,dX3::k12,k21,k23,k32,k30)={ //函数定义，连分式法对微分方程组积分一步函数pbs1中要用到
&&dX1=-k12*X1+k21*X2,
&&dX2=(-k21-k23)*X2+k12*X1+k32*X3,
&&dX3=(-k32-k30)*X3+k23*X2
};
t_i(hf,a,step,eps,t1,t2,x1,x2,x3:h,i)=&&//用于验证初值
{
& & h=(t2-t1)/step,
& & {& &pbs1[hf,t1,a,h,eps],&&//连分式法对微分方程组积分一步函数pbs1，hf为函数f的句柄
& && &&&t1=t1+h
& & }.until[abs(t1-t2)&h/2],
& & a.getra(0,&x1,&x2,&x3)
};
t_i_2(hf,a,step,eps,t1,t2,x_1,x_2,x_3:x1,x2,x3,h,i)=& & //用于计算目标函数
{
& & h=(t2-t1)/step,
& & {& &pbs1[hf,t1,a,h,eps],&&//连分式法对微分方程组积分一步函数pbs1，hf为函数f的句柄
& && &&&t1=t1+h
& & }.until[abs(t1-t2)&h/2],
& & a.getra(0,&x1,&x2,&x3),
& & (x1-x_1)^2+(x2-x_2)^2+(x3-x_3)^2
};
J(_k12,_k21,_k23,_k32,_k30:t1:hf,Array,step,eps,k12,k21,k23,k32,k30)={& & //目标函数定义
&&k12=_k12,k21=_k21,k23=_k23,k32=_k32,k30=_k30,
&&t1=0, Array.setra(0,100,0,0),
&&t_i_2(hf,Array,step,eps: &t1, 1 : 90,&&8,&&2)+
&&t_i_2(hf,Array,step,eps: &t1, 3 : 73, 19,&&7)+
&&t_i_2(hf,Array,step,eps: &t1, 5 : 60, 14, 23)
};
S(_k12,_k21,_k23,_k32,_k30,c0,c1,c2,d0,d1,d2,w0,w1,w2:t1,x1,x2,x3:hf,Array,step,eps,k12,k21,k23,k32,k30)=& & //约束条件定义
{
&&c0=0,& && && &c1=0,& && && &c2=0,
&&d0=100,& && & d1=100,& && & d2=100,
&&k12=_k12,k21=_k21,k23=_k23,k32=_k32,k30=_k30,
&&t1=0, Array.setra(0,100,0,0),
&&t_i(hf,Array,step,eps: &t1, 1 :&&&x1, &x2, &x3), w0=x1+x2+x3,
&&t_i(hf,Array,step,eps: &t1, 3 :&&&x1, &x2, &x3), w1=x1+x2+x3,
&&t_i(hf,Array,step,eps: &t1, 5 :&&&x1, &x2, &x3), w2=x1+x2+x3
};
main(:d,u,v,x,_eps,k,xx,g,i:hf,Array,step,eps)=
{
& & hf=HFor(&f&),& && && && && && && && && && &&&//模块变量hf保存函数f的句柄，预先用函数HFor获得该句柄
& & Array=new[rtoi(real_s),rtoi(45)],& && && && &//申请工作数组
& & step=20,eps=1e-6,& && && && && && && && && & //积分步数step越大，积分精度eps越小越精确，用于对微分方程组积分一步函数pbs1
& & x=new[rtoi(real_s),rtoi(6)],& && && && && &&&//申请工作数组
& & xx=new[rtoi(real_s),rtoi(5),rtoi(6)],& && &&&//申请工作数组
& & g=new[rtoi(real_s),rtoi(6)],& && && && && &&&//申请工作数组
& & _eps=1e-50, d=1,u=1.6,v=0.4,k=800,& && && &&&//变换d、u、v进一步求解，k为允许的最大迭代次数
& & i=jsim[HFor(&J&),d,u,v,x,_eps,k,xx,g],& && & //求n维极值的单形调优法
& & printff{&\r\n实际迭代次数={1,r}\r\n&,i},& && &//输出实际迭代次数
& & OutVector[x],& && && && && && && && && && &&&//输出最优参数值及目标函数终值
& & delete[x],delete[xx],delete[g],delete[Array] //销毁申请的对象
};
复制代码
实际迭代次数=743.
& && & 0.156e-002& && &0.107414& &&&-0.375394& &&&-0.115403& && & 22.7783
最优参数值：0.156e-002& && &0.107414& &&&-0.375394& &&&-0.115403& && &
目标函数终值：22.7783
<h1 style="color:# 麦片财富积分
非常感谢！！
<h1 style="color:#2 麦片财富积分
关注者： 4
Matlab高手能否给出Matlab解法，比较验证一下。
会1stOpt的朋友也验证一下。
谢谢各位。
<h1 style="color:# 麦片财富积分
forcal好喜欢用forcal语言啊,能介绍下这个语言不?俺还不是很懂.
我在CSDN上也见过你~~~嘿嘿
<h1 style="color:#2 麦片财富积分
关注者： 4
谢谢SWJTU才子的关注！
从这里了解Forcal：http://xoomer.virgilio.it/forcal/sysm/forcal9/forcal.htm
Forcal下载：http://xoomer.virgilio.it/forcal/xiazai/forcal9/forcal32w.rar
<h1 style="color:# 麦片财富积分
如果这个系统需要调整的话，比如添加一个房室， 在你的代码里面有哪些参数需要更改？
<h1 style="color:#2 麦片财富积分
关注者： 4
//这里是代码窗口，请将Forcal代码写在下面
//若长时间没有结果，按Ctrl+Alt+q退出Forcal运行
i: OutVector(p:k,i)= k=FCDLen(p),printff{&\r\n&},i=0,(i&k).while{printff{&{1,r,14.6}&,get[p,i]},i++},printff
{&\r\n&};& & //输出一维数组
!using[&XSLSF&]; //使用命名空间XSLSF
f(t,X1,X2,X3,dX1,dX2,dX3::k12,k21,k23,k32,k30)={ //函数定义，连分式法对微分方程组积分一步函数pbs1中要用到
&&dX1=-k12*X1+k21*X2,& && && && && & //添加房室似乎就是增加微分方程个数，从这里添加。看一下pbs1的说明
&&dX2=(-k21-k23)*X2+k12*X1+k32*X3,& &//如果增加优化参数个数，从k12,k21,k23,k32,k30处添加
&&dX3=(-k32-k30)*X3+k23*X2& && && &&&//这里有了修改，后面相应部分也要修改
};
t_i(hf,a,step,eps,t1,t2,x1,x2,x3:h,i)=&&//用于约束条件定义及验证初值
{
& & h=(t2-t1)/step,
& & {& &pbs1[hf,t1,a,h,eps],&&//连分式法对微分方程组积分一步函数pbs1，hf为函数f的句柄
& && &&&t1=t1+h
& & }.until[abs(t1-t2)&h/2],
& & a.getra(0,&x1,&x2,&x3)
};
t_i_2(hf,a,step,eps,t1,t2,x_1,x_2,x_3:x1,x2,x3,h,i)=& & //用于计算目标函数
{
& & h=(t2-t1)/step,
& & {& &pbs1[hf,t1,a,h,eps],&&//连分式法对微分方程组积分一步函数pbs1，hf为函数f的句柄
& && &&&t1=t1+h
& & }.until[abs(t1-t2)&h/2],
& & a.getra(0,&x1,&x2,&x3),
& & (x1-x_1)^2+(x2-x_2)^2+(x3-x_3)^2
};
J(_k12,_k21,_k23,_k32,_k30:t1:hf,Array,step,eps,k12,k21,k23,k32,k30)={& & //目标函数定义
&&k12=_k12,k21=_k21,k23=_k23,k32=_k32,k30=_k30,
&&t1=0, Array.setra(0,100,0,0),
&&t_i_2(hf,Array,step,eps: &t1, 1 : 90,&&8,&&2)+&&//如果数据量较大，用数组存放数据，然后用循环取数据计算
&&t_i_2(hf,Array,step,eps: &t1, 3 : 73, 19,&&7)+
&&t_i_2(hf,Array,step,eps: &t1, 5 : 60, 14, 23)
};
S(_k12,_k21,_k23,_k32,_k30,c0,c1,c2,d0,d1,d2,w0,w1,w2:t1,x1,x2,x3:hf,Array,step,eps,k12,k21,k23,k32,k30)=& & //约束条件定义
{
&&c0=0,& && && &c1=0,& && && &c2=0,
&&d0=100,& && & d1=100,& && & d2=100,
&&k12=_k12,k21=_k21,k23=_k23,k32=_k32,k30=_k30,
&&t1=0, Array.setra(0,100,0,0),
&&t_i(hf,Array,step,eps: &t1, 1 :&&&x1, &x2, &x3), w0=x1+x2+x3,
&&t_i(hf,Array,step,eps: &t1, 3 :&&&x1, &x2, &x3), w1=x1+x2+x3,
&&t_i(hf,Array,step,eps: &t1, 5 :&&&x1, &x2, &x3), w2=x1+x2+x3
};
main(:a,b,alpha,_eps,k,x,xx,dx,i:hf,Array,step,eps)=
{
& & hf=HFor(&f&),
& & Array=new[rtoi(real_s),rtoi(45)],
& & step=20,eps=1e-6,& & //step越大，eps越小越精确，用于积分一步的变步长欧拉方法
& & a=new[rtoi(real_s),rtoi(5),rtoi(EndType),1e-50,1e-50,1e-50,1e-50,1e-50],&&//存放常量约束条件中的变量的下界
& & b=new[rtoi(real_s),rtoi(5),rtoi(EndType), 1, 1, 1, 1, 1],&&//存放常量约束条件中的变量的上界
& & x=new[rtoi(real_s),rtoi(6),rtoi(EndType), 0..1..514e-007,9.],&&//其中前n个分量存放初始复形的第一个顶点坐标值（要求满足所有的约束条件），返回时存放极小值点各坐标值。最后一个分量返回极小值。
& & xx=new[rtoi(real_s),rtoi(6),rtoi(10)],
& & _eps=1e-10, alpha=1.01,k=800,dx=1e-4,& &//_eps控制精度要求；反射系数alpha一般取1.3左右；k为允许的最大迭代次数；dx为微小增量。需选择合适的alpha,_eps,k,dx值求解。
& & i=cplx[HFor(&J&),HFor(&S&),a,b,alpha,_eps,x,xx,k,dx],&&//cplxd的用法参考其说明
& & printff{&\r\n实际迭代次数={1,r}\r\n&,i},
& & OutVector[x],
& & delete[Array],delete[a],delete[b],delete[x],delete[xx]
};
验证初值(_k12,_k21,_k23,_k32,_k30:c0,c1,c2,d0,d1,d2,w0,w1,w2:hf,Array,step,eps)=
{
& & hf=HFor(&f&),
& & Array=new[rtoi(real_s),rtoi(45)],
& & step=20,eps=1e-6,& & //step越大，eps越小越精确，用于积分一步的变步长欧拉方法
& & S(_k12,_k21,_k23,_k32,_k30,c0,c1,c2,d0,d1,d2,&w0,&w1,&w2),
& & printff{&\r\n验证总量：d0={1,r}, d1={2,r}, d2={3,r}, \r\n&,w0,w1,w2},
& & printff{&\r\n验证目标函数：J={1,r}\r\n&,J(_k12,_k21,_k23,_k32,_k30)},
& & delete[Array]
};
验证初值(0.1,0.1,0.1,0.1,0.1);&&//第一次选取的初值
验证初值(0..1..514e-007,9.); //第二次选取的初值复制代码
重新下载Forcal，看你所用的是否是最新版本。
要熟练使用该程序，（1）看FORCAL新手参考，熟悉Forcal的语法；（2）看XSLSF库中的pbs1和cplx两个函数。
看FORCAL新手参考时如有疑难，及时与我联系。我想知道新手看Forcal的语法有哪些困难，谢谢！
站长推荐 /2
活动地点已更新
Powered by最优化方法报告_百度文库
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
最优化方法报告
上传于||暂无简介
阅读已结束，如果下载本文需要使用1下载券
想免费下载本文？
定制HR最喜欢的简历
下载文档到电脑，查找使用更方便
还剩5页未读，继续阅读
定制HR最喜欢的简历
你可能喜欢最优化建模方法及matlab实现_图文_百度文库
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
评价文档：
最优化建模方法及matlab实现
上传于||暂无简介
大小：890.00KB
登录百度文库，专享文档复制特权，财富值每天免费拿！
你可能喜欢 上传我的文档
 下载
 收藏
该文档贡献者很忙，什么也没留下。
 下载此文档
正在努力加载中...
优化设计实例
下载积分：30
内容提示：机械专业研究生理论课教程
文档格式：PDF|
浏览次数：340|
上传日期： 00:13:52|
文档星级：
该用户还上传了这些文档
优化设计实例
官方公共微信
