# Euler's method for matrices

**URL:** <https://fenicsproject.discourse.group/t/eulers-method-for-matrices/4766>\
**Category:** variational formulation\
**Created:** [December 18, 2020, 4:04pm UTC](https://fenicsproject.discourse.group/t/eulers-method-for-matrices/4766 "2020-12-18T16:04:10Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Cou](https://avatars.discourse-cdn.com/v4/letter/c/919ad9/32.png) [@Cou](https://fenicsproject.discourse.group/u/Cou)\
**Post date:** [December 18, 2020, 4:04pm UTC](https://fenicsproject.discourse.group/t/eulers-method-for-matrices/4766/1 "2020-12-18T16:04:11Z")

</div>

Dear community,

In order to calculate a concentration at the borders of the domain, named c\_s, by the Euler’s method, I wrote the following code

```auto
# dc_s / dt = integral^{Gamma} (u - c_s) dA
def func_c_s(t, c_s, u_prev ):    
    return ((u_prev-c_s)*ds)
   
def euler_c_s( t0, c_s, dt, t, u_prev ):
    temp = -0

    # Iterating till the point at which we 
    # need approximation 
    while t0 < t: 
        temp = c_s
        c_s = c_s + dt * func_c_s(t0, c_s, u_prev) 
        t0 = t0 + dt
    return c_s
       
def value_shape(self):
    return ()

```

, which is run by

```auto
while t < T:
    t = t+dt
    c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)

```

,whereas `c_s_0` is the intial value for t=0 of c\_s and described by

```auto
class My_c_s_t0(UserExpression):
    def eval(self, value, x):
                
        # boundaries
        if x[0] == 0:
            value[0] = 0.2
        
        elif x[0] == 1:
            value[0] = 0.2

        elif x[1] == 0:
            value[0] = 0.2
        
        elif x[1] == 1:
            value[0] = 0.2
                
        # domain without boundaries
        else:
            value[0] = 0

def value_shape(self):
    return ()

```

```auto
c_s_t0 = My_c_s_t0()
c_s_0 = interpolate(c_s_t0, V)

```

and `u_prev`, which is usually determined by `solve(a == L, u)`, but here, in order to reduce the complexity of the MWCE given by

```auto
class My_u_t0(UserExpression):
    def eval(self, value, x):
                
        if x[0] >= 0 and x[0] <= 0.4:
            value[0] = u_avg_domain_t0

        elif x[0] >= 0.6 and x[0] <= 1.0:
            value[0] = u_avg_domain_t0
                
        else:
            value[0] = 0
        
def value_shape(self):
    return ()

```

```auto
u_t0 = My_u_t0()
u_prev = interpolate(u_t0, V)

```

When running the MWCE, I get the following error code

```auto
  File "20201217_c_s.py", line 150, in <module>
    c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)
  File "20201217_c_s.py", line 43, in euler_c_s
    c_s = c_s + dt * func_c_s(t0, c_s, u_prev) 
TypeError: unsupported operand type(s) for +: 'Function' and 'Form'

```

**I´m unsure, if the error message is a result of a wrong implementation of the Euler´s method** (Before, I only used Euler´s method for problems with scalar initial values) **or** `func_c_s` **can’t handle the matrix of** `c_s` **in its current form?**

Thanks a lot for your time,  
Cou

> **Full MWCE**
>
> ```auto
> from fenics import *
> import numpy as np
> 
> ############### define concentration c_s (on boundaries) for t=0 ##############
> 
> class My_c_s_t0(UserExpression):
> def eval(self, value, x):
>                 
> # boundaries
> if x[0] == 0:
> value[0] = 0.2
>         
> elif x[0] == 1:
> value[0] = 0.2
> 
> elif x[1] == 0:
> value[0] = 0.2
>         
> elif x[1] == 1:
> value[0] = 0.2
>                 
> # domain without boundaries
> else:
> value[0] = 0
> 
> def value_shape(self):
> return ()
> 
> ############################## Euler method ###################################
> # dc_s / dt = integral^{Gamma} (u - c_s) dA
> def func_c_s(t, c_s, u_prev ):    
> return ((u_prev-c_s)*ds)
>    
> def euler_c_s( t0, c_s, dt, t, u_prev ):
> temp = -0
> 
> # Iterating till the point at which we 
> # need approximation 
> while t0 < t: 
> temp = c_s
> c_s = c_s + dt * func_c_s(t0, c_s, u_prev) 
> t0 = t0 + dt
> return c_s
>        
> def value_shape(self):
> return ()  
> 
> ######################## define u for t=0 #####################################
> 
> class My_u_t0(UserExpression):
> def eval(self, value, x):
>                 
> if x[0] >= 0 and x[0] <= 0.4:
> value[0] = u_avg_domain_t0
> 
> elif x[0] >= 0.6 and x[0] <= 1.0:
> value[0] = u_avg_domain_t0
>                 
> else:
> value[0] = 0
>         
> def value_shape(self):
> return ()
> 
> ######################## define mesh and function space #######################
> 
> nx = ny = 16
> mesh = UnitSquareMesh(nx, ny)
> space_dim = mesh.geometry().dim()
> V = FunctionSpace(mesh, 'CG', 1)
> 
> ################################ geometry #####################################
> 
> class Left(SubDomain):
> def inside(self, x, on_boundary):
> return near(x[0], 0.0)
> 
> class Right(SubDomain):
> def inside(self, x, on_boundary):
> return near(x[0], 1.0)
> 
> class Bottom(SubDomain):
> def inside(self, x, on_boundary):
> return near(x[1], 0.0)
> 
> class Top(SubDomain):
> def inside(self, x, on_boundary):
> return near(x[1], 1.0)
> 
> class Obstacle_left(SubDomain):
> def inside(self, x, on_boundary):
> return (between(x[0], (0.0, 0.4)))
>     
> class Obstacle_right(SubDomain):
> def inside(self, x, on_boundary):
> return (between(x[0], (0.6, 1.0)))
>     
> # Initialize sub-domain instances
> left = Left()
> top = Top()
> right = Right()
> bottom = Bottom()
> obstacle_left = Obstacle_left()
> obstacle_right = Obstacle_right()
> 
> # Initialize mesh function for interior domains
> domains = MeshFunction("size_t", mesh, space_dim-1)
> domains.set_all(0)
> obstacle_left.mark(domains, 1)
> obstacle_right.mark(domains, 2)
> 
> # Initialize mesh function for boundary domains
> boundaries = MeshFunction("size_t", mesh, space_dim-1)
> boundaries.set_all(0)
> left.mark(boundaries, 1)
> top.mark(boundaries, 2)
> right.mark(boundaries, 3)
> bottom.mark(boundaries, 4)
> 
> # Define new measures associated with the interior domains and
> # exterior boundaries
> dx = Measure('dx', domain=mesh, subdomain_data=domains)
> ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
> 
> ##################### define constants #########################################
> t = 0
> t0 = 0
> T = 2.0                     
> num_steps = 10              
> dt = T / num_steps          
> 
> u_avg_domain_t0 = 1   
> u_t0 = My_u_t0()
> c_s_t0 = My_c_s_t0()
> 
> u_prev = interpolate(u_t0, V)
> c_s_0 = interpolate(c_s_t0, V)
> 
> c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)
> print('t = %.2f, c_s(0,0) = %.3f' % (t, c_s(0,0)))
> 
> while t < T:
> t = t+dt
> c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)
> print('t = %.2f, c_s(0,0) = %.3f' % (t, c_s(0,0)))
> 
> ```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [December 21, 2020, 8:03am UTC](https://fenicsproject.discourse.group/t/eulers-method-for-matrices/4766/2 "2020-12-21T08:03:30Z")

</div>

As you have not supplied a minimal working code example (one that can be executed and reproduce your error by copy-pasting it into a text editor and running it), it slightly tricky to debug.  
My best guess is that if you remove the outer brackets on the following term:

> [@Cou](#):
>
> `((u_prev-c_s)*ds)`

your code will work.

---

<div class="post-metadata">

**Author:** ![Cou](https://avatars.discourse-cdn.com/v4/letter/c/919ad9/32.png) [@Cou](https://fenicsproject.discourse.group/u/Cou)\
**Post date:** [January 6, 2021, 1:31pm UTC](https://fenicsproject.discourse.group/t/eulers-method-for-matrices/4766/3 "2021-01-06T13:31:59Z")

</div>

Thanks @dokken and happy New Year to everyone!

I removed the outer brackets, but unfortunately I´m faced with the same error message.

```auto
  File "20210106_c_s.py", line 150, in <module>
    c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)
  File "20210106_c_s.py", line 43, in euler_c_s
    c_s = c_s + dt * func_c_s(t0, c_s, u_prev) 
TypeError: unsupported operand type(s) for +: 'Function' and 'Form'

```

That´s my MWCE:

```auto
from fenics import *
import numpy as np

############### define concentration c_s (on boundaries) for t=0 ##############

class My_c_s_t0(UserExpression):
    def eval(self, value, x):
                
        # boundaries
        if x[0] == 0:
            value[0] = 0.2
        
        elif x[0] == 1:
            value[0] = 0.2

        elif x[1] == 0:
            value[0] = 0.2
        
        elif x[1] == 1:
            value[0] = 0.2
                
        # domain without boundaries
        else:
            value[0] = 0

def value_shape(self):
    return ()

############################## Euler method ###################################
# dc_s / dt = integral^{Gamma} (u - c_s) dA
def func_c_s(t, c_s, u_prev ):    
    return (u_prev-c_s)*ds
   
def euler_c_s( t0, c_s, dt, t, u_prev ):
    temp = -0

    # Iterating till the point at which we 
    # need approximation 
    while t0 < t: 
        temp = c_s
        c_s = c_s + dt * func_c_s(t0, c_s, u_prev) 
        t0 = t0 + dt
    return c_s
       
def value_shape(self):
    return ()  

######################## define u for t=0 #####################################

class My_u_t0(UserExpression):
    def eval(self, value, x):
                
        if x[0] >= 0 and x[0] <= 0.4:
            value[0] = u_avg_domain_t0

        elif x[0] >= 0.6 and x[0] <= 1.0:
            value[0] = u_avg_domain_t0
                
        else:
            value[0] = 0
        
def value_shape(self):
    return ()

######################## define mesh and function space #######################

nx = ny = 16
mesh = UnitSquareMesh(nx, ny)
space_dim = mesh.geometry().dim()
V = FunctionSpace(mesh, 'CG', 1)

################################ geometry #####################################

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)

class Obstacle_left(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[0], (0.0, 0.4)))
    
class Obstacle_right(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[0], (0.6, 1.0)))
    
# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
obstacle_left = Obstacle_left()
obstacle_right = Obstacle_right()

# Initialize mesh function for interior domains
domains = MeshFunction("size_t", mesh, space_dim-1)
domains.set_all(0)
obstacle_left.mark(domains, 1)
obstacle_right.mark(domains, 2)

# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, space_dim-1)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

##################### define constants #########################################
t = 0
t0 = 0
T = 2.0                     
num_steps = 10              
dt = T / num_steps          

u_avg_domain_t0 = 1   
u_t0 = My_u_t0()
c_s_t0 = My_c_s_t0()

u_prev = interpolate(u_t0, V)
c_s_0 = interpolate(c_s_t0, V)

c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)
print('t = %.2f, c_s(0,0) = %.3f' % (t, c_s(0,0)))

while t < T:
    t = t+dt
    c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)
    print('t = %.2f, c_s(0,0) = %.3f' % (t, c_s(0,0)))

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 6, 2021, 1:39pm UTC](https://fenicsproject.discourse.group/t/eulers-method-for-matrices/4766/4 "2021-01-06T13:39:02Z")

</div>

> [@Cou](#):
>
> ```auto
> c_s = c_s + dt * func_c_s(t0, c_s, u_prev) 
> TypeError: unsupported operand type(s) for +: 'Function' and 'Form'
> 
> ```

The initial `c_s` is a function, with no integration domain, thus the error.  
I am not sure what you are trying to achieve in this loop.

---

<div class="post-metadata">

**Author:** ![Cou](https://avatars.discourse-cdn.com/v4/letter/c/919ad9/32.png) [@Cou](https://fenicsproject.discourse.group/u/Cou)\
**Post date:** [January 21, 2021, 3:34pm UTC](https://fenicsproject.discourse.group/t/eulers-method-for-matrices/4766/5 "2021-01-21T15:34:02Z")

</div>

> [@dokken](#):
>
> > [@Cou](#):
> >
> > ```auto
> > > c_s = c_s + dt * func_c_s(t0, c_s, u_prev) 
> > [...]
> > 
> > ```
> 
> […]  
> I am not sure what you are trying to achieve in this loop.

Thanks for pointing that out, because after some thinking I´m unsure, if my approach is logical sound.

Anyway, my goal is to calculate c\_s, which is the concentration outside of the domain, described by  
\frac{c\_s(t)}{dt}=\int\_{\Gamma}^{}(u(\mathbf{x},t)-c\_s(t))dA  
and by implementing the Euler´s method  
c\_{s}^{n+1}=c\_s^{n}+\Delta t\int\_{\Gamma}^{}(u(\mathbf{x},t)-c\_s(t))dA  
like shown [here](https://www.geeksforgeeks.org/euler-method-solving-differential-equation/).

Because c\_s is a scalar value by definition, but I have to substract c\_s from u, I decided to “scale” c\_s “up” by using (`c_s_t0 = df.zero((nx,ny))` was not working)

```auto
class My_c_s_t0(UserExpression):
    def eval(self, value, x):
                
        # boundaries
        if x[0] == 0:
            value[0] = 0
        
        elif x[0] == 1:
            value[0] = 0

        elif x[1] == 0:
            value[0] = 0
        
        elif x[1] == 1:
            value[0] = 0
                
        # domain without boundaries
        else:
            value[0] = 0

def value_shape(self):
    return ()

c_s_t0 = My_c_s_t0()
c_s_0 = interpolate(c_s_t0, V)

```

in order to obtain the same number of rows and columns as u.  
By the way, this is an update to the previous MWCE, because c\_s(t=0)=0.

So, besides the eventual lack of logical soundness, how am I able to add an integration domain to c\_s?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 22, 2021, 12:06am UTC](https://fenicsproject.discourse.group/t/eulers-method-for-matrices/4766/6 "2021-01-22T00:06:49Z")

</div>

> [@Cou](#):
>
> `dt * func_c_s(t0, c_s, u_prev)`

You would Need to assemble this (so that it becomes a scalar and can be added to a float). However, in your code, `c_s` is a function (that is spatially varying, as it is defined by a function. I guess you could use the “R” space for “c”: see for instance: [https://bitbucket.org/fenics-project/dolfin/src/b495043d6b3914383eb939ea6c6794442080a3a5/python/demo/documented/neumann-poisson/demo\_neumann-poisson.py.rst?at=master#demo\_neumann-poisson.py.rst-48,71,134:135](https://bitbucket.org/fenics-project/dolfin/src/b495043d6b3914383eb939ea6c6794442080a3a5/python/demo/documented/neumann-poisson/demo_neumann-poisson.py.rst?at=master#demo_neumann-poisson.py.rst-48,71,134:135)