I am trying to solve :

$$-u''(x) + u(x) = \sin(2\pi x)\, ,\quad 0<x<1\, ,$$ $t>0$, with $u(0) = u(1) = 0$. That has as exact solution

$$u(x) = \frac{\sin(2\pi x)}{1 + 4\pi^2}\, .$$

But the Forward Euler approximation solution does not match the exact solution.

Any help?

import numpy as np
import matplotlib.pyplot as plt
L   = 1
Nx  = 19
Nt  = 800 
T   = 0.1
x   = np.linspace(0, L, Nx+1)    # mesh points in space
dx  = x[1] - x[0]
t   = np.linspace(0, T, Nt+1)    # mesh points in time
dt  = t[1] - t[0]
a   = 1
F   = a*dt/dx**2
u   = np.zeros(Nx+1)
u_n = np.zeros(Nx+1)

def I(x):
    return(np.sin( 2*x*np.pi ))

# Set initial condition u(x,0) = I(x)  
for i in range(0, Nx+1):
    u_n[i] = I(x[i])  

for n in range(0, Nt):
    # Compute u at inner mesh points
    for i in range(1, Nx):
        u[i] = u_n[i] + F*(u_n[i-1] - 2*u_n[i] + u_n[i+1])

    # Insert boundary conditions
    u[0]  = 0  
    u[Nx] = 0

    # Update u_n before next step
    u_n[:]= u
exact = np.sin(2*np.pi *x ) / (1+4*np.pi**2)
  • 1
    $\begingroup$ The code refers to initial conditions but the PDE is not time-dependent. Why? $\endgroup$
    – knl
    Sep 20 '20 at 19:38
  • $\begingroup$ I typed your equations, please check that they are OK. Also, please use MathJax for future posts. $\endgroup$
    – nicoguaro
    Sep 20 '20 at 23:02
  • 1
    $\begingroup$ You mention that your FE solution does not match the analytic solution. You could describe how they are different or show a plot. $\endgroup$
    – nicoguaro
    Sep 20 '20 at 23:03
  • $\begingroup$ i want to to take different values of like t =0,1,10,100 and see how it the temperature distribution diffuses on these values. $\endgroup$
    – user37062
    Sep 20 '20 at 23:52
  • 2
    $\begingroup$ At the least then, the question is poorly explained. You have a code that references a time variable, your description shows something else, and you don't tell us how the solution doesn't match your expectations. How about reducing the code to its minimal form that only solves the problem as described, and show us a picture with what it yields? $\endgroup$ Sep 21 '20 at 3:56

Comments above are right: it seems that you are also integrating in time (and indeed you also set the numer of points in time in your code), but the equation is only in variable $x$.

The following snippet produce the correct solution to your problem with linear elements in Python.

To compute $\int_0^1 \phi_i(x) f(x)dx$ I used integrate.quad from scipy, which performs Gaussian quadrature. That integral has been computed splitting the interval in $[x_{i-1},x_{i}]$ and $[x_{i},x_{i+1}]$, since the hat functions are non-differentiable at $x=x_i$. Btw, notice that this integral may be computed analytically, but quadrature is usually the choice in a fem solver.

    import numpy as np
    import matplotlib.pyplot as plt
    from scipy import integrate
    M = 10 #points in x
    L = 1 #endpoint
    x = np.linspace(0,L,M+1)
    h = x[1]-x[0]
    def uex(x):
        return(np.sin( 2*x*np.pi ))/(1+4*np.pi**2)
    def stiffassembly(M):
        diag = np.zeros(M-1) #x_1,...,x_{M-1} (M-1)
        offd = np.zeros(M-2) #off diagonal terms
        for i in range(1,M):
            diag[i-1] = 1/h +1/h
        for k in range(1,M-1):
            offd[k-1] = -1/h
        A = np.diag(offd,-1) + np.diag(diag,0) + np.diag(offd,+1)
        return A
    def massmatrix(N):
        diag = np.zeros(N-1) #x_1,...,x_M-1 (M-1)
        subd = np.zeros(N-2) 
        supr = np.zeros(N-2)
        for i in range(1,N):
            diag[i-1] = (h + h)/3
        for k in range(1,N-1):
            supr[k-1] = h/6
            subd[k-1] = h/6
        M = np.diag(subd,-1) + np.diag(diag,0) + np.diag(supr,+1)
        return M
    def f(x):
        return np.sin(2*np.pi*x)
    def load(M):
        load = np.zeros(M-1)
        for k in range(1,M):
            load[k-1] = integrate.quad(lambda w: f(w)*(1/h)*(w-x[k-1]),x[k-1],x[k])[0] + integrate.quad(lambda w: f(w)*(1/h)*(x[k+1]-w),x[k],x[k+1])[0]
        return load
    u = np.linalg.solve(+stiffassembly(M) + massmatrix(M),load(M))
    U = np.r_[0,u,0]

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy