PDA

查看完整版本 : [MATLAB混合编程] 在饮酒驾车的离散模型中,如何用最小二乘法求参数k1,k2,v


fyy
2010-05-15, 13:38
1/n:记开始饮酒的时刻为0,并以 1/n小时为步长离散化。
k1:酒精从吸收室向中心室的转移速率。
k2:中心室中的酒精消耗(向外排出、分解或吸收)速率。
p:每次喝下酒精的总量(毫克)。
v:人体液所占的体积(百毫升)。
N: 饮酒次数。
T: 饮酒周期。
k: k=0,1,2……
f(k/n):饮酒后k/n 小时吸收室的酒精量。
g(k/n):饮酒后k/n小时中心室的酒精量。
c(k/n):饮酒后k/n 小时体液的酒精浓度。
由公式上传的缩略图,已知p=11250mg,只要求出系数k1,k2,v,就可确定血液中酒精浓度与实践的函数关系,取n=10,
应该怎么编程得出k1,k2,v
有关参考数据(附件里上传的参考数据)有如下:

anbcjys
2010-05-16, 09:55
有些变量 比如p是多少?

slgu
2010-05-16, 11:55
请运行下列小程序:
clear,clc
x=[0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0];
y=[55.0, 83.0, 97.0, 102.6, 101.9, 95.0, 88.0, 80.7, 73.3, 67.0, 61.0, 56.0, 46.7, 38.9, 32.0, 27.0, 22.5, 18.7, 15.6, 13.0, 10.0, 9.0, 7.5];
b=[0.0168483978389, 0.214620048639, 107.517916808, 0.0000692159359995, 0.112844163331];
fx=@(b,x)b(5)*b(1)*((1-b(2)*x/10).^b(3)-(1-b(1)*x/10).^b(3))./b(4)/(b(1)-b(2));
%[k1,k1,k,v,p]=b
b=lsqcurvefit(fx,b,x,y)
b=nlinfit(x,y,fx,b)
n=length(x)
SSy=var(y')*(n-1)
y1=fx(b,x);
q=(y-y1)*(y-y1)'
rsquare=(SSy-q)/SSy
figure(1),clf
plot(x,y,'o','markersize',10,'markerfacecolor','k')
hold on
x1=0.2:.1:16.5;
y1=fx(b,x1);
plot(x1,y1,'linewidth',2)
legend('data','fit','location','north')
axis tight

fyy
2010-05-16, 12:50
有些变量 比如p是多少?

p=11250mg,n可以取10:)

fyy
2010-05-16, 12:55
请运行下列小程序:
clear,clc
x=[0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0];
y=[55.0, 83.0, 97.0, 102.6, 101.9, 95.0, 88.0, 80.7, 73.3, 67.0, 61.0, 56.0, 46.7, 38.9, 32.0, 27.0, 22.5, 18.7, 15.6, 13.0, 10.0, 9.0, 7.5];
b=[0.0168483978389, 0.214620048639, 107.517916808, 0.0000692159359995, 0.112844163331];
fx=@(b,x)b(5)*b(1)*((1-b(2)*x/10).^b(3)-(1-b(1)*x/10).^b(3))./b(4)/(b(1)-b(2));
%[k1,k1,k,v,p]=b
b=lsqcurvefit(fx,b,x,y)
b=nlinfit(x,y,fx,b)
n=length(x)
SSy=var(y')*(n-1)
y1=fx(b,x);
q=(y-y1)*(y-y1)'
rsquare=(SSy-q)/SSy
figure(1),clf
plot(x,y,'o','markersize',10,'markerfacecolor','k')
hold on
x1=0.2:.1:16.5;
y1=fx(b,x1);
plot(x1,y1,'linewidth',2)
legend('data','fit','location','north')
axis tight
万分感谢您的帮助!再次请教:怎样估计血液中的酒精含量在什么时间最高?因为要求是要离散,所以不可以求导来解.那要怎么办呢?

fyy
2010-05-16, 12:55
有些变量 比如p是多少?
p=11250mg,n可以取10

fyy
2010-05-16, 13:01
请运行下列小程序:
clear,clc
x=[0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0];
y=[55.0, 83.0, 97.0, 102.6, 101.9, 95.0, 88.0, 80.7, 73.3, 67.0, 61.0, 56.0, 46.7, 38.9, 32.0, 27.0, 22.5, 18.7, 15.6, 13.0, 10.0, 9.0, 7.5];
b=[0.0168483978389, 0.214620048639, 107.517916808, 0.0000692159359995, 0.112844163331];
fx=@(b,x)b(5)*b(1)*((1-b(2)*x/10).^b(3)-(1-b(1)*x/10).^b(3))./b(4)/(b(1)-b(2));
%[k1,k1,k,v,p]=b
b=lsqcurvefit(fx,b,x,y)
b=nlinfit(x,y,fx,b)
n=length(x)
SSy=var(y')*(n-1)
y1=fx(b,x);
q=(y-y1)*(y-y1)'
rsquare=(SSy-q)/SSy
figure(1),clf
plot(x,y,'o','markersize',10,'markerfacecolor','k')
hold on
x1=0.2:.1:16.5;
y1=fx(b,x1);
plot(x1,y1,'linewidth',2)
legend('data','fit','location','north')
axis tight
请问这里的k1,k2,v,p各是多少呢?请指教!

slgu
2010-05-16, 15:36
你运行程序了吗?那个最高点一方面可以通过求导得到,另外当然可以通过图形中最高点对应的x坐标得到。

b=[b1,b2,b3,b4,b5]=[k1,k1,k,v,p]

fyy
2010-05-31, 14:45
你运行程序了吗?那个最高点一方面可以通过求导得到,另外当然可以通过图形中最高点对应的x坐标得到。

b=[b1,b2,b3,b4,b5]=[k1,k1,k,v,p]
b=[0.0168483978389, 0.214620048639, 107.517916808, 0.0000692159359995, 0.112844163331];
这里的参数是怎么得来的?是自己赋值的吗?
有没有过程讲解?