0
$\begingroup$

I want to plot the 100,200 and 400 iterations of this function of non homogeneous parabolic pde

def solver(L, Nx, Nt, T, theta=0.5):
    

    x = np.linspace(-L,L, Nx+1)   # mesh points in space
    dx = (L-(-L))/(Nx-1)
    dt = 1/Nt

    Nt = int(round(T/float(dt)))
    t = np.linspace(0, T, Nt+1)   # mesh points in time

    h = (L-(-L))/(Nx+1)
    t = T / Nt
    m = t/h**2
    #print("m =", round(m,2))

    u   = np.zeros(Nx+1)   # solution array at t[n+1]
    u_n = np.zeros(Nx+1)   # solution at t[n]
    a   = np.zeros(Nx+1) + (1+x**2)
    c   = np.ones(Nx+1) * -1

    Dl = 0.5*theta
    Dr = 0.5*(1-theta)

   
    diagonal = np.zeros(Nx+1)
    lower    = np.zeros(Nx)
    upper    = np.zeros(Nx)
    b        = np.zeros(Nx+1)

 
    diagonal[1:-1] = 1 + Dl*(a[2:] + 2*a[1:-1] + a[:-2])
    lower[:-1] = -Dl*(a[1:-1] + a[:-2])
    upper[1:]  = -Dl*(a[2:] + a[1:-1])
    
    
    diagonal[0] = 1
    upper[0] = 0
    diagonal[Nx] = 1
    lower[-1] = 0

    A = scipy.sparse.diags(
        diagonals=[diagonal, lower, upper],
        offsets=[0, -1, 1],
        shape=(Nx+1, Nx+1),
        format='csr')
    
    def I(x):
        return(np.sin(np.pi*x))

   
    for i in range(0,Nx+1):
        u_n[i] = I(x[i])

   

    # Time loop
    for n in range(0, Nt):
        b[1:-1] = u_n[1:-1] + Dr*(
            (a[2:] + a[1:-1])*(u_n[2:] - u_n[1:-1]) -
            (a[1:-1] + a[0:-2])*(u_n[1:-1] - u_n[:-2])) + dt*theta*c[1:-1] + dt*(1-theta)*c[1:-1]
        
        # Boundary conditions
        b[0]  = 0
        b[-1] = 0
        # Solve
        u[:] = scipy.sparse.linalg.spsolve(A, b)

        u_n,u = u, u_n 
    
        plt.plot(u[:,100])

I receive this error:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-91-41c85be859a3> in <module>
----> 1 Back = solver(L=1, Nx=39, Nt=400, T=1, theta=1)

<ipython-input-90-ef7886a9d8f0> in solver(L, Nx, Nt, T, theta)
     68         u_n,u = u, u_n
     69 
---> 70         plt.plot(u[:,100])
     71 

IndexError: too many indices for array


Can someone help me?

$\endgroup$
10
  • $\begingroup$ U a vector. Just use plt.plot(u) I guess. $\endgroup$ – Laurent90 Dec 1 '20 at 18:18
  • $\begingroup$ it plots all the iterations if do this $\endgroup$ – user37062 Dec 1 '20 at 18:19
  • $\begingroup$ No, as far as I understand u is a vector that has the same dimension as your number of mesh points. U is then the value of your field at time t. $\endgroup$ – Laurent90 Dec 1 '20 at 19:13
  • $\begingroup$ The same I thought.But You can run it.You will see. $\endgroup$ – user37062 Dec 1 '20 at 20:08
  • 1
    $\begingroup$ For the question at hand: Laurent90 is of course right -- u is a vector, and u[:,100] correctly points out that you can't subscript it with two indices. The syntax you use would extract the 100'th column from a matrix, but it is not the right syntax for extracting a subset of a vector. $\endgroup$ – Wolfgang Bangerth Dec 1 '20 at 20:17
2
$\begingroup$

Alright I've taken a closer look a the code. One first advice that I could give you is to try to understand how the code works. It is very easy to understand that the for-loop is a time loop, i.e. the system is advanced forward one time step at each iteration of the loop. The time marching equation uses a theta-scheme for the time discretization (theta=0.5 --> Crank-Nicolson scheme), and the resulting linear system is solved with the sparse linear solver from Scipy. The resolution of this linear system gives the vector $u$ of discrete temperatures (or the actual physical value whose evolution you're computing) on every mesh point. The way you code is written, it store the past values of u. SO either you save (say, in a list for example) or you just use the value of $u$ on the fly when you want to.

By the way, the method to update $Nt$ was wrong, as well as the method to compute $dt$, so I've modified your code accordingly.

Also, try to post a code that is directly usable ;) I had to import some modules, and also setup the call to the solver function to make it work.

Here is the final result with the call solver(L=1, Nx=30, Nt=400, T=200., theta=1) (so 30 elements):

enter image description here

By the way, something seems wrong with the model, as it does not converge in space. For example, with 3000 points instead of 30 I get the next figure. And the situation only worsens as the number of points $N_x$ is increased... I don't know what is is you are modelling, but that should be checked. I also see that the mesh spacing dx is not used in your code, which seems rather surprising.

enter image description here

And the code is:

import numpy as np
import matplotlib.pyplot as plt
import scipy.sparse.linalg

