[Math] PYTHON RK2 (Midpoint Method)

numerical methodsordinary differential equationspartial differential equationspythonrunge-kutta-methods

Good evening, I am writing code for a Numerical Analysis Project, but I am having difficulty iterating the RK2 (Midpoint Method) Correctly. I am using Python to do it, could anyone take a look at my code and see what I am doing wrong?

from my understanding when running the code something, is wrong with the array y that I create to put the iteration points into that will be later used for plotting.

Meaning in the array the initial point y[0] and y[1] are populated but after that, all values just go to zero. I'm not sure what is exactly wrong.

Here is my Code:

    import packages needed
    from __future__ import division
    `import numpy as np`
#===============================


def Midpoint_meth(def_fn, a, b, N, ya):

    """
        Test the Midpoint Method to solve
        initial value problem y'=f(t,y) with t in [a,b] and y(a) = ya.
        Step size h is computed according to input number of mesh points N
    """

    f = def_fn #input definining function

    h = (b-a)/N # developing step size h, from input values divided by N

    t = np.arange(a, b+h, h) #array intialized to hold mesh points t
    y = np.zeros((N+1,))      #array to hold Midpoint Method approximated y values

    y[0] = ya   #intial condition 

    #iterative method
    for  i in range(0, N):

        tau = t[i]      #current mesh point t 
        w = y[i]        #current value y(t)

        k1 = f(tau,w)
        h1 = h/2.0

        # next iteration using midpoint method 
        y[i + 1] = y[i] + h * f(tau + h1, w + (h1*k1) )


        return (t, y)

trying to show:

$$\ w_i=\alpha$$
$$\ w_{i+_1}= \ w_i + hf\biggl((\ t_i + {\frac h2}, \ w_i + {\frac h2} f(\ t_i, \ w_i)\biggr) $$

Best Answer

You need to be careful with indentations in python. It does not matter how many empty lines you put between them, lines with the same indentation level belong to the same code block.

In your case, the return statement is on the same indentation level as the body of the loop. Thus you insert y[0] on initialization, y[1] during the first execution of the loop body and directly return after that computation.

In general it helps to use a debugger or insert print statements (that you later remove) at important way points to follow the way of the execution. For instance in the loop 

print 'loop',i, t[i], y[i]
Related Question