[FIZYKA] Lab 3

KOD PONIŻEJ (midPoint) TO TEN SAM KOD CO W NECIE TYLKO TROCHĘ INACZEJ UŁOŻONY, ŻEBY BYŁ CZYTELNIEJSZY – NIE KOPIUJCIE JAK MACIE JUŻ Z NETA

// p - punkt, dt - czas od ostatnich obliczen ruchu [reguluje szybkosc symulacji]
void solveMidPoint(Punkt *p, float dt)
{
	// KOMENTARZ
	// wzory po przeksztalceniu na choc troche prostsza matematyke sa w 'LAB03 metody.pdf'
	// glownie trzeba sie wzorowac na kodzie z 'solveEuler'
	// i dorobic jego rozwiniecie w midPoint
	// i jeszcze dluzsze rozwiniecie w RK4
	// KOMENTARZ END


	// y = pozycja w przestrzeni [zwana w kodzie 'r']
	// t = predkosc i kierunek ruchu [zwana w kodzie 'v']

	// troche zmiennych do trzymania wynikow obliczen
	Wektor k1[2],k2[2];
	Wektor y_k[2];
	Wektor y_h[2];
	Wektor yout[2];
	
	// obliczanie k1 czyli: f(y , t)
	y_h[0] = p->r; // ustawiam 'y'
	y_h[1] = p->v; // ustawiam 't'
	derivatives(y_h, yout, p); // magicznie sie liczy
	k1[0] = yout[0]; // zapisuje wynik 'y' do kolejnych obliczen
	k1[1] = yout[1]; // zapisuje wynik 't' do kolejnych obliczen
	
	// obliczanie k2 czyli: f(y+(h/2)*k1 , t +  h/2)
	y_k[0] = y_h[0] + dt * 0.5 * k1[0]; // ustawiam 'y' uwzgledniajac wynik obliczen k1
	y_k[1] = y_h[1] + dt * 0.5 * k1[1]; // ustawiam 't' uwzgledniajac wynik obliczen k1
	derivatives(y_k, yout, p); // magicznie sie liczy
	k2[0] = yout[0]; // zapisuje wynik 'y' do kolejnych obliczen
	k2[1] = yout[1]; // zapisuje wynik 't' do kolejnych obliczen

	// obliczam nowa pozycje i predkosc 'punktu'
	// ze wzoru: y[n+1] = y[n] + k2 * h + O(h^3)
	// to ostatnie '+ O(h^3)' znacz tyle,
	// ze blad obliczen moze wyniesc maksymalnie 'h do 3' - dosc malo
	// ale tego nie ma jak uwzglednic w obliczeniach
	p->r = p->r + k2[0] * dt;
	p->v = p->v + k2[1] * dt;	
}

Kod RK4 jaki znalazłem w necie ma dość dużo błędów i naprawienie ich zajęło mi dużo czasu [o ile w ogóle poprawiłem je dobrze…], więc tym podzielę się dopiero jak oddam sam [możecie sami kombinować o co w tym chodzi lub oddać za tydzień za 50% pkt.].