湖南人文科技学院-袁野.docVIP

湖南人文科技学院 机电系 机械优化设计 班级：机械三班 姓名：袁野 学号： 教师：禹翼 吸收室c1（t）中心室 吸收室c1（t） 中心室 c（t） k1 k 口服给药方式相当于先有一个将药从肠胃吸收入血液的过程，这个过程可简化为在药物进入中心室之前有一个吸收室（如图），记中心室和吸收室的容积分别为V，V1，而t时刻的血药浓度分别为c（t），c1（t）；中心室的排除速率为k，吸收速率为k1（这里k和k1分别是中心室和吸收室血药浓度变化率与浓度本身的比例系数）。设t=0时刻口服剂量为d的药物，容易写出吸收室的血药浓度c1（t）的微分方程为 中心室血药浓度c（t）的变化率由两部分组成：与c成正比的排除（比例系数ｋ）；与ｃ１成正比的吸收（比例系数ｋ１）。 再考虑到中心室和吸收室的容积分别为Ｖ，Ｖ１，得到ｃ（ｔ）的微分方程为 由以上两个微分方程不难解出中心室血药浓度 在制定给药方案时必须知道这种药物的3个参数k，k1，b(=d/v)，实际中通常通过实验数据确定。设t=0时刻口服一定剂量的药物，下表是实验数据c(t)，请由此确定k，k1，b。 t 0.083 0.167 0.25 0.5 0.75 1 1.5 c(t) 10.9 21.1 27.3 36.4 35.5 38.4 34.8 t 2.25 3 4 6 8 10 12 c(t) 24.2 23.6 15.7 8.2 8.3 2.2 1.8 【模型建立】 问题可以转换为利用数据拟合方程 b=d/V 这个方程是非线性的。由于这是一个非线性最小二乘拟合的问题，故可以选用MATLAB中的Isqnonlin命令或Isqcurvefit命令。 【程序】 function f=medicine(x,t)%x(1)--b,x(2)--k1,x(3)--k f = x(1)*x(2)/(x(2) - x(3))*(exp(-x(3)*t)-exp(-x(2)*t)); 方法一 用Isqcurvefit命令 x0=[1,1,0]; t = [0.083 0.167 0.25 0.50 0.75 1.0 1.5 2.25 3.0 4.0 6.0 8.0 10.0 12.0]; c = [10.9 21.1 27.3 36.4 35.5 38.4 34.8 24.2 23.6 15.7 8.2 8.3 2.2 1.8]; opt = optimset; [x,norm,res,ef,out,lam] = lsqcurvefit(@medicine,x0,t,c,[],[],opt); 【输出结果】 x =46.8275 3.6212 0.2803 norm =34.2317 res = 1.1090 -0.3885 -0.5052 -0.5820 2.2767 -1.4071 -1.6859 2.8020 -1.7058 0.8438 1.2452 -2.9075 0.8787 -0.0423 ef =4 out = iterations: 34 funcCount: 255 stepsize: 1 cgiterations: [] firstorderopt: [] algorithm: medium-scale: Levenberg-Marquardt, line-search message: [1x57 char] lam = upper: [] lower: [] 方法二 用Isqnonlin命令 x0 = [50,1,0]; t = [0.083 0.167 0.25 0.50 0.75 1.0 1.5 2.25 3.0 4.0 6.0 8.0 10.0 12.0]; c = [10.9 21.1 27.3 36.4 35.5 38.4 34.8 24.2 23.6 15.7 8.2 8.3 2.2 1.8]; opt=optimset(LargeScale,on,MaxFunEvals,1000,MaxIter,200); opt=optimset(opt,tolx,1e-16,tolf,1e-16); [x,norm,res,ef,out,lam]=lsqnonlin(@medicine,x0,[],[],[],t,c); 【输出结果】 x =46.8217 3.6219 0.2802 norm =34.2318 res =1.10