Can one solve system of ODE with fenics?

1

Howdy

I was just interested in solving
x''=x-2sin(x)
x(0)=0
x'(0)=1

I broke it into a dynamical system with x_1=x,x_2=x'

x_1' = x_2
x_2'=x_1 - 2sin(t)
x_1(0)=0
x_2(0)=1

Solution in my mind is sin(t)

I tried to solve it, and it seems not even initial condition is holding…

from dolfin import *
# from dolfin_adjoint import *
import matplotlib 
matplotlib.use("TkAgg")
import matplotlib.pyplot as plt
import numpy as np
import math
import time

#going to solve the ode in two ways 
#x'' = x +f
#In system form x_1=x, x_2=x'
#so system
#x_2'=x_1 + f
#x_1' = x_2 




#establish where initial condition will be placed i.e 
#vector solution  [x1(shift),x2(shift)] = [x10,x20] 
def bd_func(shift):

	def boundary(x):
		return x[0] - shift < DOLFIN_EPS
	return boundary


#We solve the ode on [0,t_final]
t_final = 2*math.pi
#we have a system of ODE - two unknown functions
number_of_odes = 2
#divide the interval into 1000 elements
mesh = IntervalMesh(500, 0, t_final)
#for solving system of first order ODE
W = VectorFunctionSpace(mesh, 'P', 10, dim = number_of_odes)
Ws = FunctionSpace(mesh, 'P', 10)



bcsys = DirichletBC(W, [Constant(0.0),Constant(1.0)] ,bd_func(0.0))


x = Function(W)
v = TestFunction(W)
f = Expression('-2.0*sin(x[0])', degree=2)
#assign(x,interpolate(Expression('sin(x[0])', degree=2), ))


x1,x2 = split(x)
v1,v2 = split(v)

v1der = v1.dx(0)
v2der = v2.dx(0)

weak_form = -x2*v2der*dx - x1*v2*dx -f*v2*dx -x1*v1der*dx -x2*v1*dx
Jac = derivative(weak_form, x, TrialFunction(W))


solve(weak_form==0,x,J=Jac,bcs=bcsys)
x1,x2 = split(x)
plot(x1, title="state")
plt.savefig("state.pdf")
plt.clf()
plot(x2, title="derivative of state")
plt.savefig("derv.pdf")
plt.clf()
2

Hi, the following works

from dolfin import *


def main(ncells):
    mesh = IntervalMesh(ncells, 0, 2*pi)

    Welm = MixedElement([FiniteElement('Lagrange', interval, 2),
                         FiniteElement('Lagrange', interval, 1)])
    W = FunctionSpace(mesh, Welm)

    bcsys = [DirichletBC(W.sub(0), Constant(0.0), 'near(x[0], 0)'),
             DirichletBC(W.sub(1), Constant(1.0), 'near(x[0], 0)')]

    up = Function(W)
    u, p = split(up)
    v, q = split(TestFunction(W))
    source = Expression('2.0*sin(x[0])', degree=2)

    weak_form = u.dx(0)*v*dx - p*v*dx + p.dx(0)*q*dx - u*q*dx  + source*q*dx
    Jac = derivative(weak_form, up, TrialFunction(W))

    solve(weak_form == 0, up, J=Jac, bcs=bcsys)

    u, p = split(up)
        
    x, = SpatialCoordinate(mesh)

    print 'u error', sqrt(abs(assemble(inner(u - sin(x), u - sin(x))*dx)))
    print 'p error', sqrt(abs(assemble(inner(p - cos(x), p - cos(x))*dx)))


map(main, [10**i for i in range(2, 5)])
4 Likes