RK4 i matlab

2010-09-27, 18:30 #1

Medlem

Reg: Nov 2009

Inlägg: 197

Om jag vill tillämpa runge-kutte 4 på en andra ordningens diffekvation i matlab -- jag vill ha ut läget ur en ekvation för accelerationen (alltså integrera två gånger) -- hur gör jag lättast? Kan någon ge mig lite ledning?

Tack, och ledsen om frågeställningen är lite vag.

Citera

2010-09-27, 18:53 #2

Medlem

Reg: Jan 2008

Inlägg: 3 806

Citat:

Ursprungligen postat av Astronom

Om jag vill tillämpa runge-kutte 4 på en andra ordningens diffekvation i matlab -- jag vill ha ut läget ur en ekvation för accelerationen (alltså integrera två gånger) -- hur gör jag lättast? Kan någon ge mig lite ledning?

Tack, och ledsen om frågeställningen är lite vag.

Separera den som ett system av två första ordningens ODEs genom att sätta y = dx/dt.

Alltså, säg att vi har ekvationen t.ex.

d²x/dt² + p(x, t) dx/dt + q(x, t) x + r(x, t)= 0.

Om vi då sätter y = dx/dt, så kan vi skriva detta som systemet

dy/dt + p(x, t) y + q(x, t) x + r(x, t) = 0
dx/dt = y

vilket alltså man kan se som en diffekvation för vektorn (x, y), och sen lösa denna med Runge-Kutta.

Citera

2010-09-27, 19:23 #3

Medlem

Reg: Nov 2009

Inlägg: 197

Bingo, det är exakt det där jag vill. Jag ska fippla lite och se om jag fattar, annars skriver jag igen.

Tack så mycket sålänge.

Citera

2010-09-27, 19:45 #4

Medlem

Reg: Nov 2009

Inlägg: 197

Jag är inte riktigt på det klara över hur man applicerar runge kutta på en vektor. Vilken form skulle koden ta?

Kan posta min algoritm som den ser ut nu, om det hade hjälpt.

Kod:

 for n=1:(N-1)
     t(n+1)=t(n) + dt;
     k1=(g-1/2*C*A*rho/m*v(n)^2);
     k2=(g-1/2*C*A*rho/m*(v(n)+1/2*dt*k1)^2);
     k3=(g-1/2*C*A*rho/m*(v(n)+1/2*dt*k2)^2);
     k4=(g-1/2*C*A*rho/m*(v(n)+dt*k3)^2);
     v(n+1)=v(n)+1/6*dt*(k1+2*k2+2*k3+k4);     
 end

__________________
Senast redigerad av Astronom 2010-09-27 kl. 19:59.

Citera

2010-09-27, 20:33 #5

Medlem

Reg: Jan 2008

Inlägg: 3 806

Citat:

Ursprungligen postat av Astronom

Jag är inte riktigt på det klara över hur man applicerar runge kutta på en vektor. Vilken form skulle koden ta?

Kan posta min algoritm som den ser ut nu, om det hade hjälpt.

Kod:

 for n=1:(N-1)
     t(n+1)=t(n) + dt;
     k1=(g-1/2*C*A*rho/m*v(n)^2);
     k2=(g-1/2*C*A*rho/m*(v(n)+1/2*dt*k1)^2);
     k3=(g-1/2*C*A*rho/m*(v(n)+1/2*dt*k2)^2);
     k4=(g-1/2*C*A*rho/m*(v(n)+dt*k3)^2);
     v(n+1)=v(n)+1/6*dt*(k1+2*k2+2*k3+k4);     
 end

Tja, säg att vi har ekvationen

dv/dt = f(v, t).

Då så måste du ändra i koden så att v är en vektor. Lämpligast är väl att göra v till en matris där den n:te kolonnen är värdet på v vid steg n. DÅ skulle sista raden bli något som t.ex.

v(:, n+1) = v(:, n) + 1/6 * dt * (k1 + 2*k2 + 2*k3 + k4).

Sen så måste du ändra så att k1, k2, k3, k4 är kolonnvektorer, dessa ska alltså vara olika värden för den vektorvärda funktionen f, så exakt hur det blir beror på exakt vilken ekvation det handlar om.

Citera

2010-09-27, 20:36 #6

Medlem

Reg: Nov 2009

Inlägg: 197

Mhm, jag ska leka mer med det och se vad som händer. Du är till stor hjälp, tack.

Min ekvation är för övrigt

dv/dt = g - 1/(2m)*C*A*rho*v^2

__________________
Senast redigerad av Astronom 2010-09-27 kl. 20:43.

Citera

2010-09-27, 23:52 #7

Medlem

Reg: Nov 2009

Inlägg: 197

Jag har kollat mer på det nu, och jag vet vad jag vill göra, men jag har svårt att se hur jag ska lägga upp koden.

Skulle någon kunna visa hur RK4-algoritmen ser ut till exempel för dv/dt = g - 1/(2m)*C*A*rho*v^2, om jag vill ha ut både läget och hastigheterna?

Citera

2010-09-28, 01:59 #8

Medlem

Reg: Jan 2008

Inlägg: 3 806

Låt s(t) vara sträckan som färdats, s(0) = s0; v(t) vara hastigheten, v(0) = v0.

Diffekvationsystemet är alltså

s'(t) = v(t)
v'(t) = g - kv(t)²

där jag har tagit mig friheten att införa konstanten k = C*A*ρ/(2m).

Så, låt oss betrakta vektorn x = (s, v)^T. Då är

x' = (x_2, g - k(x_2)²)^T

Om vi ska Runge-Kutta-lösa detta i Matlab blir det t.ex.

Kod:

function dx_dt = f(x)
	% givet en vektor x, ger dx/dt för detta x
	dx_dt = [x(2); g - k*x(2)^2];
end;

% Lös från steg 1 till N, varje tidssteg dt
x = zeros(2, N);
x(:,1) = [s0; v0];

for i=1:N
     k1=f(x(i));
     k2=f(x(i)+1/2*dt*k1);
     k3=f(x(i)+1/2*dt*k2);
     k4=f(x(i)+dt*k3);
     x(:,n+1) = x(:, n) + (k1 + 2*k2 + 2*k3 + k4)/6;   
end;

Har inte kört koden, så den funkar troligtvis inte, men men.

Citera

RK4 i matlab

Stöd Flashback