I am experimenting with using diffeqpy to solve a PDE by discretization of the spacial dimension, while I treat the time dimension as a set of ordinary differential equations. I managed to solve a very simple problem using a for loop. However, when I try to use matrixes, to scale the problem up, the solver provides incorrect answers.
The following piece of code works:
from diffeqpy import de
import numpy as np
def f(du,u,p,t):
#define shape of matrix
s = (6,7)
cc = np.matrix((np.zeros(s)))
for j in range(0,6):
for i in range(0,6):
if (j == i):
cc[j,i] = -1.0
cc[j,i+1] = 1.0
for j in range(0,6):
du[j] = cc[j,0]*u[0] + cc[j,1]*u[1] + cc[j,2]*u[2] + cc[j,3]*u[3] + cc[j,4]*u[4] + cc[j,5]*u[5] + cc[j,6]*u[6]
u0 = [0.1,0.0,0.0,0.0,0.0,0.0,1.0]
tspan = (0., 20.)
prob = de.ODEProblem(f, u0, tspan)
sol = de.solve(prob)
This codes is similar to the following piece of code that also works:
from diffeqpy import de
def f(du,u,p,t):
du[0] = -u[0]+u[1]
du[1] = -u[1]+u[2]
du[2] = -u[2]+u[3]
du[3] = -u[3]+u[4]
du[4] = -u[4]+u[5]
du[5] = -u[5]+u[6]
u0 = [0.1,0.0,0.0,0.0,0.0,0.0,1.0]
tspan = (0., 20.)
prob = de.ODEProblem(f, u0, tspan)
sol = de.solve(prob)
However, when I try and use a matrix operation, the problem just does not solve correctly. I don't have a background in computer science. However, I would like to learn more. Why is the following piece of code not working? Has it got to do with mutable vs. immutable object? How can I utilize a matrix to make this problem scale to bigger discretisation steps?
from diffeqpy import de
import numpy as np
def f(du,u,p,t):
#define shape of matrix
s = (6,7)
cc = np.matrix((np.zeros(s)))
for j in range(0,6):
for i in range(0,6):
if (j == i):
cc[j,i] = -1.0
cc[j,i+1] = 1.0
x = np.matrix(u).T
du = (cc*x).T
u0 = [0.1,0.0,0.0,0.0,0.0,0.0,1.0]
tspan = (0., 20.)
prob = de.ODEProblem(f, u0, tspan)
sol = de.solve(prob)
I would appreciate any guidance on this problem.
from diffeqpy import de
import numpy as np
def f(u,p,t):
#define shape of matrix
s = (6,6)
cc = np.matrix((np.zeros(s)))
for j in range(0,6):
for i in range(0,5):
if (j == i):
cc[j,i] = -1.0
cc[j,i+1] = 1.0
cc[-1,-1] = -1.0
x = (np.matrix(u).T)
du = (list(((cc*x).T).flat))
return du
u0 = [0.1,0.0,0.0,0.0,0.0,0.1]
tspan = (0., 20.)
prob = de.ODEProblem(f, u0, tspan)
sol = de.solve(prob)