《生物医学数据分析课件医用高等数学实例分析》由会员分享，可在线阅读，更多相关《生物医学数据分析课件医用高等数学实例分析（45页珍藏版）》请在金锄头文库上搜索。
1、第第3章医用数学基础及其章医用数学基础及其MATLAB实现实现 知识结构知识结构第第3章医用数学基础及其章医用数学基础及其MATLAB实现实现 The Tower of Babel 沟通的重要性第第3章医用数学基础及其章医用数学基础及其MATLAB实现实现医用高等数学实例分析医用高等数学实例分析 函数分析微分方程定义：设在某个变化过程中存在两个变量x和y，若对于变量x的每一个允许取的值，变量y按照某一对应关系，都有唯一确定的值与之对应，则称变量y为变量x的函数，记为其中：x为自变量，y为因变量，f表示二者间的对应函数关系，D为x的定义域，即允许取值的范围。函数分析函数分析 函数函数【例3.1】
2、 静脉注射G钠盐100000单位后，血清中的药物浓度C为时间t的函数C(t)：其中时间t的单位为小时(h)，C的单位为单位每毫升(单位/mL)。试绘制C随时间t变化的曲线。 函数分析函数分析函数分析函数分析 函数文件的建立函数文件的建立%func3f1.mfunctiony=func3f1(x)if(x=0)&(x=0.25)&(x=1)&(x3)|(xy=func3f1(2)y=0.2u间接调用间接调用函数句柄是用于间接调用一个函数的Matlab值或数据类型。在调用其它函数时可以传递函数句柄，也可在数据结构中保存函数句柄备用。通过命令形式fhandle=functionname可以创建函数句
3、柄trigFun=sin，或匿名函数sqr=(x)x.2;。trigFun(1)y=feval(func3f1,2)函数分析函数分析%exam31.m%方法1t=0:0.05:3;fork=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,03);title(血清中药物浓度随时间变化图)xlabel(t/h);ylabel(C/unit/ml);函数分析
4、函数分析 常用重要的函数分析命令常用重要的函数分析命令result=func(x0)%求函数值fplot(func,xlxd)%在区间xlxd作图xsolve=fzero(func,x0)%求x0附近的零点Xmin=fminbnd(func,xl,xd)%求区间xlxd内极小值ezplot%函数作图【例3.2】已知口服一定剂量的药物后，血药浓度C与时间t的关系为其中时间t的单位为小时(h)，C的单位为单位每毫升(单位/ml)。求最高的血药浓度以及达到最高血药浓度的时间，并做出C-t曲线。 函数分析函数分析解题思路：解题思路：判断其单调性：一阶导数图形寻找极大值：可通过对原函数乘以-1，将最大值
5、变为最小值求解。 函数分析函数分析%exam32.mtmax=fminbnd(x)-40*(exp(-0.2*x)-exp(-2.3*x),0,24);%找到对应于最大血药浓度的时刻tmaxfmax=feval(x)(40*(exp(-0.2*x)-exp(-2.3*x),tmax);%得到最大血药浓度fplot(x)(40*(exp(-0.2*x)-exp(-2.3*x),024,k-)title(血药浓度随时间变化图);%加标题xlabel(t/h);%x坐标轴标注ylabel(C/unit/ml);%y坐标轴标注holdonfplot(x)(-8*exp(-0.2*x)+92*exp(-
6、2.3*x),024,k-.)t2test=(x)(-8*exp(-0.2*x)+92*exp(-2.3*x),5);legend(C(t),C(t);plot(tmax,fmax,*,MarkerSize,10);text(tmax-0.1,fmax+4,strcat(t=,num2str(tmax),Cmax=,num2str(fmax);holdoff【例3.2】函数分析函数分析函数分析函数分析 多元函数分析多元函数分析【例3.3】 研究机体对某种药物的反应，设给予药量x单位，经过t小时后机体产生某种反应E，且有：其中a为常量(可允许给予的最大药量)，本例中取a10单位。试画出该函数图形
7、并进行简单的分析。函数分析函数分析【例3.3】x=0:0.5:10;t=0:0.5:10;xx,tt=meshgrid(x,t);E=xx.2.*(10-xx).*tt.2.*exp(-tt);figure,mesh(xx,tt,E,EdgeColor,black);alpha(0.6);%settransparencypropertiestitle(药物反应E(x，t);xlabel(给药量/unit),ylabel(时间/hour),zlabel(药物反应);view(55,18);%寻找E的最大值及对应的x，tEmax=max(E(:);i,j=find(E=Emax);xxmax=xx
8、(i,j);ttmax=tt(i,j);text(xxmax,ttmax,Emax+5,strcat(,num2str(xxmax),num2str(ttmax),num2str(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=
10、kx+b函数分析函数分析 最小二乘拟合最小二乘拟合(Least Squares Fitting) 函数polyfit用于最小二乘拟合如下形式的线性多项式：其用法为 p=polyfit(x，y，n)函数polyval用于计算线性多项式的值其用法为 y=polyval(p,x)函数分析函数分析【例3.4】xo=74668869917366965873;yo=13101311169714510;%将xo按升序排列，便于画连线图。A=sortrows(xo;yo)x=A(:,1);y=A(:,2);%进行不同函数形式的最小二乘拟合p=polyfit(x,y,1);%最小二乘拟合一阶直线方程pp=pol
11、yfit(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(实测值,拟合直线,拟合二次多项式);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+,num2
12、str(pp(2),x,num2str(pp(3),err=,num2str(err2);xlabel(发硒x);ylabel(血硒y);title(x-y散点图及其拟合直线);函数分析函数分析需要特别注意的是，进行拟合时首先要确定待拟合函数的形式，最小二乘法只能保证在拟合函数形式确定的情况下残差平方和最小，并不能保证拟合函数是对原始数据集的最佳拟合函数形式。 关于函数拟合关于函数拟合 函数分析函数分析幂函数型 关于函数拟合关于函数拟合 指数函数型 函数分析函数分析 关于函数拟合关于函数拟合 曲线拟合工具箱，可拟合各种形式的函数 cftool 微分方程微分方程 微分方程数学模型微分方程数学模型
13、 (1) (1) 寻找改变量，通常是描述方程文字中的变化率寻找改变量，通常是描述方程文字中的变化率 ( (微商、单位增加量或单位减少量微商、单位增加量或单位减少量) )。 (2) (2) 对问题中的特征进行数学刻画。对问题中的特征进行数学刻画。(3) (3) 用微元法建立微分方程。用微元法建立微分方程。 (4) (4) 确定微分方程的定解条件确定微分方程的定解条件( (初、边值条件初、边值条件) )。 (5) (5) 求解或讨论方程求解或讨论方程( (数值解或定性理论数值解或定性理论) )。 (6) (6) 模型和结果的讨论与分析模型和结果的讨论与分析。微分方程微分方程数学模型的建立流程数学模
14、型的建立流程 微分方程微分方程 微分方程求解微分方程求解 u对于简单的微分方程或微分方程组，可以采用dsolve命令来获得解析解；微分方程微分方程 微分方程求解微分方程求解 u对于大部分常微分方程，可使用ode45命令来进行数值求解。 其中：(1)ode45表示采用4、5阶的龙格库塔法来求解常微分方程。(2) tspan表示数值求解的时间范围，如0，10表示0到10秒。(3) y0表示待求变量的初始值。(4)返回值T为数值求解时间范围内的一系列采样点，Y为对应时间点的待求变量值。(5) odefun表示待求解的微分方程。需要特别注意的是，odefun总是要写成n个一阶微分方程组的形式。T，Y=
15、ode45(odefun，tspan，y0)微分方程微分方程 微分方程求解微分方程求解 u对于大部分常微分方程，可使用ode45命令来进行数值求解。 其中：(1)ode45表示采用4、5阶的龙格库塔法来求解常微分方程。(2) tspan表示数值求解的时间范围，如0，10表示0到10秒。(3) y0表示待求变量的初始值。(4)返回值T为数值求解时间范围内的一系列采样点，Y为对应时间点的待求变量值。(5) odefun表示待求解的微分方程。需要特别注意的是，odefun总是要写成n个一阶微分方程组的形式。T，Y=ode45(odefun，tspan，y0)微分方程微分方程例如：例如：Van der
16、 Pol Equation对于任意阶的常微分方程：利用总可以写为如下的形式：按标准形式改写为：function dydt = DifferentialCoe(t,y)dydt = y(2);(1-y(1)2)*y(2)-y(1);则odefun可定义为：微分方程微分方程【例3.5】细菌繁殖模型：在理想环境中，细菌的繁殖率与细菌的数目成正比。若 时细菌的数目为 ，求细菌的繁殖规律。functiondxdt=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=
17、10;%解析解画图fplot(t)10*exp(0.5*t),tspan,k-);%数值解t,x_n=ode45(func3f3,tspan,x0);holdon,plot(t,x_n,ro)legend(解析解,数值解);title(细菌繁殖模型的解);xlabel(时间t)ylabel(细菌数目);holdoff微分方程微分方程【例3.6】细菌繁殖的非理想环境模型：考虑除系统本身的繁殖外有的细菌向系统外迁移或死亡，损耗速率是时间t的线性函数，即At+B，系统内细菌的繁殖率与细菌的数目成正比，若 时细菌的数目为 ，求系统内的细菌繁殖规律。解：设t时刻的细菌数目为 ，则根据题意有 微分方程微分
18、方程【例3.6】解析解： 数值解： functiondxdt=func3f4(t,x,k,A,B)dxdt=k*x-(A*t+B);应掌握带调用参数的函数的微分方程解法 微分方程微分方程【例3.6】k=0.5;tspan=0,10;x0=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),t
19、span,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.5,A=1.8,B=1);title(细菌繁殖模型的解(非理想环境);xlabel(时间t)ylabel(细菌数目);微分方程微分方程【例3.6】图不同参数对细菌繁殖系统行为的影响 微分方程微分方程【例3.7】细菌自然生长模型(Logistic增长模型)：t时刻的细菌数目 由细菌的繁殖率m和死亡率n决定，实践证明：若 时细菌的数目为 ，求系统内的细菌繁殖规律。解
20、：由题意可知，单位时间内细菌的净增加数目应为当前的细菌数目乘以(繁殖率死亡率),即微分方程微分方程【例3.7】解析解： 其各阶导数的求法symsrktx0f=r*x0/(k*x0-(k*x0-r)*exp(-r*t);f1=diff(f,t);%一阶导数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=
21、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,plot(t,x1,-,LineWidth,2);text(t(jj),x1(jj)+0.2,r=1,k=0.2);holdon,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
22、=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)ylabel(细菌数目);holdoff;微分方程微分方程【例3.7】图细菌繁殖的Logistic模型解 微分方程微分方程【例3.8】Gompertz模型。设v表示t时刻肿瘤的体积(也可以表示肿瘤的大小、肿瘤细胞数目等)，由实际经验得知，某一时刻肿瘤体积v的增长速率与当时的v值成正比，比例系数k随时间t减小，其减小速率与当时k的大小成正比，比例系
23、数 。试确定肿瘤生长规律。 解：根据题意，可建立肿瘤的生长模型：微分方程微分方程【例3.8】解析解： 微分方程微分方程【例3.8】数值解： %func3f5.mfunctionddt=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】%绘图figure,subplot(2,2,1)plot(tn1,vn1,*,tn1,va1,-);title(v(t),a=0.2);xlabel(时间t)ylabel(肿瘤体积);subplot(2,2,2)plot(t
24、n2,vn2,*,tn2,va2,-);title(v(t),a=1);xlabel(时间t)ylabel(肿瘤体积);subplot(2,2,3)plot(tn1,kn1,*,tn1,ka1,-);title(k(t),a=0.2);xlabel(时间t)ylabel(系数k);subplot(2,2,4)plot(tn2,kn2,*,tn2,ka2,-);legend(数值解,解析解)title(k(t),a=1);xlabel(时间t)ylabel(系数k);tspan=0,50;v0=5;k0=2;%数值解a=0.2;tn1,xn=ode45(t,x)func3f5(t,x,a),ts
25、pan,v0,k0);vn1=xn(:,1);%xn的第1列为v(t)的值,a=0.2kn1=xn(:,2);%xn的第2列为k(t)的值,a=0.2%解析形式对应的值va1=v0*exp(k0/a*(1-exp(-a*tn1);ka1=k0*exp(-a*tn1);%数值解a=1;tn2,xn=ode45(t,x)func3f5(t,x,a),tspan,v0,k0);vn2=xn(:,1);%xn的第1列为v(t)的值,a=1kn2=xn(:,2);%xn的第2列为k(t)的值,a=1va2=v0*exp(k0/a*(1-exp(-a*tn2);ka2=k0*exp(-a*tn2);微分方程微分方程【例3.8】图肿瘤生长Gompertz模型的解
