Using finite difference in python

Question

I am trying to use Python with Numpy to solve a basic equation using the finite difference method. The code gives me the correct first value for a i.e it gives me a[1]; however, every other value after that is just zero? I don't know what I'm doing wrong because it obviously works for the first value so how do I fix this? Any ideas would be very helpful.

from numpy import *
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
from scipy.integrate import odeint

def solver(omega_m, dt):
    #t_0, H_0, a_0, dt, n, T; always the same; omega's change
    t_0 = 0.0004
    a_0 = 0.001
    H_0 = 1./13.7

    T = 13.7
    dt = float(dt)
    n = int(round((T - t_0)/dt))
    x = zeros(n+1)
    t = linspace(t_0, T, n+1)

    x[0] = a_0

    for i in range (0, n):
        x[i+1] = x[i] + (H_0 * ((omega_m)**(1./2.)) * ((x[i])**(-1./2.)) * dt)
        return x, t

a, t = solver(omega_m =1, dt=0.001)
print a, t

Moses Koledoye · Accepted Answer · 2017-03-17 22:35:42Z

4

Your function returns after the first iteration because your return statement is inside the for loop. You should dedent the return statement so your loop does not terminate prematurely:

for i in range (0, n):
    x[i+1] = x[i] + (H_0 * ((omega_m)**(1./2.)) * ((x[i])**(-1./2.)) * dt)
return x, t

answered Mar 17, 2017 at 22:35

Moses Koledoye

78.8k8 gold badges139 silver badges141 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

rcss92 Over a year ago

Thank you very much, really should have noticed that myself! But thats a big help, working as it should now.

rcss92 Over a year ago

Done now. Sorry didn't know how to do that.

Collectives™ on Stack Overflow

Using finite difference in python

1 Answer 1

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related