生物医学数据分析讲义--医用高等数学实例分析

1、生物医学数据分析讲义-医用高等数学实例分析第3章医用数学根底及其MATLAB实现 The Tower of Babel 沟通的重要性第3章医用数学根底及其MATLAB实现医用高等数学实例分析 函数分析 微分方程定义：设在某个变化过程中存在两个变量x和y，假设对于变量x的每一个允许取的值，变量y按照某一对应关系，都有唯一确定的值与之对应，那么称变量y为变量x的函数，记为其中：x为自变量，y为因变量，f表示二者间的对应函数关系，D为x的定义域，即允许取值的范围。函数分析 函数【例3.1】 静脉注射G钠盐100000单位后，血清中的药物浓度C为时间t的函数C(t)：其中时间t的单位为小时(h)，C的

2、单位为单位每毫升(单位/mL)。试绘制C随时间t变化的曲线。 函数分析函数分析 函数文件的建立% func3f1.mfunction y=func3f1(x)if (x=0)&(x=0.25)&(x=1)&(x3)|(xy=func3f1(2) y=0.2间接调用 函数句柄是用于间接调用一个函数的Matlab值或数据类型。在调用其它函数时可以传递函数句柄，也可在数据结构中保存函数句柄备用。 通过命令形式 fhandle = functionname 可以创立函数句柄 trigFun=sin，或匿名函数sqr = (x) x.2;。 trigFun(1) y=feval(func3f1,2) 函

3、数分析%exam31.m%方法1t=0:0.05:3; for k=1:length(t)C(k)=feval(fun1,t(k);ndfigure,plot(t,C,o-),title(血清中药物浓度随时间变化图);xlabel(t/h);ylabel(C/unit/ml);函数分析%exam31.m%方法2figure,fplot(func3f1,0 3);title(血清中药物浓度随时间变化图)xlabel(t/h);ylabel(C/unit/ml); 函数分析 常用重要的函数分析命令result = func(x0) %求函数值fplot(func,xl xd) %在区间xl xd作

4、图xsolve = fzero(func,x0) %求x0附近的零点Xmin = fminbnd(func,xl,xd)%求区间xl xd内极小值ezplot%函数作图【例3.2】口服一定剂量的药物后，血药浓度C与时间t的关系为其中时间t的单位为小时(h)，C的单位为单位每毫升(单位/ml)。求最高的血药浓度以及到达最高血药浓度的时间，并做出C-t曲线。 函数分析解题思路：判断其单调性：一阶导数图形寻找极大值：可通过对原函数乘以-1，将最大值变为最小值求解。 函数分析%exam32.mtmax=fminbnd(x) -40*(exp(-0.2*x)-exp(-2.3*x),0,24);%找到对

5、应于最大血药浓度的时刻tmaxfmax=feval(x) (40*(exp(-0.2*x)-exp(-2.3*x),tmax); %得到最大血药浓度fplot(x) (40*(exp(-0.2*x)-exp(-2.3*x),0 24,k-)title(血药浓度随时间变化图);%加标题xlabel(t/h);%x坐标轴标注ylabel(C/unit/ml);%y坐标轴标注hold onfplot(x) (-8*exp(-0.2*x)+92*exp(-2.3*x),0 24,k-.)t2test=(x) (-8*exp(-0.2*x)+92*exp(-2.3*x),5);legend(C(t),C

6、(t);plot(tmax,fmax,*,MarkerSize,10);text(tmax-0.1,fmax+4,strcat(t=,num2str(tmax),Cmax=,num2str(fmax);hold off【例3.2】函数分析函数分析 多元函数分析【例3.3】 研究机体对某种药物的反响，设给予药量x单位，经过t小时后机体产生某种反响E，且有：其中a为常量(可允许给予的最大药量)，本例中取a10单位。试画出该函数图形并进行简单的分析。函数分析【例3.3】x=0:0.5:10;t=0:0.5:10;xx,tt=meshgrid(x,t);E=xx.2.*(10-xx).*tt.2.*e

7、xp(-tt);figure,mesh(xx,tt,E,EdgeColor,black);alpha(0.6);%set transparency propertiestitle(药物反响E(x，t);xlabel(给药量/unit),ylabel(时间/hour),zlabel(药物反响);view(55,18);%寻找E的最大值及对应的x，tEmax=max(E(:);i,j=find(E=Emax);xxmax=xx(i,j);ttmax=tt(i,j);text(xxmax,ttmax,Emax+5,strcat(,num2str(xxmax),num2str(ttmax),num2s

8、tr(Emax),);banana = (x) -1*(x(1)2*(10-x(1)*x(2)2*exp(-x(2); x,fval = fminsearch(banana,5, 5)函数分析【例3.3】函数分析 函数拟合没有解析形式的函数关系怎么办？实际问题的处理过程刚好与前面函数作图的过程相反，实验中得到的只有观测数据表，即自变量和因变量的具体值，而函数关系就蕴含于这些观测数据之中。通过观测数据表可以绘出因变量随自变量变化的关系图，这种直观的观察有助于确定应采用哪种形式的函数对数据进行拟合处理或经验估计，通过一定条件下的近似，从而建立解析形式的函数关系。如何根据实验数据来建立变量间函数关系

9、的最正确解析表达式呢？ 函数分析 最小二乘拟合(Least Squares Fitting) 【例3.4】 测得某克山病区10名健康儿童头发与全血中的硒含量，见表3-1，试用最小二乘法建立由发硒x推算血硒y的经验公式。函数分析 最小二乘拟合(Least Squares Fitting) y=kx+b函数分析 最小二乘拟合(Least Squares Fitting) 函数polyfit用于最小二乘拟合如下形式的线性多项式：其用法为 p=polyfit(x，y，n)函数polyval用于计算线性多项式的值其用法为 y=polyval(p,x)函数分析【例3.4】xo=74 66 88 69 91

10、 73 66 96 58 73;yo=13 10 13 11 16 9 7 14 5 10;%将xo按升序排列，便于画连线图。A=sortrows(xo;yo)x=A(:,1);y=A(:,2);%进行不同函数形式的最小二乘拟合p=polyfit(x,y,1); %最小二乘拟合一阶直线方程pp=polyfit(x,y,2);%最小二乘拟合二次多项式ye=polyval(p,x); %计算由拟合直线得到的因变量y值ye2=polyval(pp,x);%计算由拟合直线得到的因变量y值figure,plot(x,y,k*,x,ye,ko-,x,ye2,k-.);legend(实测值,拟合直线,拟合二

11、次多项式);err1=sum(ye-y).2);err2=sum(ye2-y).2);text(70,8,strcat(y=,num2str(p(1),x,num2str(p(2), err=,num2str(err1);text(70,6,strcat(y=,num2str(pp(1),x2+,num2str(pp(2),x,num2str(pp(3),err=,num2str(err2);xlabel(发硒x);ylabel(血硒y);title(x-y散点图及其拟合直线);函数分析 需要特别注意的是，进行拟合时首先要确定待拟合函数的形式，最小二乘法只能保证在拟合函数形式确定的情况下残差平

12、方和最小，并不能保证拟合函数是对原始数据集的最正确拟合函数形式。 关于函数拟合 函数分析幂函数型 关于函数拟合 指数函数型 函数分析 关于函数拟合 曲线拟合工具箱，可拟合各种形式的函数 cftool 微分方程 微分方程数学模型 (1) 寻找改变量，通常是描述方程文字中的变化率 (微商、单位增加量或单位减少量)。 (2) 对问题中的特征进行数学刻画。(3) 用微元法建立微分方程。 (4) 确定微分方程的定解条件(初、边值条件)。 (5) 求解或讨论方程(数值解或定性理论)。 (6) 模型和结果的讨论与分析。微分方程数学模型的建立流程 微分方程 微分方程求解 对于简单的微分方程或微分方程组，可以采

13、用dsolve命令来获得解析解；微分方程 微分方程求解 对于大局部常微分方程，可使用ode45命令来进行数值求解。 其中：ode45表示采用4、5阶的龙格库塔法来求解常微分方程。 tspan表示数值求解的时间范围，如0，10表示0到10秒。 y0表示待求变量的初始值。返回值T为数值求解时间范围内的一系列采样点，Y为对应时间点的待求变量值。 odefun表示待求解的微分方程。需要特别注意的是，odefun总是要写成n个一阶微分方程组的形式。T，Y=ode45(odefun，tspan，y0)微分方程 微分方程求解 对于大局部常微分方程，可使用ode45命令来进行数值求解。 其中：ode45表示采

14、用4、5阶的龙格库塔法来求解常微分方程。 tspan表示数值求解的时间范围，如0，10表示0到10秒。 y0表示待求变量的初始值。返回值T为数值求解时间范围内的一系列采样点，Y为对应时间点的待求变量值。 odefun表示待求解的微分方程。需要特别注意的是，odefun总是要写成n个一阶微分方程组的形式。T，Y=ode45(odefun，tspan，y0)微分方程例如：Van der Pol Equation对于任意阶的常微分方程：利用总可以写为如下的形式：按标准形式改写为：function dydt = DifferentialCoe(t,y)dydt = y(2);(1-y(1)2)*y(2

15、)-y(1);那么odefun可定义为：微分方程【例3.5】细菌繁殖模型：在理想环境中，细菌的繁殖率与细菌的数目成正比。假设 时细菌的数目为 ，求细菌的繁殖规律。function dxdt=func3f3(t,x)dxdt=0.5*xT，Y =ode45(func3f3,0,10,10) 式中的D表示 微分方程【例3.5】当k=0.5时k=0.5;tspan=0,10;x0=10;%解析解画图fplot(t) 10*exp(0.5*t),tspan,k-);%数值解t,x_n=ode45(func3f3,tspan,x0);hold on,plot(t,x_n,ro)legend(解析解,数值

16、解);title(细菌繁殖模型的解);xlabel(时间t)ylabel(细菌数目);hold off微分方程【例3.6】细菌繁殖的非理想环境模型：考虑除系统本身的繁殖外有的细菌向系统外迁移或死亡，损耗速率是时间t的线性函数，即At+B，系统内细菌的繁殖率与细菌的数目成正比，假设 时细菌的数目为 ，求系统内的细菌繁殖规律。解：设t时刻的细菌数目为 ，那么根据题意有 微分方程【例3.6】解析解： 数值解： function dxdt=func3f4(t,x,k,A,B)dxdt=k*x-(A*t+B);应掌握带调用参数的函数的微分方程解法 微分方程【例3.6】k=0.5;tspan=0,10;x

17、0=10;%数值解t1,x_1n=ode45(t,x)func3f4(t,x,0.5,2,1),tspan,x0);%k=0.5,A=2,B=1;t2,x_2n=ode45(t,x)func3f4(t,x,0.5,2.2,1),tspan,x0);%k=0.5,A=2.2,B=1;t3,x_3n=ode45(t,x)func3f4(t,x,0.5,1.8,1),tspan,x0);%k=0.5,A=1.8,B=1;figure,plot(t1,x_1n,r-o,t2,x_2n,b-+,t3,x_3n,k*-);legend(k=0.5,A=2.0,B=1,k=0.5,A=2.2,B=1,k=0

18、.5,A=1.8,B=1);title(细菌繁殖模型的解(非理想环境);xlabel(时间t)ylabel(细菌数目);微分方程【例3.6】图 不同参数对细菌繁殖系统行为的影响 微分方程【例3.7】细菌自然生长模型(Logistic增长模型)：t时刻的细菌数目 由细菌的繁殖率m和死亡率n决定，实践证明：若 时细菌的数目为 ，求系统内的细菌繁殖规律。解：由题意可知，单位时间内细菌的净增加数目应为当前的细菌数目乘以(繁殖率死亡率),即微分方程【例3.7】解析解： 其各阶导数的求法syms r k t x0f=r*x0/(k*x0-(k*x0-r)*exp(-r*t);f1=diff(f,t);%一

19、阶导数f2=diff(f1,t); %二阶导数微分方程【例3.7】tspan=0,10;x0=2;r=1;k=0.2;t=linspace(tspan(1),tspan(2),100);x1=r*x0./(k*x0-(k*x0-r).*exp(-r.*t);r=0.7;x2=r*x0./(k*x0-(k*x0-r).*exp(-r.*t);r=0.4;x3=r*x0./(k*x0-(k*x0-r).*exp(-r.*t);r=0.01;x4=r*x0./(k*x0-(k*x0-r).*exp(-r.*t);jj=find(t5);jj=jj(1);%找到t=5附近的时刻对应的位置figure,

20、plot(t,x1,-,LineWidth,2);text(t(jj),x1(jj)+0.2,r=1,k=0.2);hold on,plot(t,x2,-,LineWidth,2);text(t(jj),x2(jj)+0.2,r=0.5,k=0.2);plot(t,x3,-,LineWidth,2);text(t(jj),x3(jj)+0.2,r=0.4,k=0.2);plot(t,x4,-,LineWidth,2);text(t(jj),x4(jj)+0.2,r=0.01,k=0.2);axis(0,10,0,6)title(不同参数条件下的logistic模型的解);xlabel(时间t)

21、ylabel(细菌数目);hold off;微分方程【例3.7】图 细菌繁殖的Logistic模型解 微分方程【例3.8】Gompertz模型。设v表示t时刻肿瘤的体积(也可以表示肿瘤的大小、肿瘤细胞数目等)，由实际经验得知，某一时刻肿瘤体积v的增长速率与当时的v值成正比，比例系数k随时间t减小，其减小速率与当时k的大小成正比，比例系数 。试确定肿瘤生长规律。 解：根据题意，可建立肿瘤的生长模型：微分方程【例3.8】解析解： 微分方程【例3.8】数值解： %func3f5.mfunction ddt=func3f5(t,x,a)ddt=x(1)*x(2);-a*x(2);tn1,xn=ode45(t,x) func3f5(t,x,a),tspan,v0,k0); 微分方程【例3.8】%绘图fig