Runge-Kutta numerisk løsning i python

[tool-nes]

Har et problem med å sette opp en algoritme for Runge-Kutta i python.
Er ikke sikker på om dette er rette stedet å spørre, men prøver..

Jeg har en pendel som består av en ball som henger i et tau som igjen henger ifra taket.
Jeg har en linje ifra det laveste stedet hvor ballen kan befinne seg og opp til der tauet er festet i tauet.
Vinkelen mellom denne linja og tauet er [tex]\theta[/tex]. og tauet er [tex]r_0[/tex] langt.

Jeg skal sette opp en algoritme som skal vise posisjonen til ballen i [tex]xy-planet[/tex] i de første 10sekundene.

hvordan kan jeg få til dette i python?

[tisstrange]

Du kan, helt generellt, skrive en programm som løser diffligninger med en 4.ordens runge-kutta. Bare se på http://en.wikipedia.org/wiki/Runge_kutta og implementer dette i en enkel for løkke. Jeg har en gammel kode, men det er nok lurt å skrive sånt selv. Du bruker bare masse tid på å skjønne koden min uansett...

I forhold til oppgaven din: Dette ser ut som en matematisk pendel, som har en difflingning som du fint kan løse analytisk (i hvert fall for små vinkler). Er du sikkert på at du ikke skal ha noe friksjon/andre bevegelser/krefter?

[Bogfjellmo]

Hvor langt har du kommet selv? Har du satt opp en difflikning du vi python skal løse? Hvordan ser den i så fall ut?

Som tisstrange sier, kan du løse denne analytisk for små utslag om du gjør en liten forenkling, men det mest interessante her er kanskje å implementere en løser i python...

Du bør i alle fall bruke en eksplisitt RK-metode, med implisitte må du ofte løse ikkelineære systemer, og å gjøre det med data er nesten et helt fagfelt i seg selv.

[tool-nes]

slik ser koden min ut nå:

from scitools.all import *
from numpy import *

k = 2000
L = 1
m = 0.1
g = array([0, 9.81])
r0 = array([cos(pi/6), sin(pi/6)])
v0 = array([0, 0])

time = 10
dt = 0.001
n = int(round(time/dt))

r = zeros((n, 2), float)
v = zeros((n, 2), float)
t = zeros(n, float)

r[0] = r0
v[0] = v0

for i in range(n-1):
rr = sqrt(r[0]**2 + r[1]**2)
a = -g - (k/m) * (rr - L) * (r/rr)
v[i+1] = v + a*dt
r[i+1] = r + v[i+1]*dt
t[i+1] = t + dt

plot(r[:,0], r[:,1])
xlabel("x [m]")
ylabel("y [m]")
axis("equal")


men dette git ikke riktig plot. Det plotet jeg får ut er en nesten full sirkel, mangler kanskje 30grader. I denne oppgaven så vet jeg at tauet skal fungere som en fjær, og med den k'en som er oppgitt til 200, vil fjæra slakkes og strammes, slik at det skal bli en del av en sirkel som er bølgete.
Vet ikke hvor jeg gjør galt

[Bogfjellmo]

Du har snora opphengt i (0,0), ikke sant? Nå begynner du med positiv y, jeg vet ikke om det var meningen. For det andre er snora veldig stiv, så dette bør bli som på en sirkel, tilnærmet. Skal du få mer framtredende bølger bør du prøve med mindre k.

For det tredje ville jeg personlig brukt polarkoordinater, men det er ikke så farlig.

[tool-nes]

snora henger ifra et punkt i "taket", snora danner vinkelen [tex]\theta[/tex] med en vertikal linje fra det punktet snora henger. Snora skal starte med [tex]\theta = 30grader[/tex]. Er litt usikker på om snora er opphengt i (0,0). Tror egentlig ikke det. Står ikke noe om det i oppgaven. k er oppgitt i oppgaven som 200, en annen oppgave er å prøve med k=2000, slik at vi skal se forskjellen. Altså at snora strekkes og strammes mindre enn med k=200.

[Bogfjellmo]

Sånn som du har skrevet koden din, har du snora opphengt i (0,0), og kula starter i [tex](\cos(\frac \pi 6), \sin(\frac \pi 6)[/tex], som danner en vinkel [tex]30^\textdegree[/tex] med horisontalen, altså inne i taket, slik som du beskriver det. Jeg foreslår at du tegner en figur for å finne ut hvilke startverdier du egentlig vil ha.

Koden ser ellers helt fin ut, men du har ikke helt oversikt, virker det som. Har du skrevet koden selv, eller har du kopiert fra noe?

[tool-nes]

oki.. hvis jeg bytter om sin og cos i startverdien så blir det 30grader med vertikalen. godt mulig jeg har bytta om de.. figuren er tegna.
det er programmet jeg sliter med..