Implizite Runge Kutta DGL-System
Verfasst: Dienstag 26. Mai 2020, 11:04
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.
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-" )