Система ODEs с использованием четвертого порядка Runge-kutta в Python

Question

Я пишу программу python для решения системы 2x2 дифференциальных уравнений первого порядка, учитывая оба начальных условия. Мой код:

from math import * 
import numpy as np 

np.set_printoptions(precision=6) ## control numpy's decimal output :) 

form1 = raw_input('Enter the 1st function of u & v only >> ') 
form2 = raw_input('Enter the 2nd function of u & v only >> ') 
start = input("Enter the lower limit of the interval >> ") 
stop = input("Enter the upper limit of the interval >> ") 
h = input("Using step size? ") 

N = int(np.round((stop - start)/h, decimals = 4)) ##calculate the  number of times to iterate 
## np.round fixes a python bug returning 2.99999997 instead of 
## 3 which reduced the number of times iterated by -1 
k = np.zeros((2, 4)) 

u = np.zeros((N +1,)) ## fill our u's first with N+1 0's 
v = np.zeros((N +1,)) ## fill our v's first with N+1 0's 


u[0] = input("Give me an initial value for u'>> ") 
v[0] = input("Give me the second initial value for v' >> ") 
t = np.arange(start, stop + h, h) 


def f1(t, u, v): 
    return eval(form1) 

def f2(t, u, v): 
    return eval(form2) 
##for u now 

def trialu(): 
    for j in range(0, N): 
     k[0, 0] = h * f1(t[j], u[j], v[j]) 
     k[1, 0] = h * f2(t[j], u[j], v[j]) 
     k[0, 1] = h * f1(t[j] + h/2, u[j] + 0.5*k[0, 0], v[j] + 0.5*k[1, 0]) 
     k[1, 1] = h * f2(t[j] + h/2, u[j] + 0.5*k[0, 0], v[j] + 0.5*k[1, 0]) 
     k[0, 2] = h * f1(t[j] + h/2, u[j] + 0.5*k[0, 1], v[j] + 0.5*k[1, 1]) 
     k[1, 2] = h * f2(t[j] + h/2, u[j] + 0.5*k[0, 1], v[j] + 0.5*k[1, 1]) 
     k[0, 3] = h * f1(t[j], u[j] + k[0, 2], v[j] + k[1, 2]) 
     k[1, 3] = h * f2(t[j], u[j] + k[0, 2], v[j] + k[1, 2]) 
     u[j+1] = u[j] + (k[0, 0] + 2*k[0, 1] + 2*k[0, 2] + k[0, 3])/6 
     v[j+1] = v[j] + (k[1, 0] + 2*k[1, 1] + 2*k[1, 2] + k[1, 3])/6 
    return u 

Я знаю, что могу return (u, v), но я хочу print "u ~ ", u по-разному, а также v на другой линии, но мне кажется, что я просто должен повторить trialu вернуться V, Есть ли лучший способ? без повторения def trialu()? Я также хочу напечатать k для j = 1, 2, ... отдельно?
Примером является решение u' = -3*u + 2*v и v' = 3*u - 4*v от [0,0.4] с использованием шагового размера h = 0.1 и u[0] = 0 и v = 0.5.

Answer 1

Я не уверен, что вы спрашиваете, может быть, это полезно:

# example u and v 
u = [1.003, 1.002, 1.001] 
v = [0, 0, 0] 

for i, (uapprox, vapprox) in enumerate(zip(u, v)): 
    print "on iteration", i, "u ~", uapprox, "and v ~", vapprox 

Выход:

on iteration 0 u ~ 1.003 and v ~ 0 
on iteration 1 u ~ 1.002 and v ~ 0 
on iteration 2 u ~ 1.001 and v ~ 0