Аппроксимация данных. Построение модели

4. Потемкин В. Введение  в MATLAB CHM. М.: Диалог-МИФИ, 2000. - 256 с.

5. Пыльнов Ю.В. Регрессионный  анализ полиномиальных моделей. – М.: МИРЭА, 1994, 56 с.

    6. Михальский А.И. Лекционные материалы по курсу КТ в МБС, 2013 год.

 

 

 

 

 

 

 

 

 

 

 

 

 

Приложение

Текст программы в MATLAB:

Y1=[1.35,1.8,2.25,2.7,4.1,4.55,7.3,7.75,8.2,9.1,11.8,15,15.45,18.2];  % Значения объёма по точкам.

T1=[6,7,9,10,11,12,13,14,15,17,18,19,21,22]; % Значения времени по точкам.

V01=1.35;

 

n=14;

P=fminsearch(@G1,1,[],Y1,T1,V01) % Экспоненциальная.

Y2=V01.*exp(T1*P(1));

SKO1 = (1/n)* sum ((Y1 - Y2).^2);

 

P1=fminsearch(@G2,[1,1],[],Y1,T1,V01) % Степенная.

Y3=V01.*(1+T1/(P1(1))).^(P1(2));

SKO2 = (1/n)* sum ((Y1 - Y3).^2);

figure

plot (T1,Y1,'ko-',T1,Y2,'k+:',T1,Y3,'k-.'),grid;

legend('Исходные данные','Экспоненциальная','Степенная');

xlabel('Время');

ylabel('Объём');

 

figure

m1=[];

m2=[];

for V01=0.0001:0.01:3;

    [P3,f]=fminsearch(@G1,1,[],Y1,T1,V01);

    m1=[m1;V01];

    m2=[m2;f];

end

plot(m1,m2,'k-.','LineWidth',2),grid;

xlabel ('Мин. Объём');

ylabel ('СКО');

 

figure

n1=[];

n2=[];

for V01=0.001:0.01:2;

    [P4,f2]=fminsearch(@G2,[1,1],[],Y1,T1,V01);

    n1=[n1;V01];

    n2=[n2;f2];

end

plot(n1,n2,'k-.','LineWidth',2),grid;

xlabel ('Мин. Объём');

ylabel ('СКО');   

 

% До лечения.

V01=0.8200;

a1=0.1440;

Y2=V01*exp(T1.*a1);

 

V02=0.0755;

a2=0.029;

b=2.065;

h=0.1; %Промежуток развития

Y3=V02*diag(1+(T1/a2))^b; % Создадим массив объёма созданного по формуле с учётом оценки.

Y3=diag(Y3);

Y3=Y3';

 

 

figure

plot (T1,Y1,'ko',T1,Y2,'k-.',T1,Y3,'k'); % Сравниваем получившиеся данные, на графике без лечения. В качестве проверки строим экспоненциальную оценку.

grid on

xlabel('Время');

ylabel('Объём');

legend('Исходные данные', 'Экспоненциальная функция','Степенная функция');

 

% Без лечения.

T2m=[0:h:32]; % Время жизни без лечения.

 

Y8=V02*diag(1+(T2m/a2))^b; % Создадим массив объёма созданного по формуле с учётом оценки.

Y8=diag(Y8);

Y8=Y8';

zp3=sum(Y8.*h); % Считаем примерный запас жизненных сил, на примере без лечения.

 

% С лечением.

E=[0:1:5]; % Зависимость эффекта, от максимально переносимой дозы.

MPD=[0:0.2:1];

figure

plot(MPD,E,'k-.');

grid on

xlabel('Переносимая доза');

ylabel('Эффект');

 

 

%Найдем эффект, при дозе n от переносимой.

D=0.8;

tgEMPD=E/MPD; %Находим тангенс при макс. дозе.

MPD=1;

ED=D*MPD*tgEMPD; %Находим эффект при дозе n.

 

t1=7; %Время начала лечения.

n1=3; %Количество инъекций.

td=6; %Промежуток между инъекциями.

tj1=t1+(n1-1)*td; %День последней инъекции.

 

Y41=[];

Y42=[];

Y43=[];

Y44=[];

Tx1=[];

Tx2=[];

Tx3=[];

Tx4=[];

Tj=0; % Введём время развития, максимальное значение которого, на выходе из последнего цикла и будет временем жизни с лечением.

Taux1=t1; % Время 1-ой инъекции.

Taux2=t1+(n1-2)*td; % Время 2-ой инъекции.

Taux3=tj1; % Время последней инъекции.

Tl=31; % Время жизни.

 

 

  while Tj<Taux1

   

    Y4=V02*diag(1+(Tj/a2))^b;

    Y41=[Y41;Y4];

    Tx1=[Tx1;Tj];

Tj=Tj+h;

   

   end;

   Tj1=Tj-ED; % Ввод дозы 1.

   while Tj<Taux2

      Y4=V02*diag(1+(Tj1/a2))^b;

       Y41=[Y41;Y4];

       Tx1=[Tx1;Tj];

       Tj1=Tj1+h;

Tj=Tj+h;

   end;

   Tj2=Tj1-ED; % Ввод дозы 2.

   while Tj<Taux3

   

    Y4=V02*diag(1+(Tj2/a2))^b;

    Y41=[Y41;Y4];

    Tx1=[Tx1;Tj];

    Tj2=Tj2+h;

Tj=Tj+h;

   end;

   Tj3=Tj2-ED; % Ввод дозы 3.

   Zp=sum(Y41.*h);

 

    while Zp<zp3; % Рассчитаем график, до летального объёма.

    Zp=Zp+Y4*h;

    Y4=V02*diag(1+(Tj3/a2))^b;

    Y41=[Y41;Y4];

   

    Tx1=[Tx1;Tj];

    Tj3=Tj3+h;   

Tj=Tj+h;

   end; 

figure

plot(T2m,Y8,'k-.',Tx1,Y41,'k');

grid on

xlabel('Время');

ylabel('Объём');

legend('Без лечения','С лечением');

M-файл G1 для расчета экспоненциальной функции:

function out=G1(P,V1,T1,V01)

a=P(1);

Y=V01.*exp(T1*a);

out=V1-Y;

out=sum(out.*out);

M-файл G2 для расчета степенной функции:

function out=G2(P1,V1,T1,V01)

a=P1(1);

b=P1(2);

Y=V01.*(1+T1/a).^b;

out=V1-Y;

out=sum(out.*out);