生物医学数据分析课件--医用高等数学实例分析

第3章医用数学基础及其MATLAB实现知识结构第3章医用数学基础及其MATLAB实现“Hey,noproblem”TheTowerofBabel沟通的重要性第3章医用数学基础及其MATLAB实现医用高等数学实例分析函数分析微分方程定义：设在某个变化过程中存在两个变量x和y，若对于变量x的每一个允许取的值，变量y按照某一对应关系，都有唯一确定的值与之对应，则称变量y为变量x的函数，记为其中：x为自变量，y为因变量，f表示二者间的对应函数关系，D为x的定义域，即允许取值的范围。函数分析函数()Dyfxx，111(,,)DDnnnyfxxxx，，，【例3.1】静脉注射G钠盐100000单位后，血清中的药物浓度C为时间t的函数C(t)：其中时间t的单位为小时(h)，C的单位为单位每毫升(单位/mL)。试绘制C随时间t变化的曲线。函数分析14.500.25()4.664.260.2510.60.213ttCttttt，≤，≤，≤≤函数分析函数文件的建立%%func3f1.mfunctiony=func3f1(x)if(x=0)&&(x0.25)y=14.5*x;endif(x=0.25)&&(x1)y=4.66-4.26*x;endif(x=1)&&(x=3)y=0.6-0.2*x;endif(x3)||(x0)error('自变量不在该函数定义域内取值')end%%endfunc3f1.m函数分析函数文件的调用直接调用输出参数＝函数名(输入参数)y=func3f1(2)y=0.2间接调用函数句柄是用于间接调用一个函数的Matlab值或数据类型。在调用其它函数时可以传递函数句柄，也可在数据结构中保存函数句柄备用。通过命令形式fhandle=@functionname可以创建函数句柄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');函数分析常用重要的函数分析命令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曲线。函数分析0.22.3()40(ee)ttCt解题思路：判断其单调性：一阶导数图形寻找极大值：可通过对原函数乘以-1，将最大值变为最小值求解。函数分析%%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(-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为常量(可允许给予的最大药量)，本例中取a＝10单位。试画出该函数图形并进行简单的分析。22()etExaxt函数分析【例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(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】函数分析函数拟合没有解析形式的函数关系怎么办？实际问题的处理过程刚好与前面函数作图的过程相反，实验中得到的只有观测数据表，即自变量和因变量的具体值，而函数关系就蕴含于这些观测数据之中。通过观测数据表可以绘出因变量随自变量变化的关系图，这种直观的观察有助于确定应采用哪种形式的函数对数据进行拟合处理或经验估计，通过一定条件下的近似，从而建立解析形式的函数关系。如何根据实验数据来建立变量间函数关系的最佳解析表达式呢？函数分析最小二乘拟合(LeastSquaresFitting)【例3.4】测得某克山病区10名健康儿童头发与全血中的硒含量，见表3-1，试用最小二乘法建立由发硒x推算血硒y的经验公式。函数分析最小二乘拟合(LeastSquaresFitting)y=kx+b2211,min()nniiiiiEykxbE2211,min()nniiiiiEykxbE1120210niiiiniiiEykxbxkEykxbb22iiiiiinxyxyknxx222iiiiiiixyxxybnxx函数分析最小二乘拟合(LeastSquaresFitting)函数polyfit用于最小二乘拟合如下形式的线性多项式：1121()nnnnyxpxpxpxp其用法为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=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('实测值','拟合直线','拟合二次多项式');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)),'x^2+',num2str(pp(2)),'x',num2str(pp(3)),',err=',num2str(err2)));xlabel('发硒x');ylabel('血硒y');title('x-y散点图及其拟合直线');函数分析需要特别注意的是，进行拟合时首先要确定待拟合函数的形式，最小二乘法只能保证在拟合函数形式确定的情况下残差平方和最小，并不能保证拟合函数是对原始数据集的最佳拟合函数形式。关于函数拟合函数分析幂函数型关于函数拟合byaxlnln()lnln()lnlnbbyaxaxabx''yabxebxya指数函数型lnln(e)lnln(e)lnbxbxyaaabx''yabx函数分析关于函数拟合曲线拟合工具箱，可拟合各种形式的函数cftool微分方程微分方程数学模型(1)寻找改变量，通常是描述方程文字中的变化率(微商、单位增加量或单位减少量)。(2)对问题中的特征进行数学刻画。(3)用微元法建立微分方程。(4)确定微分方程的定解条件(初、边值条件)。(5)求解或讨论方程(数值解或定性理论)。(6)模型和结果的讨论与分析。微分方程数学模型的建立流程微分方程微分方程求解对于简单的微分方程或微分方程组，可以采用dsolve命令来获得解析解；d()()0dxtkxtkt，0()ektxtxdsolve('*','(0)0')fDxkxxx微分方程微分方程求解对于大部分常微分方程，可使用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)微分方程微分方程求解对于大部分常微分方程，可使用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)微分方程例如：VanderPolEquation0)1(222ydtdyydtyd1212221)1(yyydtdyydtdy()(1)=(')nnyftyyy-，，，，(1)12=='=nnyyyyyy-，，，122312'='='()nnyyyyyftyyy，，