Міністерство освіти та науки України
Національний університет “Львівська політехніка”
Контрольна робота 
З курсу
 «Алгоритмізація і програмування»
на тему:
 «Відносна зміна тиску в процесі автоматичного регулювання»
                      
Львів-2008
Мета: побудувати графік зміни тиску в процесі автоматичного регулювання; отримати розв’язок різними методами розв’язування нелінійних рівнянь.
Завдання:  
                      0 c ≤ t ≤ 30 c
Знайти в який момент часу t після зміни тиску відносна зміна вихідної величини (t)=0,6.
Відносна зміна тиску в процесі автоматичного регулювання
Табулювання
Блок-схема алгоритму 
Програма мовою С
//tabyljyvannja
 #include<stdio.h>
 #include<conio.h>
 #include<math.h>
 float fi(float t)
   {return 1/25*t*exp(-(1+0.04)/2/25*t);}
main()
   {clrscr();
    for(float t=0;t<=30;t+=3)
	printf("\n t=%.f fi=%6.4f",t,fi(t));
    getchar();
   }
Результат виконання
 t=0   fi=0.0                    t=12   fi=0.3740               t=24   fi=0.5827
 t=3   fi=0.1127               t=15   fi=0.4392               t=27   fi=0.6159
 t=6   fi=0.2118               t=18   fi=0.4951               t=30   fi=0.6430
 t=9   fi=0.2985               t=21   fi=0.5427
Програма в середовищі MATLAB
%rozraxa2
clear,clc
t=[0:3:30];
fi=roz2(t);
disp([t(:) fi(:)]);
figure(1);
plot(t,fi,'k');
xlabel('t');ylabel('fi(t)');
grid;
Функція-процедура
%roz2
function y=roz2(t)
         T0=25; Kp=0.04; V0=1;
         y=V0/T0.*t.*exp(-(1+Kp)/2./T0.*t);
Результат виконання програми
         0         0               12.0000    0.3740         24.0000    0.5827
    3.0000    0.1127         15.0000    0.4392         27.0000    0.6159 
    6.0000    0.2118         18.0000    0.4951         30.0000    0.6430
    9.0000    0.2985         21.0000    0.5427
Графік залежності густини ртуті від температури
      
