Pagina 1 van 1
ODE matlab
Geplaatst: 24 feb 2015, 11:19
door 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')
Re: ODE matlab
Geplaatst: 24 feb 2015, 22:18
door 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');
Re: ODE matlab
Geplaatst: 25 feb 2015, 18:19
door samcam
de constanten zijn A = 0.001 en B = 0.000577
y'' + y' + y = 0
alvast bedankt,
met vriendelijke groet,
Sam
Re: ODE matlab
Geplaatst: 26 feb 2015, 12:18
door 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)
Re: ODE matlab
Geplaatst: 26 feb 2015, 15:14
door 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
Re: ODE matlab
Geplaatst: 26 feb 2015, 18:22
door 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:
^2 + (y_2-z_2)^2 + ... + (y_n-z_n)^2})
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