[Tutor] hi

Oscar Benjamin oscar.j.benjamin at gmail.com
Mon Jul 22 13:59:07 CEST 2013


On 12 July 2013 11:08, Vick <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/

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


More information about the Tutor mailing list