Виділення кореня
fi(t)=1/25∙t∙exp(-(1+0.04)/2/25*t)-0,6
t([0;30].                        fi(0)=0;     fi(30)= 0.6430   
Блок-схема алгоритму
Програма мовою С
//vudilennja_korenja
 #include<stdio.h>
 #include<conio.h>
 #include<math.h>
 float fi(float t)
   {return 1/25*t*exp(-(1+0.04)/2/25*t)-0,6;}
main()
   {clrscr();
    for(float t=0;t<=30;t+=3)
	printf("\n t=%.f fi=%6.4f",t,fi(t));
    getchar();
   }
 
Результат виконання
 t=0     f=-0.6000
 t=3     f=-0.4872
 t=6     f=-0.3881
 t=9     f=-0.3014
 t=12   f=-0.2260
 t=15   f=-0.1608
 t=18   f=-0.1048
 t=21   f=-0.0572
 t=24   f=-0.0172
 t=27   f= 0.0159
 t=30   f= 0.0429
Програма в середовищі MATLAB
Функція-процедура
%vud_k
function y=vud_k(t)
         T0=25; Kp=0.04; V0=1;
         y=V0/T0.*t.*exp(-(1+Kp)/2./T0.*t)-0.6;
%vud_kor
clear,clc
t=[0:3:30];
fi=vud_k(t);
disp([t(:) fi(:)]);
format long;
z=fzero('vud_k',0)
figure(1);
plot(t,fi,'k',z,0,'*k');
xlabel('t');ylabel('fi(t)');
grid;
Результат виконання програми
   0             -0.6000          12.0000     -0.2260          24.0000     -0.0172
   3.0000     -0.4872          15.0000     -0.1608          27.0000       0.0159
   6.0000     -0.3881          18.0000     -0.1048          30.0000       0.0429
   9.0000     -0.3014          21.0000     -0.0572
z =25.48761202714193
       
Метод половинного ділення
Блок-схема алгоритму
Програма мовою С
//pol_dil
 #include<stdio.h>
 #include<conio.h>
 #include<math.h>
 float fi(float t)
   {return 1/25*t*exp(-(1+0,04)/2/25*t);}
 main()
   {clrscr();
    float a=0,b=30,t0,eps=1e-4;
    int k=0;
  A:t0=(a+b)/2;
    if(fabs(a-b)>eps)
     {if(fi(a)*fi(t0)>0) a=t0; else b=t0;
      k++; goto A;
     }
    printf("\n t0=%.7f fi(t0)=%.7f",t0,fi(t0));
    printf("  k=%d eps=%.4f a-b=%.7f",k,eps,fabs(a-b));
    getchar();
   }
Результат виконання програми
  t0=0.00002548        fi(t0)=0.0      k=19      eps=0.0001          a-b=0.0000572
Програма в середовищі MATLAB
%pol_dil
clear;clc
a=0; b=30;
k=0; eps=1e-4;
t0=(a+b)/2
while abs(a-b)>eps
    if vud_k(a)*vud_k(t0)>0  a=t0;
       else b=t0;
    end;    
    t0=(a+b)/2;
    k=k+1;
end;
format short;
disp('  t0  ')
disp([  t0  ])
disp('  k   ')
disp([  k   ])
disp('  abs(a-b) ')
disp([abs(a-b)]) 
Результат виконання програми
  t0=25.4876                       k=19               abs(a-b)=5.7220e-005
Метод простих ітерацій
fi(t)=1/25∙t∙exp(-(1+0.04)/2/25*t)-0.6;
fi´(t)=1/25*exp(-(1+0.04)/2/25*t)*(1-(1+0.04)/2/25*t);
φ(t)=t-f(t)∕k  , де |k|≥Q∕2
Q=max|f´(t)|
    [0;30]
fi´(0)=4∙10-2=Q
fi´(30)=0.81∙10-2
fi´(t)>0      =>   k>0
k=2*10-2
φ(t)=t-fi(t)/k; 
q=max|φ´(t)|
    [0;30]
φ´(t)=1- fi´(t)/k;
φ´(0)=0.6777;     φ´(30)=-0.6      =>  q=0.6777<1
Умова виходу:|t1-t|<∆(1-q)∕q;  ε=∆∙(1-q)∕q
Блок-схема алгоритму
 
Програма мовою С
//prostux_iteracij
 #include<stdio.h>
 #include<conio.h>
 #include<math.h>
 float f(float t)
   {return 1/25*t*exp(-(1+0.04)/2/25*t)-0.6;}
 main()
   {clrscr();
    float a=0,b=30,t,t1,det=1e-4,eps=det*(1-0.6777)/0.6777;
    int k=0;
    t=(a+b)/2;
  A:t1= t-(1/25*t*exp(-(1+0.04)/2/25*t)-0.6)/2/1e-2;
    if(fabs(t1-t)>eps)
     {t=t1;   k++; goto A;}
      else
     {printf("\n t1=%.4f f(t1)=%.7f",t1,f(t1));
     printf("  k=%d eps=%.7f ",k,eps);}
     getchar();
   }
Результат виконання програми
 t1=25.4876           f(t1)=0.0000000              k=14            eps=0.0001
Програма в середовищі MATLAB
%pro_itr
clear;clc
a=0; b=30;
k=0; eps=1e-4;
t=(a+b)/2;
t1=t-(1/25*t*exp(-(1+0.04)/2/25*t)-0.6)/2/1e-2;
while abs(t1-t)>eps*(1-0.6777)/0.6777
    t=t1;
    t1=t-(1./25.*t.*exp(-(1+0.04)./2./25*t)-0.6)/2/1e-2;
    k=k+1;
end;
disp('  a   ')  disp([  a   ])
disp('  t1  ')  disp([  t1  ])
disp('  b   ')  disp([  b   ])
disp('  k   ')  disp([  k   ])
disp('  eps ')  disp([  eps ])
Результат виконання програми
  a=0        t1=25.4876         b=30        k=14           eps=1.0000e-004
Метод хорд
fi(a)<0                  fi(b)>0             fi´(t)>0                fi´´(t)>0
t=a-fi(a)∙(b-a)/(fi(b)-fi(a)) ;        γ=fi(t0)∙(t0-tfi)/(fi(t0)-fi(tfi)) ;
t=t0-γ; 
Умова виходу: |γ|≤ε;     ε=∆∙m/(M-m)
m=min |f´(t)|=1.83675∙10-4;       M=max |f´(t)|=1.94287∙10-4;
Програма мовою С
//metod_xord
 #include<stdio.h>
 #include<conio.h>
 #include<math.h>
 float fi(float t)
   {return 1/25∙t∙exp(-(1+0.04)/2/25*t)-0.6;}
 main()
   {clrscr();
    float a=0,b=30,t0,t,tf,ha,det=1e-4,
    m=81*1e-4,M=400*1e-4,eps=m/(M-m)*det;
    int k=0;
    t=a-fi(a)/(fi(b)-fi(a))*(b-a);
    if(fi(a)*fi(t)>0) {t0=a;tf=b;} else {t0=b;tf=a;}
  A:ha=fi(t0)*(t0-tf)/(fi(t0)-fi(tf));
    t=t0-ha;
    if(fabs(ha)>eps)
     {t0=t; k++; goto A;}
    else
    {printf("\n t=%.4f fi(t)=%.7f",t,fi(t));
    printf("  k=%d eps=%.7f fabs(ha)=%.7f",k,eps,fabs(ha));}
    getchar();
   }
Результат виконання програми
 t=25.4876   f(t)=-0.0000    k=8     eps=0.00002539   abs(ha)=0.0000089559
Блок-схема алгоритму
Програма в середовищі MATLAB
%Metod hord
clear;clc
k=0; t0=0; tf=30;
eps=1e-4; M=400*1e-4; m=81*1e-4;
det=eps*m/(M-m);
ha=vud_k(t0)*(t0-tf)/(vud_k(t0)-vud_k(tf));
t=t0-ha;
 while abs(ha)>eps;
    t0=t;
    ha=vud_k(t0)*(t0-tf)/(vud_k(t0)-vud_k(tf));
    t=t0-ha;
    k=k+1; 
 end;
disp(' t '); disp([ t ]);
disp(' k '); disp([ k ]);
format long; 
disp(' eps '); disp([ eps ]);
disp(' abs(ha) '); disp([ abs(ha) ]);
Результат виконання програми
t=25.4876    k=8     eps=2.53918e-005   abs(ha)=8.955989e-006
Метод дотичних
fi´(t)=1/25*exp(-(1+0.04)/2/25*t)*(1-(1+0.04)/2/25*t);
t=b-fi(b)/ fi´(b)
γ=fi(t0)/ fi´(t0)
t=t0-γ;
ε=√2∙m1∙∆/M2
Умова виходу: γ≤√2∙m1∙∆/M2
m1=min|fi´(t)|     =>   m1=fi´(30)=81∙10-4
M2=max|fi´´(t)|
fi´´(t)=-(1+0.04)/2/25/25*exp(-(1+0.04)/2/25*t)*(2-(1+0.04)/2/25*t);
fi´´(0)= -16.64∙10-4 
fi´´(30)= -6.13∙10-4
M2=16.64∙10-4
Програма мовою С
//metog_dotu4nux
 #include<stdio.h>
 #include<conio.h>
 #include<math.h>
 float fi(float t)
   {return 1/25*t*exp(-(1+0.04)/2/25*t)-0.6;}
 main()
   {clrscr();
    float a=0,b=30,t0,t,tf,ha,det=1e-4,m1=81*1e-4,
    M2=16.64*1e-4,eps=sqrt(2*m1*det/M2);
    int k=0;
    t=b-fi(b)/(1/25*exp(-(1+0.04)/2/25*t)*(1-(1+0.04)/2/25*t));
    if(fi(a)*fi(t)>0) {t0=b;tf=a;} else {t0=a;tf=b;}
  A:ha=fi(t0)/(1/25*exp(-(1+0.04)/2/25*t)*(1-(1+0.04)/2/25*t));
    t=t0-ha;
    if(fabs(ha)>eps)
      {t0=t; k++; goto A;}
    else
    {printf("\n t=%.7f fi(t)=%.7f",t,fi(t));
    printf("  k=%d eps=%.7f ha=%.7f ",k,eps,ha);}
    getchar();
}
Результат виконання програми
 t=25.4876            k=2             eps=0.03120005            ha=-0.02104984
Блок-схема алгоритму
Програма в середовищі MATLAB
Функція-процедура
%poxidna
function y=poxidna(t)
         y=1./25.*exp(-(1+0.04)./2./25.*t).*(1-(1+0.04)./2./25.*t);
%metod_dotuchnux
clear;clc
k=0; tf=0; t0=30;
M2=16.64*1e-4; m1=81*1e-4;
det=1e-4;
eps=sqrt(2*m1*det/M2);
ha=vud_k(t0)/poxidna(t0);
t=t0-ha;
 while abs(ha)>eps;
    t0=t;
    k=k+1;
    ha=vud_k(t0)/poxidna(t0);
    t=t0-ha;
 end;
format long;
disp(' t ');   disp([ t ]);
disp(' k ');   disp([ k ]);
disp(' eps '); disp([ eps ]);
disp(' ha ');  disp([ ha ]);
Результат виконання програми
t=25.4876            k=2             eps=0.0312018            ha=-0.02104984
Комбінований метод
Наближення зліва (метод хорд):
γ=fi(a)∙(a-b)/(fi(a)-fi(b)) ;    t1=a-γ.
Наближення справа (метод дотичних):
γ=fi(b)/ fi´(b) ;  t2=b-γ.
Блок-схема алгоритму
Програма мовою С
//kombinovanuj_metod
 #include<stdio.h>
 #include<conio.h>
 #include<math.h>
 float fi(float t)
   {return 1/25*t*exp(-(1+0.04)/2/25*t)-0.6;}
 main()
   {clrscr();
    float a=0,b=30,t2,t1,t,ha,det=1e-4;
    int k=0;
  A:ha=fi(a)*(a-b)/(fi(a)-fi(b));
    t1=a-ha;
    ha=fi(b)/(1/25*exp(-(1+0.04)/2/25*t)*(1-(1+0.04)/2/25*t));
    t2=b-ha;
    if(fabs(b-a)>det)
     {a=t1;b=t2; k++; goto A;}
     else
     {t=(t1+t2)/2;
     printf("\n  t=%.7f",t);
     printf("  k=%d   eps=%.7f ",k,fabs(b-a));}
    getchar();
   }
Результат виконання роботи
t=25.4876         k=3          eps= 0.0000000
Програма в середовищі MATLAB
%Kombinovanui metod
clear;clc
a=0; b=30; eps=1e-4; k=0;
b1=b-vud_k(b)/poxidna(b);
a1=a-(vud_k(a)*(a-b1))/(vud_k(a)-vud_k(b1));
while (abs(b-a)>eps)
    b1=b-vud_k(b)/poxidna(b);
    a1=a-(vud_k(a)*(a-b1))/(vud_k(a)-vud_k(b1));
    b=b1;
    a=a1;
    k=k+1;
end;
t=(a+b)/2;
format long;
disp(' t ');      disp([ t ]);
disp(' k ');      disp([ k ]);
disp('abs(b-a)'); disp([ abs(b-a) ]);
Результат виконання роботи
t=25.4876       k=3
abs(b-a)=1.442278010799214e-005
Порівняльна таблиця результатів обчислень
Метод уточнення кореня
Задана похибка
Значення кореня
Кількість ітерацій
половинного ділення
10-4
   
25.4876
19
простих ітерацій
10-4
   
25.4876
14
хорд
10-4
  
25.4876
8
дотичних
10-4
   
25.4876
2
комбінований
10-4
   
25.4876
3