生物醫(yī)學(xué)數(shù)據(jù)分析講義--醫(yī)用高等數(shù)學(xué)實(shí)例分析

1、生物醫(yī)學(xué)數(shù)據(jù)分析講義-醫(yī)用高等數(shù)學(xué)實(shí)例分析第3章醫(yī)用數(shù)學(xué)根底及其MATLAB實(shí)現(xiàn) The Tower of Babel 溝通的重要性第3章醫(yī)用數(shù)學(xué)根底及其MATLAB實(shí)現(xiàn)醫(yī)用高等數(shù)學(xué)實(shí)例分析 函數(shù)分析 微分方程定義：設(shè)在某個(gè)變化過程中存在兩個(gè)變量x和y，假設(shè)對(duì)于變量x的每一個(gè)允許取的值，變量y按照某一對(duì)應(yīng)關(guān)系，都有唯一確定的值與之對(duì)應(yīng)，那么稱變量y為變量x的函數(shù)，記為其中：x為自變量，y為因變量，f表示二者間的對(duì)應(yīng)函數(shù)關(guān)系，D為x的定義域，即允許取值的范圍。函數(shù)分析 函數(shù)【例3.1】 靜脈注射G鈉鹽100000單位后，血清中的藥物濃度C為時(shí)間t的函數(shù)C(t)：其中時(shí)間t的單位為小時(shí)(h)，C的

2、單位為單位每毫升(單位/mL)。試?yán)L制C隨時(shí)間t變化的曲線。 函數(shù)分析函數(shù)分析 函數(shù)文件的建立% func3f1.mfunction y=func3f1(x)if (x=0)&(x=0.25)&(x=1)&(x3)|(xy=func3f1(2) y=0.2間接調(diào)用 函數(shù)句柄是用于間接調(diào)用一個(gè)函數(shù)的Matlab值或數(shù)據(jù)類型。在調(diào)用其它函數(shù)時(shí)可以傳遞函數(shù)句柄，也可在數(shù)據(jù)結(jié)構(gòu)中保存函數(shù)句柄備用。 通過命令形式 fhandle = functionname 可以創(chuàng)立函數(shù)句柄 trigFun=sin，或匿名函數(shù)sqr = (x) x.2;。 trigFun(1) y=feval(func3f1,2) 函

3、數(shù)分析%exam31.m%方法1t=0:0.05:3; for k=1:length(t)C(k)=feval(fun1,t(k);ndfigure,plot(t,C,o-),title(血清中藥物濃度隨時(shí)間變化圖);xlabel(t/h);ylabel(C/unit/ml);函數(shù)分析%exam31.m%方法2figure,fplot(func3f1,0 3);title(血清中藥物濃度隨時(shí)間變化圖)xlabel(t/h);ylabel(C/unit/ml); 函數(shù)分析 常用重要的函數(shù)分析命令result = func(x0) %求函數(shù)值fplot(func,xl xd) %在區(qū)間xl xd作

4、圖xsolve = fzero(func,x0) %求x0附近的零點(diǎn)Xmin = fminbnd(func,xl,xd)%求區(qū)間xl xd內(nèi)極小值ezplot%函數(shù)作圖【例3.2】口服一定劑量的藥物后，血藥濃度C與時(shí)間t的關(guān)系為其中時(shí)間t的單位為小時(shí)(h)，C的單位為單位每毫升(單位/ml)。求最高的血藥濃度以及到達(dá)最高血藥濃度的時(shí)間，并做出C-t曲線。 函數(shù)分析解題思路：判斷其單調(diào)性：一階導(dǎo)數(shù)圖形尋找極大值：可通過對(duì)原函數(shù)乘以-1，將最大值變?yōu)樽钚≈登蠼狻?函數(shù)分析%exam32.mtmax=fminbnd(x) -40*(exp(-0.2*x)-exp(-2.3*x),0,24);%找到對(duì)

5、應(yīng)于最大血藥濃度的時(shí)刻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(血藥濃度隨時(shí)間變化圖);%加標(biāo)題xlabel(t/h);%x坐標(biāo)軸標(biāo)注ylabel(C/unit/ml);%y坐標(biāo)軸標(biāo)注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】函數(shù)分析函數(shù)分析 多元函數(shù)分析【例3.3】 研究機(jī)體對(duì)某種藥物的反響，設(shè)給予藥量x單位，經(jīng)過t小時(shí)后機(jī)體產(chǎn)生某種反響E，且有：其中a為常量(可允許給予的最大藥量)，本例中取a10單位。試畫出該函數(shù)圖形并進(jìn)行簡(jiǎn)單的分析。函數(shù)分析【例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(時(shí)間/hour),zlabel(藥物反響);view(55,18);%尋找E的最大值及對(duì)應(yīng)的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)函數(shù)分析【例3.3】函數(shù)分析 函數(shù)擬合沒有解析形式的函數(shù)關(guān)系怎么辦？實(shí)際問題的處理過程剛好與前面函數(shù)作圖的過程相反，實(shí)驗(yàn)中得到的只有觀測(cè)數(shù)據(jù)表，即自變量和因變量的具體值，而函數(shù)關(guān)系就蘊(yùn)含于這些觀測(cè)數(shù)據(jù)之中。通過觀測(cè)數(shù)據(jù)表可以繪出因變量隨自變量變化的關(guān)系圖，這種直觀的觀察有助于確定應(yīng)采用哪種形式的函數(shù)對(duì)數(shù)據(jù)進(jìn)行擬合處理或經(jīng)驗(yàn)估計(jì)，通過一定條件下的近似，從而建立解析形式的函數(shù)關(guān)系。如何根據(jù)實(shí)驗(yàn)數(shù)據(jù)來建立變量間函數(shù)關(guān)系

