Численные методы интегрирования и оптимизации сложных систем
7. Приложение 1 (Листинг скриптов для нахождения корней полинома)
function secush
clc
e=10.^-5;
x=-8.1;
xm1=-8
Asm1=8.6159999
i=0;
As=0.252*(x.^3)+1.41*(x.^2)+14.2*x+161;
x1=x-(As.*(xm1-x))./(Asm1-As);
Asm1=As;
As=0.252*(x1.^3)+1.41*(x1.^2)+14.2*x1+161;
i=i+1;
while abs(x1-x)>e
xm1=x;
x=x1;
x1=x-(As.*(xm1-x))./(Asm1-As);
Asm1=As;
As=0.252*(x1.^3)+1.41*(x1.^2)+14.2*x1+161;
i=i+1;
A(i)=x;
end
hold on
for n=1:i
plot(n,A(n),'b-o')
end
grid on
xlabel('iteraciya')
ylabel('roots')
disp('ответ');
disp(x);
8. Приложение 2 (Листинг скриптов для решения ДУ)
function Difer
clc
T=4;
a0=638.89;
a1=56.35;
a2=5.60;
b0=595.24;
h=0.0005;
A_X(1,1:3)=[0 0 0];
A=[0 1 0;
0 0 1;
a0 a1 a2];
B=[0 0 b0]';
k=0;
t=0;
while (t < (T-h))
if (t <= 3*h)
K1=A*(A_X(k+1,:))';
K2=A*(A_X(k+1,:))'+1/3*K1;
K3=A*(A_X(k+1,:))'+1/6*K1+1/6*K2;
K4=A*(A_X(k+1,:))'+1/8*K1+3/8*K2;
K5=A*(A_X(k+1,:))'+1/2*K1-3/2*K3+2*K4;
A_X(k+2,:)=(A_X(k+1,:))+h/6*(K1'+4*K4'+K5');
else
h1=h;
t=t+h1;
H=(eye(length(A_X(1,:)))-(9*h1/24)*A);
G=(eye(length(A_X(1,:)))+19*h1/24*A)*(A_X(k+1,:))'+h1/24*A*(-5*(A_X(k,:))'+(A_X(k-1,:))')
+h1/24*B*(9*1+19*1-5*1);
A_X(k+2,:)=(inv(H)*G)';
end
Otr(k+1)=t;
k=k+1;
h=-0.43496
end
plot(Otr,A_X(1:k,1),'b-');
grid on
9. Приложение 4 (Листинг скриптов для спектрального анализа)
spectr.m
syms t T;
Kx=(638.89/2)*(t-T).^2-56.35*(1./2)*(-2*(t-T))+5.6;
Ky=(595.24/2)*(t-T).^2;
F2=2*t;
L(2)=F2;
F3=4*t.^2-1;
L(3)=F3;
F4=8*t.^3-4*t;
L(4)=F4;
F5=16*t.^4-12*t.^2+1;
L(5)=F5;
F6=32*t.^5-32*t.^3+6*t;
L(6)=F6;
F7=64*t.^6-80*t.^4+24*t.^2-1;
L(7)=F7;
F8=128*t.^7-192*t.^5+80*t.^3-8*t;
L(8)=F8;
F9=256*t.^8-448*t.^6+240*t.^4-40*t.^2+1;
L(9)=F9;
F10=512*t.^9-1024*t.^7+672*t.^5-160*t.^3+10*t;
L(10)=F10;
F1=1;
L(1)=F1;
F2=2*T;
L1(2)=F2;
F3=4*T.^2-1;
L1(3)=F3;
F4=8*T.^3-4*T;
L1(4)=F4;
F5=16*T.^4-12*T.^2+1;
L1(5)=F5;
F6=32*T.^5-32*T.^3+6*T;
L1(6)=F6;
F7=64*T.^6-80*T.^4+24*T.^2-1;
L1(7)=F7;
F8=128*T.^7-192*T.^5+80*T.^3-8*T;
F9=256*T.^8-448*T.^6+240*T.^4-40*T.^2+1;
L1(9)=F9;
F10=512*T.^9-1024*T.^7+672*T.^5-160*T.^3+10*T;
L1(10)=F10;
F1=1;
L1(1)=F1;
G=L'*L1;
In=Kx*G;
r=int(In,T,0,t);
Cx=int(r,t,0,1.5);
In=Ky.*G;
r=int(In,T,0,t);
Cy=int(r,t,0,1.5);
A=((Cx+eye(10)).^-1)*Cy;
Cy=int(L,t,0,1.5);
Cx=A*Cy'
Postr.m
function H=fun(t)
Cx=[3.7672; 1.3134; 0.5181; 0.2065; 0.0819; 0.0323; 0.0127; 0.0491; 0.0189; 0.0723];
F2=2*t;
L(2)=F2;
F3=4*t.^2-1;
L(3)=F3;
F4=8*t.^3-4*t;
L(4)=F4;
F5=16*t.^4-12*t.^2+1;
L(5)=F5;
F6=32*t.^5-32*t.^3+6*t;
L(6)=F6;
F7=64*t.^6-80*t.^4+24*t.^2-1;
L(7)=F7;
F8=128*t.^7-192*t.^5+80*t.^3-8*t;
L(8)=F8;
F9=256*t.^8-448*t.^6+240*t.^4-40*t.^2+1;
L(9)=F9;
F10=512*t.^9-1024*t.^7+672*t.^5-160*t.^3+10*t;
L(10)=F10;
F1=1;
L(1)=F1;
H=(Cx'*L');
t=[0:0.01:5];
plot(t,H)
10. Приложение 5 (Листинг скриптов для оптимизации)
jivs.m
clear
clc
a=0;
b=5;
h=0.1;
Kp1=2; Kd1=1; Ki1=0;
J1=int2(Kp1, Kd1, Ki1);
while h>0.0000001
Kp2=Kp1+h;
J2=int2(Kp2, Kd1, Ki1);
if J2>J1
Kp2=Kp1-h;
J2=int2(Kp2,Kd1,Ki1);
if J2>J1
Kp2=Kp1;
end
end
Kd2=Kd1+h;
J2=int2(Kp2, Kd2, Ki1);
if J2>J1
Kd2=Kd1-h;
J2=int2(Kp2,Kd2,Ki1);
if J2>J1
Kd2=Kd1;
end
end
Ki2=Ki1+h;
J2=int2(Kp2, Kd2, Ki2,h);
if J2>J1
Ki2=Ki1-h;
J2=int2(Kp2,Kd2,Ki2,h);
if J2>J1
Ki2=Ki1;
end
end
h=fibon(a,b,h);
while J2<J1
Kp=Kp1+2*(Kp2-Kp1); Kd=Kd1+2*(Kd2-Kd1); Ki=Ki1+2*(Ki2-Ki1);
J1=J2;
J2=int2(Kp,Kd,Ki,h);
Kp1=Kp2;Kp2=Kp; Kd1=Kd2;Kd2=Kd; Ki1=Ki2;Ki2=Ki;
end
end
disp(Kp)
disp(Kd)
disp(Ki)
int2.m
function J=int2(Kp,Kd,Ki,h)
clc
T=4;
A=[0 1 0 0; 0 0 1 0; 0 0 0 1; -595.23809523*Ki -43.6507936507-595.23809523*Kp -56.34920634920635-595.23809523*Kd -5.59523809523809];
B=[0; 595.23809*Kd; 595.23809*Kp-3330.498866*Kd; 595.23809*Ki-33540.615-354308.277*(Kd)^2-3330.498*Kp-18634.934*Kd];
k=0;
t=0;
while (t < (T-h))
if (t <= 3*h)
K1=A*(A_X(k+1,:))';
K2=A*(A_X(k+1,:))'+1/3*K1;
K3=A*(A_X(k+1,:))'+1/6*K1+1/6*K2;
K4=A*(A_X(k+1,:))'+1/8*K1+3/8*K2;
K5=A*(A_X(k+1,:))'+1/2*K1-3/2*K3+2*K4;
A_X(k+2,:)=(A_X(k+1,:))+h/6*(K1'+4*K4'+K5');
else
h1=h;
t=t+h1;
H=(eye(length(A_X(1,:)))-(9*h1/24)*A);
G=(eye(length(A_X(1,:)))+19*h1/24*A)*(A_X(k+1,:))'+h1/24*A*(-5*(A_X(k,:))'+(A_X(k-1,:))')
+h1/24*B*(9*1+19*1-5*1);
A_X(k+2,:)=(inv(H)*G)';
end
Otr(k+1)=t;
k=k+1;
end
grid on
fibon.m
function h=fibon(a,b,h)
F(1)=1; F(2)=1;n=100;
for i=[1:0.1:n-2]
F(i+2)=F(i+1)+F(i);
end
j=0;
x1=a; x3=b;
L1=x3-x1;
L2=(F(n-1)/F(n))*L1+((-1)^n)/F(n)*eps;
x2=x3-L2;
x4=x1+x3-x2;
while (abs(x3-x1) > eps)
F2=x2;
F4=x4;
if ((x2 < x4)&&(norm(F2) < norm(F4)))
x1=x1; x3=x4;
x4=x1+x3-x2;
elseif ((x2 > x4)&&(norm(F2) < norm(F4)))
x1=x4; x3=x3;
x4=x1+x3-x2;
elseif ((x2 < x4)&&(norm(F2) > norm(F4)))
x1=x2; x3=x3;
x2=x1+x3-x4;
elseif ((x2 > x4)&&(norm(F2) > norm(F4)))
x1=x1; x3=x2;
x2=x1+x3-x4;
end
j=j+1;
la=x1+(x3-x1)/2;
end
l=la;