% Математичні та програмні засоби моделювання ІВС % Лабораторна робота № 1 % Дослідження чисельних методів моделювання ІВС на ЕОМ % Дифереціальне рівняння першого порядку dy(t)/dt+ay(t)=bx(t) % - - - Стандартні методи рішення ode45 ode23 ode113, % всі методи, що програмуються, та оцінка точності методів - - - % --- Початкові дані --- clear variables; % Диференціальне рівняння dy(t)/dt+ay(t)=bx(t) % Передаточна функція W(p)=K/(Tp+1) % a=1/T b=K/T T=1/a K=b/a a=1; b=1; % Початкова умова y(0) y0=0; % Вхідний сигнал системи x(t) % тип сигналу: step - ступеневий вплив; % sin - синусоїда; % sin_exp - синусоїда, що згасає Type_input_x='step'; % амплітуда сигналу, В A_input_x=1; % частота сигналу, Гц f_input_x=0.1; % Інтервал спостереження, с Tmax=10; % Крок чисельного методу, с tau=0.01; % Автоматичний вибір кроку чисельного методу для ode45 ode23 ode113, так=1/ні=0 tau_auto=0; % --- Моделювання --- % Формування вхідного сигналу Input_x, потрібен тільки для відображення на % графіку разом з рішенням диференціального рівняння та для методів, що програмуються по крокам % В функції ODE_IMS вхідний сигнал обчислюється окремо з врахуванням поточного кроку інтегрування ts=(0:tau:Tmax)'; Input_x=zeros(size(ts)); % Функція ODE_IMS для обчислення правої частини диференціального рівняння switch lower(Type_input_x) case 'step' Input_x(:)=A_input_x; ODE_IMS=@(t, y)[-a*y+b*A_input_x]; case 'sin' Input_x=A_input_x*sin(2*pi*f_input_x*ts); ODE_IMS=@(t, y)[-a*y+b*A_input_x*sin(2*pi*f_input_x*t)]; case 'sin_exp' Input_x=A_input_x*exp(-ts).*sin(2*pi*f_input_x*ts); ODE_IMS=@(t, y)[-a*y+b*A_input_x*exp(-t).*sin(2*pi*f_input_x*t)]; end; % Визначення відліків часу на інтервалі інтегрування Tspan для чисельного методу if tau_auto==1 Tspan=[0 Tmax]; else Tspan=0:tau:Tmax; end; % Задання параметрів методу чисельного рішення диференціального рівняння % відносна похибка рішення, абсолютна похибка рішення, максимальний розмір кроку Option_Lab_1 = odeset('RelTol', 0.001, 'AbsTol', 1e-6, 'MaxStep', Tmax/10); % Рішення диференціального рівняння чисельним методом % Yout - вектор-стовбець рішень диференціального рівняння першого порядку, % в якому кожен елемент відповідає моменту часу у вектор-стовбці Tout % ode45 - однокроковий явний метод Рунге-Кутта 4 та 5 порядків tic; [Tout45,Yout45] = ode45(ODE_IMS, Tspan, y0, Option_Lab_1); Time45=toc; % ode23 - однокроковий явний метод Рунге-Кутта 2 та 3 порядків tic; [Tout23,Yout23] = ode23(ODE_IMS, Tspan, y0, Option_Lab_1); Time23=toc; % ode113 - багатокроковий метод Адамса-Башфорта_Мултона tic; [Tout113,Yout113] = ode113(ODE_IMS, Tspan, y0, Option_Lab_1); Time113=toc; % метод прямокутників (Ейлера) YoutRect=zeros(size(ts)); YoutRect(1)=y0; tic; for i=2:size(ts,1) YoutRect(i)=(1-tau*a)*YoutRect(i-1)+tau*b*Input_x(i-1); end; TimeRect=toc; % метод трапецій YoutTrap=zeros(size(ts)); YoutTrap(1)=y0; tic; for i=2:size(ts,1) YoutTrap(i)=((1-tau*a/2)*YoutTrap(i-1)+tau*b*Input_x(i)/2+tau*b*Input_x(i-1)/2)/(1+tau*a/2); end; TimeTrap=toc; % метод Симпсона YoutSimp=zeros(size(ts)); YoutSimp=YoutTrap; % оцінка значень функції під інтегралом за методом трапецій tic; for i=2:(size(ts,1)-1) YoutSimp(i+1)=YoutSimp(i-1)+tau/3*(-a*YoutSimp(i-1)+b*Input_x(i-1))+... 4*tau/3*(-a*YoutSimp(i)+b*Input_x(i))+tau/3*(-a*YoutSimp(i+1)+b*Input_x(i+1)); end; TimeSimp=toc; % метод модифікований Ейлера (передбачення методом прямокутників та виправлення методом трапецій) YoutReTr=zeros(size(ts)); YoutReTr(1)=y0; tic; for i=2:size(ts,1) p=YoutReTr(i-1)+tau*(-a*YoutReTr(i-1)+b*Input_x(i-1)); YoutReTr(i)=YoutReTr(i-1)+(tau/2)*(-a*YoutReTr(i-1)+b*Input_x(i-1)-a*p+b*Input_x(i)); end; TimeReTr=toc; % метод Мілна 4-го порядку багатокроковий YoutMiln=zeros(size(ts)); YoutMiln(1)=y0; % перші значення за методом прямокутників YoutMiln(2)=YoutRect(2); YoutMiln(3)=YoutRect(3); YoutMiln(4)=YoutRect(4); tic; for i=5:size(ts,1) p=YoutMiln(i-4)+(tau/3)*(2*(-a*YoutMiln(i-1)+b*Input_x(i-1))-(-a*YoutMiln(i-2)+b*Input_x(i-2))+... 2*(-a*YoutMiln(i-3)+b*Input_x(i-3))); YoutMiln(i)=YoutMiln(i-2)+(tau/3)*((-a*p+b*Input_x(i))+4*(-a*YoutMiln(i-1)+b*Input_x(i-1))+... (-a*YoutMiln(i-2)+b*Input_x(i-2))); end; TimeMiln=toc; % метод Адамса-Мултона багатокроковий YoutAMult=zeros(size(ts)); YoutAMult(1)=y0; % перші значення за методом прямокутників YoutAMult(2)=YoutRect(2); YoutAMult(3)=YoutRect(3); YoutAMult(4)=YoutRect(4); tic; for i=5:size(ts,1) p=YoutAMult(i-1)+(tau/24)*(55*(-a*YoutAMult(i-1)+b*Input_x(i-1))-59*(-a*YoutAMult(i-2)+b*Input_x(i-2))+... 37*(-a*YoutAMult(i-3)+b*Input_x(i-3))-9*(-a*YoutAMult(i-4)+b*Input_x(i-4))); YoutAMult(i)=YoutAMult(i-1)+(tau/24)*(9*(-a*p+b*Input_x(i))+19*(-a*YoutAMult(i-1)+b*Input_x(i-1))-... 5*(-a*YoutAMult(i-2)+b*Input_x(i-2))+(-a*YoutAMult(i-3)+b*Input_x(i-3))); end; TimeAMult=toc; % метод Рунге-Кутта 4 порядку YoutRC4=zeros(size(ts)); YoutRC4(1)=y0; tic; for i=2:size(ts,1) k1=-a*YoutRC4(i-1)+b*Input_x(i-1); k2=-a*(YoutRC4(i-1)+(tau/2)*k1)+b*Input_x(i-1); k3=-a*(YoutRC4(i-1)+(tau/2)*k2)+b*Input_x(i-1); k4=-a*(YoutRC4(i-1)+(tau/2)*k3)+b*Input_x(i-1); YoutRC4(i)=YoutRC4(i-1)+(tau/6)*(k1+2*k2+2*k3+k4); end; TimeRC4=toc; % Обчислення значень точного рішення діференціального рівняння switch lower(Type_input_x) case 'step' YoutPrecise=b/a+(y0-b/a)*exp(-a*ts); YoutPrecise45=b/a+(y0-b/a)*exp(-a*Tout45); YoutPrecise23=b/a+(y0-b/a)*exp(-a*Tout23); YoutPrecise113=b/a+(y0-b/a)*exp(-a*Tout113); case 'sin' YoutPrecise=b*A_input_x*(a*sin(2*pi*f_input_x*ts)-2*pi*f_input_x*cos(2*pi*f_input_x*ts))/(a^2+(2*pi*f_input_x)^2)+... (y0+b*A_input_x*2*pi*f_input_x/(a^2+(2*pi*f_input_x)^2))*exp(-a*ts); YoutPrecise45=b*A_input_x*(a*sin(2*pi*f_input_x*Tout45)-2*pi*f_input_x*cos(2*pi*f_input_x*Tout45))/(a^2+(2*pi*f_input_x)^2)+... (y0+b*A_input_x*2*pi*f_input_x/(a^2+(2*pi*f_input_x)^2))*exp(-a*Tout45); YoutPrecise23=b*A_input_x*(a*sin(2*pi*f_input_x*Tout23)-2*pi*f_input_x*cos(2*pi*f_input_x*Tout23))/(a^2+(2*pi*f_input_x)^2)+... (y0+b*A_input_x*2*pi*f_input_x/(a^2+(2*pi*f_input_x)^2))*exp(-a*Tout23); YoutPrecise113=b*A_input_x*(a*sin(2*pi*f_input_x*Tout113)-2*pi*f_input_x*cos(2*pi*f_input_x*Tout113))/(a^2+(2*pi*f_input_x)^2)+... (y0+b*A_input_x*2*pi*f_input_x/(a^2+(2*pi*f_input_x)^2))*exp(-a*Tout113); case 'sin_exp' YoutPrecise=b*A_input_x*exp(-ts).*((a-1)*sin(2*pi*f_input_x*ts)-2*pi*f_input_x*cos(2*pi*f_input_x*ts))/((a-1)^2+(2*pi*f_input_x)^2)+... (y0+b*A_input_x*2*pi*f_input_x/((a-1)^2+(2*pi*f_input_x)^2))*exp(-a*ts); YoutPrecise45=b*A_input_x*exp(-Tout45).*((a-1)*sin(2*pi*f_input_x*Tout45)-2*pi*f_input_x*cos(2*pi*f_input_x*Tout45))/((a-1)^2+(2*pi*f_input_x)^2)+... (y0+b*A_input_x*2*pi*f_input_x/((a-1)^2+(2*pi*f_input_x)^2))*exp(-a*Tout45); YoutPrecise23=b*A_input_x*exp(-Tout23).*((a-1)*sin(2*pi*f_input_x*Tout23)-2*pi*f_input_x*cos(2*pi*f_input_x*Tout23))/((a-1)^2+(2*pi*f_input_x)^2)+... (y0+b*A_input_x*2*pi*f_input_x/((a-1)^2+(2*pi*f_input_x)^2))*exp(-a*Tout23); YoutPrecise113=b*A_input_x*exp(-Tout113).*((a-1)*sin(2*pi*f_input_x*Tout113)-2*pi*f_input_x*cos(2*pi*f_input_x*Tout113))/((a-1)^2+(2*pi*f_input_x)^2)+... (y0+b*A_input_x*2*pi*f_input_x/((a-1)^2+(2*pi*f_input_x)^2))*exp(-a*Tout113); end; % Оцінка точності методів чисельного рішення диференціального рівняння Delta45=zeros(size(Yout45)); Delta23=zeros(size(Yout23)); Delta113=zeros(size(Yout113)); DeltaRect=zeros(size(YoutRect)); DeltaTrap=zeros(size(YoutTrap)); DeltaSimp=zeros(size(YoutSimp)); DeltaReTr=zeros(size(YoutReTr)); DeltaMiln=zeros(size(YoutMiln)); DeltaAMult=zeros(size(YoutAMult)); DeltaRC4=zeros(size(YoutRC4)); Delta45=Yout45-YoutPrecise45; Delta23=Yout23-YoutPrecise23; Delta113=Yout113-YoutPrecise113; DeltaRect=YoutRect-YoutPrecise; DeltaTrap=YoutTrap-YoutPrecise; DeltaSimp=YoutSimp-YoutPrecise; DeltaReTr=YoutReTr-YoutPrecise; DeltaMiln=YoutMiln-YoutPrecise; DeltaAMult=YoutAMult-YoutPrecise; DeltaRC4=YoutRC4-YoutPrecise; Delta45Max=max(abs(Delta45)); Delta23Max=max(abs(Delta23)); Delta113Max=max(abs(Delta113)); Delta45Mean=mean(Delta45); Delta23Mean=mean(Delta23); Delta113Mean=mean(Delta113); Delta45Std=std(Delta45); Delta23Std=std(Delta23); Delta113Std=std(Delta113); DeltaRectMax=max(abs(DeltaRect)); DeltaRectMean=mean(DeltaRect); DeltaRectStd=std(DeltaRect); DeltaTrapMax=max(abs(DeltaTrap)); DeltaTrapMean=mean(DeltaTrap); DeltaTrapStd=std(DeltaTrap); DeltaSimpMax=max(abs(DeltaSimp)); DeltaSimpMean=mean(DeltaSimp); DeltaSimpStd=std(DeltaSimp); DeltaReTrMax=max(abs(DeltaReTr)); DeltaReTrMean=mean(DeltaReTr); DeltaReTrStd=std(DeltaReTr); DeltaMilnMax=max(abs(DeltaMiln)); DeltaMilnMean=mean(DeltaMiln); DeltaMilnStd=std(DeltaMiln); DeltaAMultMax=max(abs(DeltaAMult)); DeltaAMultMean=mean(DeltaAMult); DeltaAMultStd=std(DeltaAMult); DeltaRC4Max=max(abs(DeltaRC4)); DeltaRC4Mean=mean(DeltaRC4); DeltaRC4Std=std(DeltaRC4); % Оцінка швидкодії методів чисельного рішення диференціального рівняння % тільки для ступеневого вхідного сигналу, фіксується час досягнення % рівня 95% від сталого значення амплітуди вихідного сигналу if strcmp(Type_input_x,'step')==1 NState=find(Yout45>=(0.95*A_input_x)); if size(NState,1)==0, NState(1)=size(Yout45,1); end; TState45=Tout45(NState(1)); NState=find(Yout23>=(0.95*A_input_x)); if size(NState,1)==0, NState(1)=size(Yout23,1); end; TState23=Tout23(NState(1)); NState=find(Yout113>=(0.95*A_input_x)); if size(NState,1)==0, NState(1)=size(Yout113,1); end; TState113=Tout113(NState(1)); NState=find(YoutRect>=(0.95*A_input_x)); if size(NState,1)==0, NState(1)=size(YoutRect,1); end; TStateRect=ts(NState(1)); NState=find(YoutTrap>=(0.95*A_input_x)); if size(NState,1)==0, NState(1)=size(YoutTrap,1); end; TStateTrap=ts(NState(1)); NState=find(YoutSimp>=(0.95*A_input_x)); if size(NState,1)==0, NState(1)=size(YoutSimp,1); end; TStateSimp=ts(NState(1)); NState=find(YoutReTr>=(0.95*A_input_x)); if size(NState,1)==0, NState(1)=size(YoutReTr,1); end; TStateReTr=ts(NState(1)); NState=find(YoutMiln>=(0.95*A_input_x)); if size(NState,1)==0, NState(1)=size(YoutMiln,1); end; TStateMiln=ts(NState(1)); NState=find(YoutAMult>=(0.95*A_input_x)); if size(NState,1)==0, NState(1)=size(YoutAMult,1); end; TStateAMult=ts(NState(1)); NState=find(YoutRC4>=(0.95*A_input_x)); if size(NState,1)==0, NState(1)=size(YoutRC4,1); end; TStateRC4=ts(NState(1)); else TState45=0; TState23=0; TState113=0; TStateRect=0; TStateTrap=0; TStateSimp=0; TStateReTr=0; TStateMiln=0; TStateAMult=0; TStateRC4=0; end; % --- Виведення результатів моделювання --- figure; plot(ts, Input_x, 'k--', ts, YoutPrecise, 'k-', ... Tout45, Yout45(:,1), 'r-', Tout23, Yout23(:,1), 'g-.', Tout113, Yout113(:,1), 'b:'); legend('вхід', 'вихід точно', 'вихід ode45', 'вихід ode23', 'вихід ode113'); title('Моделювання систем першого порядку чис. мет. MATLAB'); xlabel('Час, с'); ylabel('Амплітуда сигналів входу і виходу, В'); grid on; figure; plot(ts, Input_x, 'k--', ts, YoutPrecise, 'k-', ... ts, YoutRect, 'r-', ts, YoutTrap, 'g-', ts, YoutSimp, 'b-', ... ts, YoutReTr, 'r--', ts, YoutMiln, 'g--', ts, YoutAMult, 'b--', ts, YoutRC4, 'r-.'); legend('вхід', 'вихід точно', 'вихід метод прямок.', 'вихід метод трап.', 'вихід метод Симпс.', ... 'вихід метод мод.Ейл.', 'вихід метод Мілна', 'вихід метод Адамс.М', 'вихід метод Р-К 4п'); title('Моделювання систем першого порядку чис. мет., що програмуються'); xlabel('Час, с'); ylabel('Амплітуда сигналів входу і виходу, В'); grid on; fprintf(1,'Лабораторна робота № 1\n'); fprintf(1,'Дослідження чисельних методів моделювання ІВС на ЕОМ\n'); fprintf(1,'Чисельні методи рішення диф рівняння ode45, ode23, ode113 та 7 методів, що програмуються по крокам\n'); fprintf(1,'Диференціальне рівняння dy/dt+ay(t)=bx(t), a = %7.3f, b = %7.3f\n', a, b); fprintf(1,'Початкова умова y(0) = %7.3f\n', y0); fprintf(1,'Вхідний сигнал системи x(t) %s\n', Type_input_x); fprintf(1,'Амплітуда сигналу, В та частота сигналу, Гц %7.3f %7.3f\n', A_input_x, f_input_x); fprintf(1,'Автоматичний вибір кроку чисельного методу для ode45 ode23 ode113, так=1/ні=0 %7.3f\n', tau_auto); fprintf(1,'Інтервал спостереження, с та Крок чисельного методу, с %7.3f %7.3f\n', Tmax, tau); fprintf(1,'Час вирішення диф рівняння для ode45 ode23 ode113, с %10.5f %10.5f %10.5f\n',... Time45, Time23, Time113); fprintf(1,'Час вирішення диф рівняння для методів, що програмуються по крокам , с %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f\n',... TimeRect, TimeTrap, TimeSimp, TimeReTr, TimeMiln, TimeAMult, TimeRC4); fprintf(1,'Оцінки амплітудних похибок (В) та швидкодії (с) вихідного сигналу для чисельних методів\n'); fprintf(1,' макс зн середне зн СКЗ, В швидкодія, с\n'); fprintf(1,'Метод ode45 %10.5f %10.5f %10.5f %10.5f\n', Delta45Max, Delta45Mean, Delta45Std, TState45); fprintf(1,'Метод ode23 %10.5f %10.5f %10.5f %10.5f\n', Delta23Max, Delta23Mean, Delta23Std, TState23); fprintf(1,'Метод ode113 %10.5f %10.5f %10.5f %10.5f\n', Delta113Max, Delta113Mean, Delta113Std, TState113); fprintf(1,'Метод прямокутників %10.5f %10.5f %10.5f %10.5f\n', DeltaRectMax, DeltaRectMean, DeltaRectStd, TStateRect); fprintf(1,'Метод трапецій %10.5f %10.5f %10.5f %10.5f\n', DeltaTrapMax, DeltaTrapMean, DeltaTrapStd, TStateTrap); fprintf(1,'Метод Симпсона %10.5f %10.5f %10.5f %10.5f\n', DeltaSimpMax, DeltaSimpMean, DeltaSimpStd, TStateSimp); fprintf(1,'Метод модиф.Ейлера %10.5f %10.5f %10.5f %10.5f\n', DeltaReTrMax, DeltaReTrMean, DeltaReTrStd, TStateReTr); fprintf(1,'Метод Мілна 4п %10.5f %10.5f %10.5f %10.5f\n', DeltaMilnMax, DeltaMilnMean, DeltaMilnStd, TStateMiln); fprintf(1,'Метод Адамса-Мултона%10.5f %10.5f %10.5f %10.5f\n', DeltaAMultMax, DeltaAMultMean, DeltaAMultStd, TStateAMult); fprintf(1,'Метод Рунге-Кутта 4п%10.5f %10.5f %10.5f %10.5f\n', DeltaRC4Max, DeltaRC4Mean, DeltaRC4Std, TStateRC4); fprintf(1,'\n');