9、的最正確解析表達(dá)式呢？ 函數(shù)分析 最小二乘擬合(Least Squares Fitting) 【例3.4】 測(cè)得某克山病區(qū)10名健康兒童頭發(fā)與全血中的硒含量，見表3-1，試用最小二乘法建立由發(fā)硒x推算血硒y的經(jīng)驗(yàn)公式。函數(shù)分析 最小二乘擬合(Least Squares Fitting) y=kx+b函數(shù)分析 最小二乘擬合(Least Squares Fitting) 函數(shù)polyfit用于最小二乘擬合如下形式的線性多項(xiàng)式：其用法為 p=polyfit(x，y，n)函數(shù)polyval用于計(jì)算線性多項(xiàng)式的值其用法為 y=polyval(p,x)函數(shù)分析【例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);%進(jìn)行不同函數(shù)形式的最小二乘擬合p=polyfit(x,y,1); %最小二乘擬合一階直線方程pp=polyfit(x,y,2);%最小二乘擬合二次多項(xiàng)式y(tǒng)e=polyval(p,x); %計(jì)算由擬合直線得到的因變量y值ye2=polyval(pp,x);%計(jì)算由擬合直線得到的因變量y值figure,plot(x,y,k*,x,ye,ko-,x,ye2,k-.);legend(實(shí)測(cè)值,擬合直線,擬合二

11、次多項(xiàng)式);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(發(fā)硒x);ylabel(血硒y);title(x-y散點(diǎn)圖及其擬合直線);函數(shù)分析 需要特別注意的是，進(jìn)行擬合時(shí)首先要確定待擬合函數(shù)的形式，最小二乘法只能保證在擬合函數(shù)形式確定的情況下殘差平

12、方和最小，并不能保證擬合函數(shù)是對(duì)原始數(shù)據(jù)集的最正確擬合函數(shù)形式。 關(guān)于函數(shù)擬合 函數(shù)分析冪函數(shù)型 關(guān)于函數(shù)擬合 指數(shù)函數(shù)型 函數(shù)分析 關(guān)于函數(shù)擬合 曲線擬合工具箱，可擬合各種形式的函數(shù) cftool 微分方程 微分方程數(shù)學(xué)模型 (1) 尋找改變量，通常是描述方程文字中的變化率 (微商、單位增加量或單位減少量)。 (2) 對(duì)問題中的特征進(jìn)行數(shù)學(xué)刻畫。(3) 用微元法建立微分方程。 (4) 確定微分方程的定解條件(初、邊值條件)。 (5) 求解或討論方程(數(shù)值解或定性理論)。 (6) 模型和結(jié)果的討論與分析。微分方程數(shù)學(xué)模型的建立流程 微分方程 微分方程求解 對(duì)于簡(jiǎn)單的微分方程或微分方程組，可以采

13、用dsolve命令來獲得解析解；微分方程 微分方程求解 對(duì)于大局部常微分方程，可使用ode45命令來進(jìn)行數(shù)值求解。 其中：ode45表示采用4、5階的龍格庫塔法來求解常微分方程。 tspan表示數(shù)值求解的時(shí)間范圍，如0，10表示0到10秒。 y0表示待求變量的初始值。返回值T為數(shù)值求解時(shí)間范圍內(nèi)的一系列采樣點(diǎn)，Y為對(duì)應(yīng)時(shí)間點(diǎn)的待求變量值。 odefun表示待求解的微分方程。需要特別注意的是，odefun總是要寫成n個(gè)一階微分方程組的形式。T，Y=ode45(odefun，tspan，y0)微分方程 微分方程求解 對(duì)于大局部常微分方程，可使用ode45命令來進(jìn)行數(shù)值求解。 其中：ode45表示采

