Python: Fourth Order Runge-Kutta Method

2
votes
from math import sin
from numpy import arange
from pylab import plot,xlabel,ylabel,show
def answer():
    print('Part a:')
    print(low(x,t))
    print('First Graph')
    print('')


def low(x,t):
    return 1/RC * (V_in - V_out)

a = 0.0
b = 10.0
N = 1000
h = (b-a)/N
RC = 0.01
V_out = 0.0

tpoints = arange(a,b,h)
xpoints = []
x = 0.0

for t in tpoints:
    xpoints.append(x)
    k1 = h*f(x,t)
    k2 = h*f(x+0.5*k1,t+0.5*h)
    k3 = h*f(x+0.5*k2,t+0.5*h)
    k4 = h*f(x+k3,t+h)
    x += (k1+2*k2+2*k3+k4)/6

plot(tpoints,xpoints)
xlabel("t")
ylabel("x(t)")
show()

So I have the fourth order runge kutta method coded but the part I'm trying to fit in is where the problem say V_in(t) = 1 if [2t] is even or -1 if [2t] is odd.

Also the I'm not sure if I'm suppose to return this equation: return 1/RC * (V_in - V_out)

Here is the problem:

Problem 8.1

It would be greatly appreciated if you help me out!

2
pythonrunge-kutta
I don't know enough about the problem domain to properly help you, but I'm pretty sure you'll have to generate a squarewave for your time points which represents your Vin (similar to how you generate the xpoints). And use the generated Vin (array) as input to the equation to generate Vout (how to get the right equation I dont know :) I'll have to revisite some math). You seem to use the timepoints directly as input, where I think you have to first generate the squarewave and use that as input.. - Emile Vrijdags

2 Answers

1
votes

So I have the fourth order runge kutta method coded but the part I'm trying to fit in is where the problem say V_in(t) = 1 if [2t] is even or -1 if [2t] is odd.

You are treating V_in as a constant. The problem says that it's a function. So one solution is to make it a function! It's a very simple function to write:

def dV_out_dt(V_out, t) :
    return (V_in(t) - V_out)/RC

def V_in(t) :
    if math.floor(2.0*t) % 2 == 0 :
        return 1
    else :
        return -1

You don't need or want that if statement in the definition of V_in(t). A branch inside of a loop is expensive, and this function will be called many times from inside a loop. There's a simple way to avoid that if statement.

def V_in(t) :
    return 1 - 2*(math.floor(2.0*t) % 2)

This function is so small that you can fold it into the derivative function:

def dV_out_dt(V_out, t) :
    return ((1 - 2*(math.floor(2.0*t) % 2)) - V_out)/RC
0
votes

The function should look something like this:

def f(x,t):
    V_out = x
    n = floor(2*t)
    V_in = (1==n%2)? -1 : 1
    return 1/RC * (V_in - V_out)