diff --git a/src/bvpy/scheme/explicit_rk.py b/src/bvpy/scheme/explicit_rk.py
index 68658b80a9bb832dfaf897cffe90063074357c07..3907b60b144c15e75f75ef0befb617b9640fee00 100644
--- a/src/bvpy/scheme/explicit_rk.py
+++ b/src/bvpy/scheme/explicit_rk.py
@@ -1,4 +1,4 @@
-class ExplicitButcherTableau():
+class ExplicitRKMethods():
def __init__(self, table):
self.table = table
@@ -11,6 +11,12 @@ class ExplicitButcherTableau():
self.extract_wheights()
self.extract_rk_matrix()
+ assert table[0][0] == 0, 'Bad format for the butcher tableau'
+ assert len(self.nodes) == len(self.wheights), 'Bad format for the butcher tableau'
+ for ind, a in enumerate(self.rk_matrix):
+ assert len(a)==ind+1, 'Bad format for the butcher tableau'
+
+
def extract_nodes(self):
self.nodes = [i[0] for i in self.table[:-1]]
self.s = len(self.nodes)
@@ -50,21 +56,77 @@ class ExplicitButcherTableau():
def step(self):
k=[]
for i in range(len(self.nodes)):
- print('compute k:', i)
sum_kj = [self.rk_matrix[i-1][j]*k[j] for j in range(i)]
ki = self.func(self.time+self.nodes[i]*self.dt, self.state +self.dt*sum(sum_kj))
k.append(ki)
- print('k',k)
res = self.state + self.dt*sum([b*k[ind] for ind, b in enumerate(self.wheights)])
self.time+=self.dt
+ self.state = res
- print('------>', res)
+# TODO: Test that
+class ImplicitRKMethods():
+ def __init__(self, table):
- self.state = res
+ self.table = table
+ self.nodes = None
+ self.rk_matrix = None
+ self.wheights = None
+ self.s = None
+ self.extract_nodes()
+ self.extract_wheights()
+ self.extract_rk_matrix()
+ def extract_nodes(self):
+ self.nodes = [i[0] for i in self.table[:-1]]
+ self.s = len(self.nodes)
+
+ def extract_wheights(self):
+ self.wheights = self.table[-1]
+
+ def extract_rk_matrix(self):
+ self.rk_matrix = []
+ for line in self.table[0:-1]:
+ self.rk_matrix.append(line[1:])
+
+ def show(self):
+ for ind_line, line in enumerate(self.table):
+ if ind_line == len(self.table)-1:
+ for i in range(len(line)):
+ print('-'*8+'--', end='')
+ print()
+ print(' |\t', end='')
+
+ for ind, val in enumerate(line):
+ print("{:.2f}\t".format(val), end='')
+ else:
+ for ind, val in enumerate(line):
+ if ind == 0:
+ print("{:.2f}|\t".format(val), end='')
+ else:
+ print("{:.2f}\t".format(val), end='')
+ if ind_line != len(self.table)-1: print()
+
+ def set_initial_state(self, state, dt, t, func):
+ self.time=0
+ self.state=state
+ self.dt=dt
+ self.func=func
+
+ def step(self):
+ k=[]
+ for i in range(len(self.nodes)):
+ sum_kj = [self.rk_matrix[i-1][j] for j in range(self.s-1)]
+ ki = self.func(self.time+self.nodes[i]*self.dt, self.state +self.dt*sum(sum_kj))
+ k.append(ki)
+ print(k)
+
+ res = self.state + self.dt*sum([b*k[ind] for ind, b in enumerate(self.wheights)])
+ self.time+=self.dt
+ self.state = res
+
if __name__ == "__main__":
from math import tan
@@ -82,31 +144,26 @@ if __name__ == "__main__":
[1/6,1/3,1/3,1/6]]
- table_RK4 = ExplicitButcherTableau(RK4)
+ table_RK4 = ExplicitRKMethods(RK4)
- table_euler = ExplicitButcherTableau(euler)
+ table_euler = ExplicitRKMethods(euler)
- table_ralston = ExplicitButcherTableau(ralston)
+ table_ralston = ExplicitRKMethods(ralston)
def f(t, y):
return tan(y) + 1
-
table_ralston.set_initial_state(1, 0.025, 1, f)
table_ralston.step()
+ print(table_ralston.state)
table_ralston.step()
+ print(table_ralston.state)
table_ralston.step()
+ print(table_ralston.state)
table_ralston.step()
+ print(table_ralston.state)
- print('------------------')
- print('------------------')
- print('------------------')
-
- table_RK4.set_initial_state(1, 0.025, 1, f)
-
- table_RK4.step()
- table_RK4.step()
- table_RK4.step()
- table_RK4.step()
+ backward_euler = [[1, 1],[1]]
+ table_be = ImplicitRKMethods(backward_euler)
diff --git a/src/bvpy/ivp.py b/src/bvpy/scheme/ivp.py
similarity index 100%
rename from src/bvpy/ivp.py
rename to src/bvpy/scheme/ivp.py
diff --git a/src/bvpy/scheme/test.py b/src/bvpy/scheme/test.py
new file mode 100644
index 0000000000000000000000000000000000000000..0504c1c763f329c5cdd9ba6366e400bf4cfe4ab4
--- /dev/null
+++ b/src/bvpy/scheme/test.py
@@ -0,0 +1,55 @@
+from fenics import *
+import time
+
+T = 2.0 # final time
+num_steps = 50 # number of time steps
+dt = T / num_steps # time step size
+
+# Create mesh and define function space
+nx = ny = 30
+mesh = RectangleMesh(Point(-2, -2), Point(2, 2), nx, ny)
+V = FunctionSpace(mesh, 'P', 1)
+
+# Define boundary condition
+def boundary(x, on_boundary):
+ return on_boundary
+
+bc = DirichletBC(V, Constant(0), boundary)
+
+# Define initial value
+u_0 = Expression('exp(-a*pow(x[0], 2) - a*pow(x[1], 2))',
+ degree=2, a=5)
+u_n = interpolate(u_0, V)
+
+# Define variational problem
+u = TrialFunction(V)
+v = TestFunction(V)
+f = Constant(0)
+
+def func(t, state):
+ F = dot(grad(state), grad(v))*dx - f*v*dx
+ return F
+
+
+# Create VTK file for saving solution
+vtkfile = File('heat_gaussian/solution.pvd')
+
+# Time-stepping
+
+t = 0
+form = (u-u_n)/dt*dx - func(t+0.5*dt, u_n + dt*func(t, u_n))
+a, L = lhs(F), rhs(F)
+u = Function(V)
+for n in range(num_steps):
+
+ # Update current time
+ t += dt
+
+ # Compute solution
+ solve(a == L, u, bc)
+
+ # Save to file and plot solution
+ vtkfile << (u, t)
+
+ # Update previous solution
+ u_n.assign(u)