14、用4、5階的龍格庫塔法來求解常微分方程。 tspan表示數(shù)值求解的時(shí)間范圍，如0，10表示0到10秒。 y0表示待求變量的初始值。返回值T為數(shù)值求解時(shí)間范圍內(nèi)的一系列采樣點(diǎn)，Y為對(duì)應(yīng)時(shí)間點(diǎn)的待求變量值。 odefun表示待求解的微分方程。需要特別注意的是，odefun總是要寫成n個(gè)一階微分方程組的形式。T，Y=ode45(odefun，tspan，y0)微分方程例如：Van der Pol Equation對(duì)于任意階的常微分方程：利用總可以寫為如下的形式：按標(biāo)準(zhǔn)形式改寫為：function dydt = DifferentialCoe(t,y)dydt = y(2);(1-y(1)2)*y(2

15、)-y(1);那么odefun可定義為：微分方程【例3.5】細(xì)菌繁殖模型：在理想環(huán)境中，細(xì)菌的繁殖率與細(xì)菌的數(shù)目成正比。假設(shè) 時(shí)細(xì)菌的數(shù)目為 ，求細(xì)菌的繁殖規(guī)律。function dxdt=func3f3(t,x)dxdt=0.5*xT，Y =ode45(func3f3,0,10,10) 式中的D表示 微分方程【例3.5】當(dāng)k=0.5時(shí)k=0.5;tspan=0,10;x0=10;%解析解畫圖fplot(t) 10*exp(0.5*t),tspan,k-);%數(shù)值解t,x_n=ode45(func3f3,tspan,x0);hold on,plot(t,x_n,ro)legend(解析解,數(shù)值

16、解);title(細(xì)菌繁殖模型的解);xlabel(時(shí)間t)ylabel(細(xì)菌數(shù)目);hold off微分方程【例3.6】細(xì)菌繁殖的非理想環(huán)境模型：考慮除系統(tǒng)本身的繁殖外有的細(xì)菌向系統(tǒng)外遷移或死亡，損耗速率是時(shí)間t的線性函數(shù)，即At+B，系統(tǒng)內(nèi)細(xì)菌的繁殖率與細(xì)菌的數(shù)目成正比，假設(shè) 時(shí)細(xì)菌的數(shù)目為 ，求系統(tǒng)內(nèi)的細(xì)菌繁殖規(guī)律。解：設(shè)t時(shí)刻的細(xì)菌數(shù)目為 ，那么根據(jù)題意有 微分方程【例3.6】解析解： 數(shù)值解： function dxdt=func3f4(t,x,k,A,B)dxdt=k*x-(A*t+B);應(yīng)掌握帶調(diào)用參數(shù)的函數(shù)的微分方程解法 微分方程【例3.6】k=0.5;tspan=0,10;x

17、0=10;%數(shù)值解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(細(xì)菌繁殖模型的解(非理想環(huán)境);xlabel(時(shí)間t)ylabel(細(xì)菌數(shù)目);微分方程【例3.6】圖 不同參數(shù)對(duì)細(xì)菌繁殖系統(tǒng)行為的影響 微分方程【例3.7】細(xì)菌自然生長(zhǎng)模型(Logistic增長(zhǎng)模型)：t時(shí)刻的細(xì)菌數(shù)目 由細(xì)菌的繁殖率m和死亡率n決定，實(shí)踐證明：若 時(shí)細(xì)菌的數(shù)目為 ，求系統(tǒng)內(nèi)的細(xì)菌繁殖規(guī)律。解：由題意可知，單位時(shí)間內(nèi)細(xì)菌的凈增加數(shù)目應(yīng)為當(dāng)前的細(xì)菌數(shù)目乘以(繁殖率死亡率),即微分方程【例3.7】解析解： 其各階導(dǎo)數(shù)的求法syms r k t x0f=r*x0/(k*x0-(k*x0-r)*exp(-r*t);f1=diff(f,t);%一

19、階導(dǎo)數(shù)f2=diff(f1,t); %二階導(dǎo)數(shù)微分方程【例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附近的時(shí)刻對(duì)應(yīng)的位置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(不同參數(shù)條件下的logistic模型的解);xlabel(時(shí)間t)

21、ylabel(細(xì)菌數(shù)目);hold off;微分方程【例3.7】圖 細(xì)菌繁殖的Logistic模型解 微分方程【例3.8】Gompertz模型。設(shè)v表示t時(shí)刻腫瘤的體積(也可以表示腫瘤的大小、腫瘤細(xì)胞數(shù)目等)，由實(shí)際經(jīng)驗(yàn)得知，某一時(shí)刻腫瘤體積v的增長(zhǎng)速率與當(dāng)時(shí)的v值成正比，比例系數(shù)k隨時(shí)間t減小，其減小速率與當(dāng)時(shí)k的大小成正比，比例系數(shù) 。試確定腫瘤生長(zhǎng)規(guī)律。 解：根據(jù)題意，可建立腫瘤的生長(zhǎng)模型：微分方程【例3.8】解析解： 微分方程【例3.8】數(shù)值解： %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