[Tutor] hi

[Vick]

Hello,

 

Thanks for taking the time to reply!

 

However, I had already solved the problem a few days ago on my own.

 

The DOPRI 8(7)13 is the Dormand Prince 8th order with 13 stages method for
solving  <http://en.wikipedia.org/wiki/Ordinary_differential_equations>
ordinary differential equations. It is a member of the Runge-Kutta family of
ODE solvers. I have been able to reproduce it in Python with a
multi-variable capability. On a particular test equation the error of the
RK4 is about 1e-4 whereas the DOPRI 8(7)13 is about 1e-13.

 

I have then used the multi-variable DOPRI to compute for the n-body problem
using 1st order ODEs which requires you to have 4 equations for 1 planet. I
have already made a 10-body problem through an n-body solver in python and
using DOPRI as my integrator.

 

It works fine.

 

Thanks again for your time and consideration.

 

Regards

Vick

 

 

-----Original Message-----
From: Oscar Benjamin [mailto:oscar.j.benjamin at gmail.com] 
Sent: Monday, 22 July, 2013 15:59
To: Vick
Cc: Tutor at python.org
Subject: Re: [Tutor] hi

 

On 12 July 2013 11:08, Vick < <mailto:vick1975 at orange.mu>
vick1975 at orange.mu> wrote:

> Hi,

 

Hi Vick,

 

Sorry for the delayed response I didn't see this email until a while after
you sent it.

 

> I have written a code to perform a numerical solution to 1st order 

> Ordinary Differential Equations (ODEs). I have written codes for the 

> famous Runge Kutta 4th order and the DOPRI 8(7)13. However the codes 

> are for 1 variable only; viz. the "y" variable.

 

It would be good to show us the code for this in your email (not in an
attachment). Also the code should be complete: the code uses a function
rungkdp8 which is not shown. Where did that come from? Until you know what
you're doing you should test your integrator on a simpler set of equations.
Also when posting to a mailing list you should simplify the problem as much
as possible which would mean using simpler ODEs. In short read this:

 <http://sscce.org/> http://sscce.org/

 

The code attached is incorrect as it applies each stage of the Runge-Kutta
method to all time-steps sequentially before moving on to the next stage. It
should apply each stage for one time step and then use the output of that
for the next time-step. In other words you need to change the order of the
loops.

 

Here is an example of how to implement an Euler stepper that should
hopefully get you started on the right track. I've deliberately avoided
using numpy in this example even though it makes most of this easier. I did
this to focus on the basics of Python which you should learn first.

 

# integrator.py

 

nan = float('nan')

 

# indices for the two variables in the system POS = 0 VEL = 1 VARS = ['pos',
'vel']

 

# system has one parameter OMEGA

OMEGA = 10  # Frequency in Hz is OMEGA / 2pi

 

# function for the derivative of the system def shm(t, x, dxdt):

    '''Equation of motion for a mass in a spring with freq OMEGA'''

    dxdt[POS] = x[VEL]     # Python uses square brackets for indexing

    dxdt[VEL] = - OMEGA**2 * x[POS]

 

# Initial conditions for the problem

x0 = [nan] * 2

x0[POS] = 1

x0[VEL] = 0

 

# This is the function that you should replace with e.g. rk4 or dopri8 def
euler(t1, x1, t2):

    '''Take 1 Euler step from x1 at t1 to x2 at t2'''

    # Create empty lists for answer

    x2 = [nan] * len(x1)

    dxdt = [nan] * len(x1)

    # Call derivative function to fill in dxdt

    shm(t1, x1, dxdt)

    # Compute x2 and return

    dt = t2 - t1

    for n in range(len(x1)):

        x2[n] = x1[n] + dt * dxdt[n]

    return x2

 

def solve(ts, x0):

    '''Compute states corresponding to the times in ts

 

    x0 is the state at ts[0] (the initial condition).

    '''

    # Create an empty list of lists for the solution

    Xt = [x0[:]]  # The initial conditions

    # Loop through timesteps computing the next value

    for n in range(len(ts) - 1):

        Xt.append(euler(ts[n], Xt[n], ts[n+1]))

    return Xt

 

def print_line(*items):

    print(', '.join(str(i) for i in items))

 

def print_solution(ts, Xt, vars):

    print_line('t', *vars)

    for t, x in zip(times, Xt):

        print_line(t, *x)

 

# Numerical parameters for the solution

DT = 0.001

t0 = 0

T  = 1

times = [t0]

while times[-1] < T:

    times.append(times[-1] + DT)

 

# Solve and print

Xt = solve(times, x0)

print_solution(times, Xt, VARS)

 

 

> I have written another code in Excel VBA for DOPRI 8 for a 

> multi-variable capability. However I have not been able to reproduce 

> it in Python. I'm having trouble in making arrays or lists, I don't 

> know which is supposed to work.

> 

> I have attached my Excel VBA code for a multi-variable numerical 

> integrator in PDF format. This does work in Excel. I have also attached my
python code.

 

I'm not familiar with Excel VBA but I think that this code suffers from the
same problem that the loop over stages takes place outside the loop over
time-steps so I believe that it is incorrect.

 

> Can anyone help me out please?

 

I happen to know quite a lot about this subject so feel free to ask more
questions.

 

 

Oscar

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/tutor/attachments/20130722/300aafac/attachment.html>