Runge-Kutta method for approximating differential equations
Methods for approximating the solutions of differential equations are useful, as most differential equations cannot be solved analytically. Sometimes it is extremely difficult to integrate the functions or provide explicit solutions.
We may then invoke approximating methods.
Such methods include the Euler method, Midpoint method and the Runge-Kutta method. The Euler method is the simplest:
It approximates the unknown function by a straight line i.e. Y=m(X-Xo) + Yo, where m = dy/dx. which is the gradient. We know from definition that dy/dx = f(x,y) (some function) which leads us to write
Y=f(Xo,Yo)(X-Xo)+Yo. So we start off at some initial point values and head off in a direction.
However we notice that the approximation would keep going in a straight line and it would not 'stick' to the true plot so we invoke an iterative method.
Yn+1= h f(Xn,Yn) + Yn
You notice that the Euler method is essentially too naive. We invoke a second method which uses the midpoint of the midpoint of the Euler method. So it is essentially correcting itself by 1 step which is a better approximation.
We then... you guessed it.... do more midpoints. We end up using a very good approximation which is known as the R-K approx.
The fundamental algorithm is:
Yn+1= Yn + k1/6+k2/3+k3/3+k4/6+O (h5)
where k1 = hf(Xn,Yn) , k2 = hf(Xn+1/2h, Yn + 1/2k1) , k3 =hf(Xn +1/2h, Yn + 1/2k2),
k4 = hf(Xn+h, Yn + k3)
Where K's are just co-efficients which aid in self correcting the algorithm in 'figuring' out the new gradient after an iteration.
We can apply this method in approximating SECOND-order differential equations if we write them in terms of first order differential equations. (If said differentials are coupled... meaning they depend on one another).
How to use coupled differentials to approx orbits
The fundamental algorithm of the R-K method is always the same. We just need to formulate 4 differential equations (2 for each dimension (X, Y)).
We can say define (X,Y,Velocity X, Velocity Y) as quantities relating to Z(1), Z(2), Z(3) and Z(4) respectively.
We say X and Y are the distances the Earth's centre is from the origin (the sun).
From general laws of classical mechanics: dX/dt = Velocity X, dY/dt = Velocity Y, d(Velocity X)/dt = acceleration X direction and d(Velocity Y)/dt = acceleration in Y direction.
This means we have reduced d^2 X/dt^2 (second order) as first order differential.
We can set the acceleration X and acceleration Y in accordance with the laws of Keplar and Newton:
- Ms*X/(R^2) (For X) and -Ms*Y/(R^2) (For Y).... In cartesian co-ords : R = Root (X^2 + Y^2)
but we have to times R^2 by another R as we are getting the magnitude (Basically using the definitions of unit vectors and changing the unit vector in the Rth direction as a value R)
So we finally end up... with our Z(4) arrays which relate to the x-y displacements and velocities of the earth
a sub function which looks like:
Dim Mass_Sun, Mass_Earth As Double
Mass_Sun = 100 : Mass_Earth = 2
dzdt(1) = z(3) '= Vel_X
dzdt(2) = z(4) '= Vel_Y
dzdt(3) = -(Mass_Sun * z(1)) / (z(1) ^ 2 + z(2) ^ 2) ^ (3 / 2) 'Accel earth X
dzdt(4) = -(Mass_Sun * z(2)) / (z(1) ^ 2 + z(2) ^ 2) ^ (3 / 2) 'Accel earth Y
With these in place and with the CORRECT initial conditions and with some program which calls plots of points... you can simulate the earth orbiting around the sun.
You can also do this with the moon. Setting the Z(4) array to an Z(8) array as we need 4 more differential functions to map the moons orbit onto some plot.
Thinking about this... the moon is orbiting around the earth but it also feels an attractive force from the sun. So its acceleration will contain one force term from contribution from the Sun's mass and another for the earth's mass. To plot it in relation to the SUN's origin we need to think about the differences in distances between earths points and moons and use these distances to define a new radius from the earths centre to the moons but still in the co-ordinate system of the sun.
We get a function like:
dzdt(5) = z(7)
dzdt(6) = z(8)
dzdt(7) = -Sun contribution in X - Earth contribution in X
dzdt(8) = '' '' for Y...
Note we have defined the Z array elements in a similar fashion to what we did with the earths.
Initial condition trial and error hence realizing the sensitivity of the systems behavior
on small deviations to initial conditions
You should have a nice initial condition for the Earth for a stable orbit. The moons orbit is much more complicated and you will find that it takes a while to see why. If you set the distance between moon and earth as too large the moon will not feel enough force to orbit the earth it will orbit the sun or fly off (if it has initial velocities). I will not say the appropriated initial conditions. This spoils the exercise.
I will however comment on the fact that, if one changes the initial velocity in a direction by a small amount.. say 1 part in 10 (+- 0.1) the system will look similar to the one before however over time the orbit will look completely different... what happened to me was that it would orbit for a small time say 1/4 of the earths orbit round the sun... then it would somehow get closer and closer to the earths centre.... its period value increases and then it just flies off at a tangent.
The system depends strongly on initial conditions which is a character of chaotic behavior (butterfly effect). Two body gravitational systems are complicated and we see in other star systems that planets orbits vary greatly and the moons orbits differ greatly too compared to others.
If the earth was slightly lighter or slightly further away from the sun we could see very different orbits from the one we are in now (over a long period of time).
It was an interesting observation. It also highlights the fact that objects orbiting other objects are actually statistically rare. Objects can fly past planets (asteroids, planet debris etc.) but not come into orbit, it takes some special cases for things to start orbiting planets and stars. Our universe and galaxy is sooooo large however that the actual occurrences of orbits will increase.
No comments:
Post a Comment