微生物生长动力学程序模拟

%% 生长曲线拟合   转载：https://blog.csdn.net/qq_33980829/article/details/113623335
clc
clear
close all
[~, ~, raw] = xlsread('C:\Users\huang\Desktop\生长.xlsx','Sheet1','C2:M16'); %导入文件
%% 创建输出变量
U_1= reshape([raw{:}],size(raw));
%% 清除临时变量
clearvars raw;
U_1(:,2:end)=U_1(:,2:end)/1000;
a_num=1;%标记数值
for i=4:size(U_1,2) %输入要拟合数据内容
    x=U_1(:,1);
    y=U_1(:,i);
    fx=@(b,x)(b(1)./(1+b(2).*exp(-b(3).*x)));
    %x=xlsread('data','a1:b111');输入第一组数据
    % data1
    
    % plot(x,y,'o','markerfacecolor','r')
    
    b=[1 0.6 4.3]; %初始迭代值 最大值 生长速率 （根据具体实验来设定，初始值在本方程拟合影响不大）
    for l=1:30 %拟合过程迭代
        b=lsqcurvefit(fx,b,x,y);
        b=nlinfit(x,y,fx,b);
    end
    n=length(y);
    disp('data1')
    SSy=var(y)*(n-1);
    y1=fx(b,x);
    RSS=(y-y1)'*(y-y1);
    Rsquare(1,a_num)=(SSy-RSS)/SSy;
    x1=linspace(min(x),max(x),300);
    y1=fx(b,x1);
    plot(x,y,'*')
    hold on
    b_num(a_num,:)=b;%保存每次拟合的系数信息
    a=strcat(num2str(b(1)),'./(1+',num2str(b(2)),'.*exp(-',num2str(b(3)),'.*x))'); %拟合的函数
    fplot(a,[min(x),max(x)],'linewidth',1)% 绘制拟合的函数
    %  plot(x1,y1,'-','linewidth',1) %绘制曲线
    %  plot(x1,y1,'b-','linewidth',2.5)
    hold on
    %  plot(x,y,'o')
    %   plot(x,y,'o','markerfacecolor','r')
    %   tex=texlabel(strcat(num2str(b(1)),'./(1+',num2str(b(2)),'.*exp(-',num2str(b(3)),'.*x))'));%公式表示
    %   text(5,1.1,tex) %文本位置
    hold on
    a_num=a_num+1;
    
end
xlabel('生长时间(h)','Interpreter','LaTex')
ylabel('OD$_{600}$','Interpreter','LaTex')
legend('0mg/L','拟合曲线0mg/L','50mg/L','拟合曲线50mg/L','100mg/L','拟合曲线100mg/L','150mg/L','拟合曲线150mg/L',...
    '200mg/L','拟合曲线200mg/L','250mg/L','拟合曲线250mg/L','300mg/L','拟合曲线300mg/L','350mg/L','拟合曲线350mg/L')
%  legend('Position',[0.75928385140126 0.341278249168987 0.0951822931443652 0.312833518330545])
xticks(min(x):2:max(x)+2); %x坐标区间范围以及间隔
hold off
%% 不同生长OD聚类
figure
% y_1=b_num(:,1);
% plot(b_num(:,1),b_num(:,3),'*')

%  ylim([0.3,max(b_num(:,3))])
%  yticks([0.3:0.1:max(b_num)]);
a=[1:8];
 [cidx, ctrs] =kmeans(b_num(:,1:2:3),3);
 plot(b_num(cidx==1,1),b_num(cidx==1,3),'r*', ...
             b_num(cidx==2,1),b_num(cidx==2,3),'b*',b_num(cidx==3,1),b_num(cidx==3,3),'y*');
         a=num2cell(a);
text(b_num(:,1),b_num(:,3)*1.04,a)
xlabel('k')
ylabel('r')

转载：https://blog.csdn.net/qq_33980829/article/details/113623335