Mathesis Wiki

Dies ist eine alte Version des Dokuments!

---

Protokoll 1.6.17:
Diffusionscode:

import matplotlib.pyplot as plt
import numpy as np
from copy import deepcopy

x0=0
x1=10
dx = 0.1
dt=0.0001
D= 10.0
def init():
	#sin
	sigma = 0.8
	a = np.linspace(x0,x1,(x1-x0)/dx)
	#sin
	#a = 10*np.sin(a)
	#gauss
	a = 1/sigma*np.exp(-0.5*((a-(x1-x0)/2)/sigma)**2)
	return a


def getnewtimestep(a):
	b=[0]*len(a)
	for i in range(1,len(a)-1):
		b[i] = a[i] + D*dt/dx**2*(a[i+1]+a[i-1]-2*a[i])
	b[0] = a[0] +D*dt/dx**2 * (a[1]-2*a[0])    #Randpunkte muessen gesondert behandelt werden
	b[-1] = a[-1] + D*dt/dx**2 * ( a[-2] - 2*a[-1])  #Randpunkte muessen gesondert behandelt werden
 	return b	




a=  init()
N = 10000

for i in range(N):
	if i % (N/10)==0:
		plt.plot(a, label="t= " + str(i))
	a = getnewtimestep(a)
plt.legend()
#plt.ylim([0,10])
plt.show()

Ansonsten: den Diffusionscode verändert (Betrag, Gerade, Anfangswerte verändert) game of life angefangen

Protokoll 08.06.17:

Lotta Volterra Pfeilmodell im iPython notebook :

import matplotlib.pyplot as plt
import numpy as np
def LotkaVolterra(N1,N2):
    eps1 = 6
    eps2 = 9
    gamma1 = 12
    gamma2 = 18

    return [ N1*(eps1-gamma1*N2),-N2*(eps2-gamma2*N1)]
    x = np.linspace(0,2,num =25)
y = np.linspace(0,2,num =25)
X,Y = np.meshgrid(x,y)
U,V = LotkaVolterra(X,Y)
plt.quiver(X,Y,U,V,linewidths=.1)
plt.axis('equal')
plt.ylabel("Jaeger")
plt.xlabel("Beute")

plt.grid()
plt.show()

Fitz-Hugh-Nagumo Modell (auch mit Pfeilen):

import matplotlib.pyplot as plt
import numpy as np

def Fitz_Hugh(a,b):
    alpha = 0.2
    beta = 5

    return [ a - a**3 - b + alpha, beta * (a -b)]
    
x = np.linspace(0,2,num =25)
y = np.linspace(0,2,num =25)
X,Y = np.meshgrid(x,y)
U,V = Fitz_Hugh(X,Y)
plt.quiver(X,Y,U,V,linewidths=0.1)
plt.axis('equal')
plt.ylabel("Produkte")
plt.xlabel("Edukte")

plt.grid()
plt.show()

Lotka-Volterra mit Populationen-Zeit-Graph:

import matplotlib.pyplot as plt
import numpy as np
from copy import deepcopy

# r fuer beute
# b fuer raeuber
def lv(r0,b0):
	epsilonR = 1
	epsilonB = 1
	gammaR = 1
	gammaB = 1
	delta = 0.1
	r=[r0]
	b=[b0]
	akt_r=r0
	akt_b=b0
	for i in range (100):
		r.append(akt_r*(epsilonR-gammaR*akt_b)*delta+akt_r)
		b.append( -akt_b*(epsilonB-gammaB*akt_r)*delta+akt_b)
		akt_r = r[-1]
		akt_b = b[-1]
	return r, b
		
r,b=lv(3,3)	

plt.plot(b,label = "Raeuber")
plt.plot(r, label ="Beute")
plt.legend()
plt.show()

Fitz-Hugh-Nagumo Konzentrations-Zeit-Diagramm:

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
from copy import deepcopy

# a für stoff 1
# b für stoff 2
def fhl(a0,b0):
	delta = 0.1
	a=[a0]
	b=[b0]
	akt_a=a0
	akt_b=b0
	alpha = 1
	beta = 1
	for i in range (100):
		a.append(((akt_a - akt_a**3 -akt_b + alpha)*delta + akt_a))
		b.append(beta * (akt_a- akt_b) * delta + akt_b)
		akt_a = a[-1]
		akt_b = b[-1]
	return a, b
		
a,b=fhl(2,4)	

plt.plot(a,label = "Stoff1")
plt.plot(b, label ="Stoff2")
plt.legend()
plt.show()

Protokoll 15.6.17 Wir haben die zwei Teile der Diffusion-Reaktionsgleichung zusammengefügt. Dazu mussten wir den Diffusionsteil etwas abändern. Zuvor haben wir immer direkt den neuen Wert von a erhalten, aber wir brauchen die Differenz zwischen zwei Werten (also zwischen a und a +1 Zeitschritt). (umbenannt und verkürzt zu laplacian) Auch den Reaktionsteil haben wir verändert (s. code)

import matplotlib.pyplot as plt
import numpy as np
from copy import deepcopy

x0 =  10
x1 = 10
dx =  1    #abstand der diskreten punkte
dt =  0.001
#D  = 10.0

def init(): 
	a = np.random.normal(loc = 0, scale = 0.05, size = 100)
	b = np.random.normal(loc = 0, scale = 0.05, size = 100)
	return a, b

def laplacian(a):
	return ( - 2 * a + np.roll(a,1,axis=0) + np.roll(a,-1,axis = 0)) * 1/dx**2

def Ra(a,b):
	alpha = -0.005
	return a-a**3-b+alpha

def Rb(a,b):
	beta =10
	return (a-b)*beta

def NextStep(a,b):
	delta_a =dt * ( 1*laplacian(a) + Ra(a,b))
	delta_b =dt*  ( 100*laplacian(b) + Rb(a,b))
	a += delta_a
	b += delta_b
	return np.array(a),np.array(b) 
	


a, b =  init()
N = 150   #Zeitschritte

for i in range(N):
	a,b = NextStep(a,b)
	if i % (N/10)==0:
		plt.plot(a,label="A", marker = "x")
		plt.plot(b, label="B " , marker = "o")
		plt.legend()
		plt.ylim(-1,1)
		plt.show()	
	a,b = NextStep(a,b)