ODE matlab

[samcam]

beste,
ik heb in matlab een exacte waarde 'y1' deze probeer ik te benaderen met euler en met runge-kutta. de runge-kutta methode zou nauwkeuriger moeten zijn dan euler, dit is bij mij echter niet het geval. heeft iemand een idee wat ik hier fout doe?

mvg. sam


clear all, close all, clc

dt = 0.10;
t=0:dt:11;
A=0.001;
B=0.000577;

y(1)=0.001;
z(1)=0;
y2(1)=0.001;
y1 =(exp(-0.5.*t).*(A.*cos((sqrt(3)/2).*t)+B.*sin((sqrt(3)/2).*t))); %exacte waarde

a(1)=0;
y3(1)=0.001;

for n=1:110 %euler
y(n+1)=y(n)+z(n).*dt;
z(n+1)=z(n)-(y(n)+z(n)).*dt;
end



for n=1:110 %runge-kutta
k1=a(n).*dt;
zk1=-(y2(n)+a(n)).*dt;
k2=(a(n)+0.5.*zk1).*dt;
zk2=(-(y2(n)+a(n))+0.5.*k1).*dt;
k3=(a(n)+0.5.*zk2).*dt;
zk3=(-(y2(n)+a(n))+0.5.*k2).*dt;
k4=(a(n)+zk3).*dt;
zk4=(-(y2(n)+a(n))+k3).*dt;

y2(n+1)=y2(n)+1/6.*(k1+2.*k2+2.*k3+k4);
a(n+1)=a(n)+1/6.*(zk1+2.*zk2+2.*zk3+zk4);

end

plot(t,y1,'r',t,y,'g',t,y2,'b')

[arie]

Hoe ziet je ODE er precies uit?

Hieronder heb ik je exacte functie iets anders geschreven:
y = R .* exp(-0.5 .* tvec) .* cos((sqrt(3)/2).*tvec - alpha);

Analytisch volgt hieruit:
dydt = -R * exp(-0.5 * t) * cos((sqrt(3)/2)*t - alpha - pi/3);
dit vindt je in de code hieronder als functie van t (de @(t) definitie).

Voor beide methoden (Euler en RK) heb ik tevens de som van de kwadraten van de verschillen tot y berekend.
Voor Euler kom ik uit op 2.0044e-011 (als ik dydx in het midden van elk interval bepaal)
Voor RK levert dit 6.0171e-021, een beduidend kleinere afwijking.

Code: Selecteer alles

clear all, close all, clc

dt = 0.10;
tvec = 0:dt:11;
A = 0.001;
B = 0.000577;

R = sqrt(A*A + B*B);
alpha = atan(B/A);


y = R .* exp(-0.5 .* tvec) .* cos((sqrt(3)/2).*tvec - alpha);

dydt = @(t) -R * exp(-0.5 * t) * cos((sqrt(3)/2)*t - alpha - pi/3);



% ---  euler:  ---
seu = 0;
eu(1) = A;
for n=1:110
eu(n+1) = eu(n) + dydt((n-1)*dt + dt/2)*dt;
seu = seu + (eu(n+1)-y(n+1))^2;
end

seu


% ---  runge-kutta:  ---
srk = 0;
rk(1) = A;
for n=1:110
k1 = dydt((n-1)*dt)*dt;
k2 = dydt((n-1)*dt+dt/2)*dt;
k3 = dydt((n-1)*dt+dt/2)*dt;
k4 = dydt((n-1)*dt+dt)*dt;

rk(n+1)=rk(n) + (k1 + 2*k2 + 2*k3 + k4)/6;
srk = srk + (rk(n+1)-y(n+1))^2;
end

srk

plot(tvec,y,'r', tvec,eu,'g', tvec,rk,'b');

[samcam]

de constanten zijn A = 0.001 en B = 0.000577

y'' + y' + y = 0

alvast bedankt,

met vriendelijke groet,

Sam

[arie]

OK.
Je analytische oplossing hiervan klopt.
Je state space klopt ook.
Het gaat mis binnen je Runge Kutta afhandeling:
In plaats van a(n)+0.5.*zk1 zou je a(n+0.5.*zk1) willen hebben, alleen dat lukt niet omdat a een matlab vector is.

Je kan dit oplossen door je waarden vanuit de state space te berekenen: z' = M * z
In dit geval:



Hieronder (in de code) zijn v en k1 t/m k4 vectoren, z1 gebruik ik om het resultaat te plotten.

De blauwe curve ligt nu vrijwel volledig over de rode.
De norm voor RK (6.7753e-007) is duidelijk beter dan die voor Euler (1.9038e-004), zoals verwacht.

Zocht je zoiets?

Code: Selecteer alles

clear all, close all, clc

dt = 0.1;
t=0:dt:11;

nmax = 11/dt;

A=0.001;
B=0.000577;

y1 =(exp(-0.5.*t).*(A.*cos((sqrt(3)/2).*t)+B.*sin((sqrt(3)/2).*t))); %exacte waarde


% ---  euler:  ---
y(1)=0.001;
z(1)=0;

for n=1:nmax
y(n+1)=y(n)+z(n).*dt;
z(n+1)=z(n)-(y(n)+z(n)).*dt;
end


% ---  runge-kutta:  ---
z1(1) = 0.001;
z2(1) = 0.0;

M = [0 1; -1 -1];

for n=1:nmax

v = [z1(n); z2(n)];
k1 = dt*(M*v);
k2 = dt*(M*(v+k1/2));
k3 = dt*(M*(v+k2/2));
k4 = dt*(M*(v+k3));

v = v + (k1+2*k2+2*k3+k4)/6;
z1(n+1) = v(1);
z2(n+1) = v(2);

end


%    exact    Euler   RungeKutta
plot(t,y1,'r',t,y,'g',t,z1,'b')
norm(y-y1)
norm(z1-y1)

[samcam]

Hey arie,

Het ziet er wel goed uit, bedankt voor de moeite en tot de volgende keer

Ik heb nog wel een vraagje. Waarvoor wordt de functie norm() gebruikt?

Met vriendelijke groet
Sam

[arie]

De norm heb ik gebruikt om te kijken hoe goed de beide benaderingen overeenkomen met de analytische oplossing.

De Euclidische norm van een n-dimensionale vector x is:



Dit is de gebruikelijk de lengte van een vector.
Voor het verschil van 2 vectoren y en z levert dit:



Hier zie je dat hoe meer vectoren y en z op elkaar lijken, hoe kleiner de norm is.
Als y identiek is aan z (dus y = z), dan is de norm van (y-z) gelijk aan nul.

Omdat we hier 2 methoden vergelijken:
- Euler ten opzichte van de analytische oplossing
- RK ten opzichte van de analytische oplossing
en de lengtes van de 3 vectoren gelijk zijn (nmax), is dit een handige foutmaat.

Zie ook:
http://en.wikipedia.org/wiki/Norm_%28mathematics%29

En hieraan gerelateerd:
- de RMS-error, nu met factor 1/n onder de wortel:
http://en.wikipedia.org/wiki/Root_mean_ ... uare_error
en
- de standaard deviatie, nu alle x waarden ten opzichte van het gemiddelde:
http://en.wikipedia.org/wiki/Standard_d ... m_variable