def solver(L, Nx, Nt, T, theta=0.5):
    x = np.linspace(-L,L, Nx+1)   # mesh points in space
    dx = (L-(-L))/(Nx-1)
       
    t = np.linspace(0, T, Nt+1)   # mesh points in time
    dt = t[1]-t[0]

    h = (L-(-L))/(Nx+1)

    u   = np.zeros(Nx+1)   # solution array at t[n+1]
    u_n = np.zeros(Nx+1)   # solution at t[n]
    a   = np.zeros(Nx+1) + (1+x**2)
    c   = np.ones(Nx+1) * -1

    Dl = 0.5*theta
    Dr = 0.5*(1-theta)
   
    diagonal = np.zeros(Nx+1)
    lower    = np.zeros(Nx)
    upper    = np.zeros(Nx)
    b        = np.zeros(Nx+1)
 
    diagonal[1:-1] = 1 + Dl*(a[2:] + 2*a[1:-1] + a[:-2])
    lower[:-1] = -Dl*(a[1:-1] + a[:-2])
    upper[1:]  = -Dl*(a[2:] + a[1:-1])
    
    diagonal[0] = 1
    upper[0] = 0
    diagonal[Nx] = 1
    lower[-1] = 0

    A = scipy.sparse.diags(
        diagonals=[diagonal, lower, upper],
        offsets=[0, -1, 1],
        shape=(Nx+1, Nx+1),
        format='csr')
    
    def I(x):
        return(np.sin(np.pi*x))

    u_n[:Nx+1] = I(x[:Nx+1])
      
    # Time loop
    for n in range(0, Nt):
        print('iter {}/{}'.format(n,Nt))
        b[1:-1] = u_n[1:-1] + Dr*(
            (a[2:] + a[1:-1])*(u_n[2:] - u_n[1:-1]) -
            (a[1:-1] + a[0:-2])*(u_n[1:-1] - u_n[:-2])) + dt*theta*c[1:-1] + dt*(1-theta)*c[1:-1]
        
        # Boundary conditions
        b[0]  = 0
        b[-1] = 0
        # Solve
        u[:] = scipy.sparse.linalg.spsolve(A, b)

        u_n,u = u, u_n 
        
        if 0 :# One plot for every selected iteration
          if np.mod(n,100)==0: # every 100 iterations
            plt.figure()
            plt.plot(x,u)
            plt.grid()
            plt.xlabel('x')
            plt.ylabel('u')
            plt.title('Solution at iteration {}'.format(n))
            plt.show()
        else: # one common figure for all the plotted iterations
          if np.mod(n,75)==0: # every 100 iterations
            if n==0: # instantiate the figure
              plt.figure()
              plt.grid()
              plt.xlabel('x')
              plt.ylabel('u')
            # plot the current iteration
            plt.plot(x,u, label='t={:.2f}'.format(t[n]))
            
            
if __name__=='__main__':
  solver(L=1, Nx=3000, Nt=400, T=200., theta=1)
  plt.legend(loc='center left', ncol=2)

EDIT: from looking at the original file you said you adapted (https://github.com/hplgit/fdm-book/blob/master/src/diffu/diffu1D_vc.py), t seems the definition of Dl and Dr does not make sense with regards to the discretisation of the diffusion operator. It should involve the diffusion coefficient and the mesh spacing, not the coefficient theta of the temporal scheme.

$\endgroup$
7
  • $\begingroup$ The problem that I try to solve is this:Consider the following inhomogeneous parabolic initial/boundary value problem: $$u_{t}(t,x) = (1+x^{2})u_{xx}(t,x)+u(t,x),$$ for $t \in [0,1]$ and $x \in [-1,1]$ $$u(0,x) = \sin(\pi x),$$ for $x \in [-1,1]$ initial condition $$u(t,-1)=u(t,1)=0,$$ $t \in [0,1]$ Dirichlet boundary conditions. Construct a Backward Euler method with $N_{x} =39$ and $Nt=400$ (i.e$\mu=1$) and show the iterations $n=100,200,400.$ In addition construct a Crank Nicolson method with $N_{x} =39$ and $Nt=20$ (i.e$\mu=1/h$) and show the iterations $n=5,10,20.$ $\endgroup$ – user37062 Dec 1 '20 at 22:28
  • $\begingroup$ This is fairly easy to do starting from the code you have. You need to compute the terms of the Jacobian A for your new system. If you don't see how to do that, start by trying to discretise the classical diffusion equation $u_t = D u_{xx}$, then add a variable diffusion ($D(u)$ or $D(x)$). Gradually increase the complexity of your model and make sure you understand each step, otherwise you might end up not understanding and not learning much from the exercise. $\endgroup$ – Laurent90 Dec 1 '20 at 22:47
  • $\begingroup$ The code above does that.The non constant coefficient $\alpha$ is being added in A and also the constant $c= -1$. $\endgroup$ – user37062 Dec 1 '20 at 22:52
  • $\begingroup$ Then your problem is solved, right ? $\endgroup$ – Laurent90 Dec 1 '20 at 23:03
  • $\begingroup$ for me yes.Also that you said " Gradually increase the complexity of your model " is also succedd $\endgroup$ – user37062 Dec 1 '20 at 23:07

Your Answer

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