《{医疗培训课件}生物医学数据分析讲义医用高等数学实例分析》由会员分享，可在线阅读，更多相关《{医疗培训课件}生物医学数据分析讲义医用高等数学实例分析（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】 静脉注射G钠盐100000单位后，血清中的药物浓度C为时间t的函数C(t)： 其中时间t的单位
2、为小时(h)，C的单位为单位每毫升(单位/mL)。试绘制C随时间t变化的曲线。,函数分析,函数分析,函数文件的建立,% func3f1.m function y=func3f1(x) if (x=0) end if (x3)|(x0) error(自变量不在该函数定义域内取值) end % end func3f1.m,函数分析,函数文件的调用,直接调用 输出参数函数名(输入参数) y=func3f1(2) y=0.2 间接调用 函数句柄是用于间接调用一个函数的Matlab值或数据类型。在调用其它函数时可以传递函数句柄，也可在数据结构中保存函数句柄备用。 通过命令形式 fhandle = fun
3、ctionname 可以创建函数句柄 trigFun=sin，或匿名函数sqr = (x) x.2;。 trigFun(1) y=feval(func3f1,2),函数分析,%exam31.m %方法1 t=0:0.05:3; for k=1:length(t) C(k)=feval(fun1,t(k);nd figure,plot(t,C,o-), title(血清中药物浓度随时间变化图); xlabel(t/h); ylabel(C/unit/ml);,函数分析,%exam31.m %方法2 figure,fplot(func3f1,0 3); title(血清中药物浓度随时间变化图) x
4、label(t/h);ylabel(C/unit/ml);,函数分析,常用重要的函数分析命令,result = func(x0) %求函数值 fplot(func,xl xd) %在区间xl xd作图 xsolve = fzero(func,x0) %求x0附近的零点 Xmin = fminbnd(func,xl,xd) %求区间xl xd内极小值 ezplot %函数作图,【例3.2】已知口服一定剂量的药物后，血药浓度C与时间t的关系为 其中时间t的单位为小时(h)，C的单位为单位每毫升(单位/ml)。求最高的血药浓度以及达到最高血药浓度的时间，并做出C-t曲线。,函数分析,解题思路： 判断
5、其单调性：一阶导数图形 寻找极大值：可通过对原函数乘以-1，将最大值变为最小值求解。,函数分析,%exam32.m tmax=fminbnd(x) -40*(exp(-0.2*x)-exp(-2.3*x),0,24); %找到对应于最大血药浓度的时刻tmax fmax=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
6、坐标轴标注 hold on fplot(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(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小时后机体产生
7、某种反应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.*exp(-tt); figure,mesh(xx,tt,E,EdgeColor,black); alpha(0.6);%set transparency properties title(药物反应E(x，t); xlabel(给药量/unit),ylabel(时间/hour),zlabel(药物反应); view(55,18
8、); %寻找E的最大值及对应的x，t Emax=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】,函数分析,函数拟合,没有解析形式的函数关系怎么办？ 实际问题的处理过程刚好与前面函数作
9、图的过程相反，实验中得到的只有观测数据表，即自变量和因变量的具体值，而函数关系就蕴含于这些观测数据之中。通过观测数据表可以绘出因变量随自变量变化的关系图，这种直观的观察有助于确定应采用哪种形式的函数对数据进行拟合处理或经验估计，通过一定条件下的近似，从而建立解析形式的函数关系。 如何根据实验数据来建立变量间函数关系的最佳解析表达式呢？,函数分析,最小二乘拟合(Least Squares Fitting),【例3.4】 测得某克山病区10名健康儿童头发与全血中的硒含量，见表3-1，试用最小二乘法建立由发硒x推算血硒y的经验公式。,函数分析,最小二乘拟合(Least Squares Fitting
10、),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 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,
11、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(7
12、0,6,strcat(y=,num2str(pp(1),x2+,num2str(pp(2),x,num2str(pp(3), err=,num2str(err2); xlabel(发硒x); ylabel(血硒y); title(x-y散点图及其拟合直线);,函数分析,需要特别注意的是，进行拟合时首先要确定待拟合函数的形式，最小二乘法只能保证在拟合函数形式确定的情况下残差平方和最小，并不能保证拟合函数是对原始数据集的最佳拟合函数形式。,关于函数拟合,函数分析,幂函数型,关于函数拟合,指数函数型,函数分析,关于函数拟合,曲线拟合工具箱，可拟合各种形式的函数,cftool,微分方程,微分方程数学模
13、型,(1) 寻找改变量，通常是描述方程文字中的变化率 (微商、单位增加量或单位减少量)。 (2) 对问题中的特征进行数学刻画。 (3) 用微元法建立微分方程。 (4) 确定微分方程的定解条件(初、边值条件)。 (5) 求解或讨论方程(数值解或定性理论)。 (6) 模型和结果的讨论与分析。,微分方程,数学模型的建立流程,微分方程,微分方程求解,对于简单的微分方程或微分方程组，可以采用dsolve命令来获得解析解；,微分方程,微分方程求解,对于大部分常微分方程，可使用ode45命令来进行数值求解。,其中： ode45表示采用4、5阶的龙格库塔法来求解常微分方程。 tspan表示数值求解的时间范围，
14、如0，10表示0到10秒。 y0表示待求变量的初始值。 返回值T为数值求解时间范围内的一系列采样点，Y为对应时间点的待求变量值。 odefun表示待求解的微分方程。需要特别注意的是，odefun总是要写成n个一阶微分方程组的形式。,T，Y=ode45(odefun，tspan，y0),微分方程,微分方程求解,对于大部分常微分方程，可使用ode45命令来进行数值求解。,其中： ode45表示采用4、5阶的龙格库塔法来求解常微分方程。 tspan表示数值求解的时间范围，如0，10表示0到10秒。 y0表示待求变量的初始值。 返回值T为数值求解时间范围内的一系列采样点，Y为对应时间点的待求变量值。
15、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)-y(1);,则odefun可定义为：,微分方程,【例3.5】细菌繁殖模型：在理想环境中，细菌的繁殖率与细菌的数目成正比。若 时细菌的数目为 ，求细菌的繁殖规律。,function dxdt=func3f3(t,x) dxdt=0.5*x T，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(解析解,数值解); title(细
