CODEWITHSUNDEEP
pythonmathmatplotliboderunge-kutta
bhjghjh
bhjghjh

Reputation: 917

Trying to solve second order DE through Euler and Runge_Kutta method

I am trying to solve a spring-mass problem with both Euler and Range-Kutta method and compare the plots. I have written functions for both Euler and Runge-Kutta but after calling the function to my problem , it seems my plot is not showing any data. please help me fix the plot and check if there is any error in my code, thanks

#function Euler
def euler ( y, t, dt, derivative):
    y_next = y + derivative(y, t) * dt
    return y_next

# function Runge-Kutta
# 2nd order Runge-Kutta method routine

def Runge_Kutta (y, time, dt, derivative):
    k0 = dt * derivative (y, time)
    k1 = dt * derivative (y + k0, time + dt)
    y_next = y + 0.5 * (k0 + k1)
    return y_next

and here is the problem I am trying to solve

[![""" A spring and mass  system. the coefficient of friction \mu is not negligible.generate a position vs. time plot for the motion of the mass, given an initial displacement x = 0.2m , spring constant k = 42 N/m , mass m =0.25 Kg, coefficient of friction \mu = 0.15 and initial velocity v = 0

F = -kx +/-mu mg """


from pylab import *
from Runge_Kutta_routine import Runge_Kutta
from eulerODE import euler

N = 500      #input ("How many number of steps to take?")
x0 = 0.2
v0 = 0.0
tau = 3.0     #input ("What is the total time of the simulation in seconds?")
dt = tau /float ( N-1)

k = 41.0      #input (" what is the spring constant?")
m = 0.25      #input ("what is the mass of the bob?")
gravity = 9.8
mu = 0.15     #input ("what is the coefficient of friction?")

""" we create a Nx2 array for storing the results of our calculations. Each 2- element row will be used for the state of the system at one instant, and each instant is separated by time dt. the first element in each row will denote position, the second would be velocity"""

y = zeros (\[N,2\])

y \[0,0\] = x0
y \[0,1\] = v0

def SpringMass (state, time):
    """ we break this second order DE into two first order DE introducing dx/ dt = v & dv/dt = kx/ m +/- mu g....

Note that the direction of the frictional force changes depending on the sign of the velocity, we handle this with an if statement."""


    g0 = state\[1\]
    if g0 > 0:
        g1 = -k/m * state \[0\] - gravity * mu
    else:
        g1 = -k/m * state \[0\] + gravity * mu

    return array (\[g0, g1\])

# Now we do the calculations
# loop only N-1 so that we don;t run into a problem addresssing y\[N+1\] on the last point

    for j in range (N-1):
        #y \[j+1\] = euler ( y\[j\] , 0, dt, SpringMass)
        y \[j+1\] = Runge_Kutta ( y\[j\], 0 , dt, SpringMass)

# Now we plot the result

time = linspace ( 0 , tau, N)
plot ( time, y\[:,0\], 'b-', label ='position')

xlabel('time')
ylabel('position')


show()][1]][1]

Upvotes: 3

Views: 959

Answers (1)

Rory Daulton
Rory Daulton

Reputation: 22564

It looks like your loop starting with for j in range (N-1): that calculates array y is indented, so Python considers those lines to be part of the function SpringMass. Since those lines come after the return statement, they are never executed.

To correct this, move those lines so the for line has no indentation and the other lines have only four-spaces of indentation. See if that solves your problem.

Note that your code as written here will still not work. You have extraneous backslashes before square brackets, you write your Euler and Runge_Kutta functions in a single unnamed module but the main code expects them to be in two different modules, etc. You also have many examples of bad style. Those problems are probably why you have not yet gotten any (other) answer. Do yourself a favor and clean your code before posting here and clean up your style.

Upvotes: 1

Related Questions