Implizite Runge Kutta DGL-System

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
Antworten
Pofi
User
Beiträge: 1
Registriert: Mittwoch 20. Mai 2020, 18:44

Hallo Leute,

bin seit Corona neu in der Programmierwelt von Python.

Erstes Problem:

Die von mir angenäherte Lösung einer DGL mit Impl. RK (Radau-IIA-2) entspricht nicht der exakten Lösung die zur Überprüfung bereits gegeben ist. Im Plot kann man sehen, dass sich die angenäherte Lösung wie die exakte verhält nur um einen gewissen Faktor 2 eben gestretcht wenn man das so beschreiben kann.

Hoffe, dass jmd von euch mal das selbe oder ähnliches Problem hatte und mir helfen kann.

Code: Alles auswählen

import numpy as m, numdifftools as nd, matplotlib.pyplot as p

# DGL System Funktion
 
def f(t,y):
    return m.array([y[1],-y[0]*m.pi**2])
   
# Aufstellung des nichtlinearen GLS und    

def psi1(k):
    return f(t + h/3 , y + h * (k[0] * 5/12 - k[1] * 1/12 ))
         
def psi2(k):
    return f(t + h , y + h * (k[0] * 3/4 + k[1] * 1/4))

def Psi(k):
    return m.array([psi1(k) , psi2(k)])

def DPsi(k):
    f_jacob = nd.Jacobian(f)
    return h * m.array([ [5/12 , -1/12],
                         [3/4  ,  1/4 ] ]) * f_jacob(t,y)

def delta_k(k):  
    return m.linalg.inv(DPsi(k) - m.array([ [1 , 0],
                                            [0 , 1 ] ])) * (Psi(k) - k)

# Startwerte                                              
                                                      
t=0

h=0.02

y=m.array([1 , 0])

k= f(t,y)

tol=0

tList=[]
y1List=[]
y2List=[]

# Newton und Iteration

while t < 2:
    
    k0  = k
    k   = k - delta_k(k)
    err = m.linalg.norm(k-k0)

    if err <= tol:
        t = t + h
        y   = y + h * (3/4 * k[0] + 1/4 * k[1])
        tList.append(t)
        y1List.append(y[0])
        y2List.append(y[1])

    
p.plot(tList,y1List,"r+")
p.plot(tList,y2List,"r+")

# Exakte Lösung

tList_ex= m.arange(0,2,0.02)

yex1 = m.cos(m.pi*tList_ex)
yex2 = -m.pi*m.sin(m.pi*tList_ex)

p.plot(tList, yex1 , "b-" )
p.plot(tList, yex2 , "b-" )
Antworten