Python: метод Рунге-Кутты четвертого порядка

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()

Итак, у меня закодирован метод runge kutta четвертого порядка, но часть, которую я пытаюсь вписать, заключается в том, что проблема говорит, что V_in (t) = 1, если [2t] четное, или -1, если [2t] нечетное.

Также я не уверен, что верну это уравнение: return 1 / RC * (V_in - V_out)

Вот в чем проблема:

Проблема 8.1

Буду очень признателен, если вы мне поможете!


python runge-kutta
person kassidylynnk    schedule 26.11.2015    source источник
comment
Я не знаю достаточно о проблемной области, чтобы должным образом помочь вам, но я почти уверен, что вам придется сгенерировать прямоугольную волну для ваших временных точек, которая представляет ваш Vin (аналогично тому, как вы генерируете xpoints). И используйте сгенерированный Vin (массив) в качестве входных данных в уравнение для генерации Vout (как получить правильное уравнение, я не знаю :) Мне придется вернуться к некоторым математическим вычислениям). Кажется, вы используете временные точки непосредственно в качестве входных данных, где, я думаю, вам нужно сначала сгенерировать квадратную волну и использовать ее в качестве входных ..   -  person Emile Vrijdags    schedule 26.11.2015

Ответы (2)


arrow_upward
1
arrow_downward

Итак, у меня закодирован метод runge kutta четвертого порядка, но часть, которую я пытаюсь вписать, заключается в том, что проблема говорит, что V_in (t) = 1, если [2t] четное, или -1, если [2t] нечетное.

Вы рассматриваете V_in как константу. Проблема говорит, что это функция. Итак, одно из решений - сделать это функцией! Это очень простая функция для написания:

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

Вам не нужен и не нужен этот if оператор в определении V_in(t). Ветвление внутри цикла стоит дорого, и эта функция будет вызываться много раз изнутри цикла. Есть простой способ избежать этого оператора if.

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

Эта функция настолько мала, что вы можете сложить ее в производную функцию:

def dV_out_dt(V_out, t) :
    return ((1 - 2*(math.floor(2.0*t) % 2)) - V_out)/RC
person David Hammen    schedule 26.11.2015

arrow_upward
0
arrow_downward

Функция должна выглядеть примерно так:

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)
person Lutz Lehmann    schedule 26.11.2015