F112019

Contents

System av differentialekvationer
Högre ordningens differentialekvationer
(Visualisering) En pendel

System av differentialekvationer

% Exempel: Volterra's ekvationer med Euler's metod
clear all; clf
f=@(t,u)[ 0.5*u(1)-0.3*u(1)*u(2)
         -0.2*u(2)+0.1*u(1)*u(2) ];

N=5000; h = (80-0)/N;  % 0 <= t <= 80
t = linspace(0,80,N+1);
U(:,1)=[5;3];
for i = 1:N
    U(:,i+1) = U(:,i) + h*f(t(i),U(:,i));
end
plot(t,U);
legend('kaniner','rävar')

% Volterra's ekvationer med ode45
clear all, clf
f=@(t,u)[ 0.5*u(1)-0.3*u(1)*u(2)
         -0.2*u(2)+0.1*u(1)*u(2) ];
[t,U] = ode45(f,[0,80],[5;3]);
plot(t,U);
legend('kaniner','r�var')

Högre ordningens differentialekvationer

% Exempel: Pendelr�relse
l = 0.1; g = 9.81;
f=@(t,u)[u(2); -g/l*sin(u(1))];
[t U]=ode45(f,[0,2],[pi/4;0]);

clf;
subplot(2,2,1);
plot(t,U(:,1)); title('Vinkeln'); xlabel('t')
subplot(2,2,2);
plot(t,U(:,2));title('Vinkelhastigheten'); xlabel('t')
subplot(2,2,3);
plot(U(:,1),U(:,2)); title('Fasporträtt');
xlabel('\theta(t)')

(Visualisering) En pendel

subplot(2,2,4)
U(:,1)=U(:,1)-pi/2;
X=l*cos(U(:,1));Y=l*sin(U(:,1));

staang=plot([0 X(1)],[0 Y(1)]);hold on
klaepp=plot(X(1),Y(1),'ro');
axis equal; axis([-1 1 -2 1]/10);
axis manual
pause(5);

for i=2:numel(t)
    staang.XData=[0 X(i)];
    staang.YData=[0 Y(i)];
    klaepp.XData = X(i);
    klaepp.YData = Y(i);
    pause((t(i)-t(i-1))*2);
end