townini 11.11.11.11

11111111111111

function [ xprim ] = rlc( t,x )

w=1;
R1=0.02;
R2=0.5;
C=0.01;
L1=1;
L2=1.5;
M=1;

D=M/(L1*L2-M^2);

A = [-R1*L2*D/M R2*L2*D/L2  -L2*D/M
     R1*D       -R2*L1*D/M  D
     1/C        0           0];

B = [L2*D/M
     -D
     0];

% e = sin(w*t);
% if t<0
%     e=0;
% else
%     e=1;
% end
if t==10
    e=1/100;
else
    e=0;
end
% if t==0 || t==1 || t==2 || t==3 || t==4 || t==5 || t==6 || t==7 || t==8 || t==9 || t==10 || t==11 || t==12
%     e=1;
% else
%     e=0;
% end
% e=sin(t)+sin(2*t)+sin(5*t);
xprim = A*x + B.*e;
end

//////for i=1:3
    for k=4:length(e)-4;
        xprim(i,k-3)=(1/25)*(48*xprim(i,k))-36*xprim(i,k-1)+16*xprim(i,k-2)-3*xprim(i,k-3)+12*dt*f(i,k+1);
    end
end

--------------------------------------------------------
2222222222222

clear all;
close all;
dt=0.001;
xx(:,1)=[0 ; 0 ; 0];
% for i=1:24000
%     xx(:,i+1)=xx(:,i)+dt*rlc(0.001*i,xx(:,i));
% end
% plot(xx(1,:))

dt=0.001;
% xx=zeros(3,1);
[tt,tty]=ode45(@rlc,[0:0.001:24],xx(:,1));
figure, plot(tt,tty(:,1));
figure, plot(tt,tty(:,2));
figure, plot(tt,tty(:,3));
% %
% tty=tty';
xx(:,2)=tty(2,:);
xx(:,3)=tty(3,:);
xx(:,4)=tty(4,:);
for k=4:(24/dt)
    xx(:,k+1) = 1/25 * (48*xx(:,k) -36*xx(:,k-1) + 16*xx(:,k-2) - 3*xx(:,k-3) + 12*dt*rlc((k+1)*dt, xx(:,k)));
    xx(:,k+1) = 1/25 * (48*xx(:,k) -36*xx(:,k-1) + 16*xx(:,k-2) - 3*xx(:,k-3) + 12*dt*rlc((k+1)*dt, xx(:,k+1)));
end
% hold on;
% t=linspace(0,24,24/dt);
% % % tty=tty';
figure
hold on
plot(tt, xx(1,:), 'r');
plot(tt, tty(:,1), 'g');
figure, plot(tt, tty(:,1)-xx(1,:)');
hold off;
figure
hold on
plot(tt, xx(2,:), 'r');
plot(tt, tty(:,2), 'g');
figure, plot(tt, tty(:,2)-xx(2,:)');
hold off;
figure
hold on
plot(tt, xx(3,:), 'r');
plot(tt, tty(:,3), 'g');
figure, plot(tt, tty(:,3)-xx(3,:)');
hold off;

[max(abs(tty(:,1)-xx(1,:)')) max(abs(tty(:,1)-xx(1,:)'))/max(abs(tty(:,1)))]
[max(abs(tty(:,2)-xx(2,:)')) max(abs(tty(:,2)-xx(2,:)'))/max(abs(tty(:,2)))]
[max(abs(tty(:,3)-xx(3,:)')) max(abs(tty(:,3)-xx(3,:)'))/max(abs(tty(:,3)))]
% % % plot(t,tty(1,:), '.-r');
% % % plot(t,tty(2,:), '.-g');
% % % figure, plot(t,xx(3,:))
% % % figure, ode45('rlc', 24, [0, 0, 0]');