@mehdi jafari 

Hi Mehdi,

sorry, I was tried last days.

@sarra 

I like your method, nice method. For me, I write three differents method without using maple, and i get the same answer...nice good. I can send a theoretical answer ( without maple) if you need.

Thank you for helping me in Maple.

@Preben Alsholm 

Hi Alsholm.

Thank you for the answer.

@Carl Love 

Thank you for remark.

@Carl Love 

Hi Carl,

Many thanks for your remark.  It's very nice. I have another, question, how can we collect phi(n*h) in your last line and put the resulat in only one sum.

second question, if I have three  terms in the sum,  can I use this:  I tried it, but doesn't working

subs([m= n, Int= sum], applyop(Change, 2, %, m= n+1),applyop(Change, 3, %, m= n+2));

@Markiyan Hirnyk 

Thank you for your answer. I tried the method in yoyr link.

restart;Test:=Sum(f(n)+f(n+1),n=1..N);
                             N                    
                           -----                  
                            \                     
                             )                    
                            /    (f(n) + f(n + 1))
                           -----                  
                           n = 1                  
map(op(0,Test),op(1,Test));
                       /-----     \   /-----         \
                       | \        |   | \            |
                       |  )       |   |  )           |
                       | /    f(n)| + | /    f(n + 1)|
                       \-----     /   \-----         /
with(student):
eq1:= expand(Test);

                       /  N       \   /  N           \
                       |-----     |   |-----         |
                       | \        |   | \            |
                       |  )       |   |  )           |
                       | /    f(n)| + | /    f(n + 1)|
                       |-----     |   |-----         |
                       \n = 1     /   \n = 1         /
collect(eq1, f(n));
                       /  N       \   /  N           \
                       |-----     |   |-----         |
                       | \        |   | \            |
                       |  )       |   |  )           |
                       | /    f(n)| + | /    f(n + 1)|
                       |-----     |   |-----         |
                       \n = 1     /   \n = 1         /

But, for me I want to rewrite the sum with using only f(n) without f(n+1)

@geischtli 

There is no error in the code, a get an animate curve.

restart;
with(PDEtools):
with(plots):
with(LinearAlgebra):
Digits := 14:


The purpose of this worksheet is to discuss the solution of nonlinear PDEs using finite difference methods.  For concreteness, we will concentrate on a diffusion equation with a potential:

     "(∂)/(∂ t)u(t,x)=d((∂)^(2))/(∂

        x^(2))u(t,x)+V(u(t,x))."

Here,
                                      d
 is the diffusion constant and
                                      V
 is a given (presumably nonlinear) function of
                                   u(t, x)
.  We take t as a time variable, hence this is an initial value problem.  We'll take Dirichlet boundary conditions

                             "u(t,-L)=u(t,L)=0,"

where
                                      L
 is a constant to be specified.
Stencils
Here is our PDE:
pde := diff(u(t,x),t)=d*diff(u(t,x),x,x)+V(u(t,x));
               d              / d  / d         \\             
              --- u(t, x) = d |--- |--- u(t, x)|| + V(u(t, x))
               dt             \ dx \ dx        //             
We will make use of the GenerateStencil procedure defined elsewhere to generate various substencils for the derivatives appearing in (1.1):
GenerateStencil := proc(F,N,{orientation:=center,stepsize:=h,showorder:=true,showerror:=false})
    local vars, f, ii, Degree, stencil, Error, unknowns, Indets, ans, Phi, r, n, phi;

    Phi := convert(F,D);
    vars := op(Phi);
    n := PDEtools[difforder](Phi);
    f := op(1,op(0,Phi));  
    if (nops([vars])<>1) then:           
      r := op(1,op(0,op(0,Phi)));
    else:
      r := 1;      
    fi:
    phi := f(vars);
    if (orientation=center) then:
      if (type(N,odd)) then:
        ii := [seq(i,i=-(N-1)/2..(N-1)/2)];
      else:
        ii := [seq(i,i=-(N-1)..(N-1),2)];    
      fi;
    elif (orientation=left) then:
      ii := [seq(i,i=-N+1..0)]:
    elif (orientation=right) then:
      ii := [seq(i,i=0..N-1)]:
    fi;
    stencil := add(a[ii[i]]*subsop(r=op(r,phi)+ii[i]*stepsize,phi),i=1..N);
    Error := D[r$n](f)(vars) - stencil;
    Error := convert(series(Error,stepsize,N),polynom);
    unknowns := {seq(a[ii[i]],i=1..N)};
    Indets := indets(Error) minus {vars} minus unknowns minus {stepsize};
    Error := collect(Error,Indets,'distributed');
    ans := solve({coeffs(Error,Indets)},unknowns);
    if (ans=NULL) then:
       print(`Failure: try increasing the number of points in the stencil`);
       return NULL;
    fi:
    stencil := subs(ans,stencil);
    Error := convert(series(`leadterm`(D[r$n](f)(vars) - stencil),stepsize,N+20),polynom);
    Degree := degree(Error,stepsize);
    if (showorder) then:
       print(cat(`This stencil is of order `,Degree));
    fi:
    if (showerror) then:
       print(cat(`This leading order term in the error is `,Error));
    fi:
    convert(D[r$n](f)(vars) = stencil,diff);

end proc:

forward_time := GenerateStencil(diff(u(t,x),t),2,orientation=right,stepsize=s);
backward_time := GenerateStencil(diff(u(t,x),t),2,orientation=left,stepsize=s);
centered_space := GenerateStencil(diff(u(t,x),x,x),3,orientation=center,stepsize=h);
                         This stencil is of order 1
                     d              u(t, x)   u(t + s, x)
                    --- u(t, x) = - ------- + -----------
                     dt                s           s     
                         This stencil is of order 1
                     d              u(t - s, x)   u(t, x)
                    --- u(t, x) = - ----------- + -------
                     dt                  s           s   
                         This stencil is of order 2
           d  / d         \   u(t, x - h)   2 u(t, x)   u(t, x + h)
          --- |--- u(t, x)| = ----------- - --------- + -----------
           dx \ dx        /        2            2            2     
                                  h            h            h      
The FTCS stencil is obtained by putting forward_time and centered_space into the PDE:
Label[1] := `FTCS`:
stencil[1] := subs(forward_time,centered_space,pde);
       u(t, x)   u(t + s, x)     /u(t, x - h)   2 u(t, x)   u(t, x + h)\
     - ------- + ----------- = d |----------- - --------- + -----------|
          s           s          |     2            2            2     |
                                 \    h            h            h      /

        + V(u(t, x))
The BTCS stencil is obtained by putting backward_time and centered_space into the PDE:
BTCS := subs(backward_time,centered_space,t=t+s,pde);
  u(t, x)   u(t + s, x)     /u(t + s, x - h)   2 u(t + s, x)   u(t + s, x + h)
- ------- + ----------- = d |--------------- - ------------- + ---------------
     s           s          |       2                2                2       
                            \      h                h                h        

  \                 
  | + V(u(t + s, x))
  |                 
  /                 
And the Crank-Nicholson stencil comes from averaging the other two:
Label[2] := `CN`:
stencil[2] := (stencil[1]+BTCS)/2;
   u(t, x)   u(t + s, x)   1   /u(t, x - h)   2 u(t, x)   u(t, x + h)\
 - ------- + ----------- = - d |----------- - --------- + -----------|
      s           s        2   |     2            2            2     |
                               \    h            h            h      /

      1              1   /u(t + s, x - h)   2 u(t + s, x)   u(t + s, x + h)\
    + - V(u(t, x)) + - d |--------------- - ------------- + ---------------|
      2              2   |       2                2                2       |
                         \      h                h                h        /

      1               
    + - V(u(t + s, x))
      2               
Forward-time centered-space solution
The FTCS stencil can be solved for
                                 u(t + s, h)
 for any form of the potential
                                 V(u(t, x))
FTCS_stencil := expand(isolate(stencil[1],u(t+s,x)));
              s d u(t, x - h)   2 s d u(t, x)   s d u(t, x + h)               
u(t + s, x) = --------------- - ------------- + --------------- + s V(u(t, x))
                     2                2                2                      
                    h                h                h                       

   + u(t, x)
In other words, we can always calculation
                                      u
 at a future timestep with knowledge of the field at a previous timestep, irregardless of the nature of the nonlinearity in the potential.  This means that the numerical FTCS algorithm to solve ?? is virtually the same as it is for the linear diffusion equation.  Here is the code:
Subs := u(t,x-h)=u1,u(t,x)=u2,u(t,x+h)=u3:
Phi := unapply(subs(Subs,rhs(FTCS_stencil)),u1,u2,u3,h,s,d):

FTCS := proc(tau,L,N,M,d,f) local s, h, XX, X, T, PlotOptions, Title, u_past, u_future, p, i, j:

        X := j -> L*(-1 + j*2/(M+1));
        T := i -> tau*i/N;
        s := evalf(T(1)-T(0));
        h := evalf(X(1)-X(0));
        XX := Array(0..M+1,[seq(X(j),j=0..M+1)],datatype=float):
        u_past := Array(0..M+1,[seq(f(XX[j]),j=0..M+1)],datatype=float):
        p[0] := Frame[1](XX,u_past,s,h,T(0));
        for i from 1 to N do:
            u_future := Array(0..M+1,datatype=float):
            u_future[0] := 0:
            u_future[M+1] := 0:
            for j from 1 to M do:
                 u_future[j] := Phi(u_past[j-1],u_past[j],u_past[j+1],h,s,d):
            od:
            ArrayTools[Copy](u_future,u_past):
            p[i] := Frame[1](XX,u_past,s,h,T(i));
        od:
        display(convert(p,list),insequence=true):

end proc:

Frame[1] := proc(xdata,ydata,s,h,T) local Title, PlotOptions:

        Title := typeset(`timestep = `,evalf[2](s),`, spacestep = `,evalf[2](h),`, `,t=evalf[2](T));
        PlotOptions := axes=boxed, labels=[x,u(t,x)], legend=["numerical solution (FTCS)"], color=red;
        plot(Matrix([[xdata],[ydata]])^%T,title=Title,PlotOptions);

end proc:
Here is some example output for a particular choice of potential:
V := u -> 1-u^2;
tau := 1;
N := 400;
M := 20;
L := 1;
d := 1;
s := evalf(2*L/(M+1));
h := evalf(tau/N);
f := x -> exp(-(4*x)^2);
movie[1] := FTCS(tau,L,N,M,d,f):
pds := pdsolve(pde,[u(t,-L)=0,u(t,+L)=0,u(0,x)=f(x)],numeric,timestep=s,spacestep=h):
movie[2] := pds:-animate(t=0..tau,frames=N+1,axes=boxed,legend=["pdsolve/numeric"],color=blue):
display([movie[1],movie[2]]);
          2
u -> 1 - u
                                      1
                                     400
                                     20
                                      1
                                      1
                                0.09523809524
                               0.002500000000
        /     2\
x -> exp\-16 x /

Our results compare favourably with those of pdsolve/numeric.

You can  test this code, there is an animation in time. 

animate_curve.mw

@Carl Love 

Hi.

First, there is no doubt that your remak helped me. Thank you very much.

Second, I modify the code, all the element of the code working fine, but i's seem for me that there is something to modify.

For example, the definition of the coeficients A[n] :

For me I want to remplace this definition of the coeficients A[n], it's seem wrong, by  another coeficients depend on n, and x.  The true coeficients A[n] is  a function of n and x.  Because later, I will take x=m*h, withm=-N..N;

In the method used, I don't think, that A[n] change with

How can I change, in the code, the two coeficients A[n] and B[n] to be function of n and x=m*h.

This the lastest version of the code

Fred.mw

Thank you.

@Carl Love 

I modify the code, It's better working now.
First thanks for your remark.

But, In Fredholm_stencil   I have sum(  "int ( f(x), x=n*h..(n+1)*h), n=-N..N-1).

But there is a problem, the dispaly show that the the lower limit of the integral is always nh, doesn't change instead of the upper limit of the integral change according to the change of n=-N..N-1.

Also, I want to simplify the presentation.

Can I add two function i(n)=

and J(n) =

and use them in all my computation.

 The code:

Fred.mw

Thank you very much for any help.

@Carl Love 

Fred.mw

Thank you for your idea.

I have a problem to:

first question. I would  like  recompute (2) by subs the coeficient An and Bn, and then compute (1).
Second question.  Puting x=m*h, in (1),

how can we generate a linear Matrix from (1).

Thank for any help.

@Preben Alsholm 

Thank you. It's okay.

@Carl Love 

Sorry, I delete the original question. Yjis is the original question.

Analytic_numerical_.mw

Please, how can I plot the error beteen the exact and the numerical solution. Thank you

@Carl Love 

Thank for the reference.

How to compute the difference between the two matrix, i.e Numerical solution and analytic solution.

Thanks for your help.

@sarra 

Thanks for all the remarks about this question. 

## The code is very nice, now it's working. For me, I note u[i,j]^k=u(x(i),y(j),t(k));

N:=2;

x := i -> 1/(N)*i;
y := j -> 1/(N)*j;
t := k -> 1/(N)*k;
x[0] = x(0),x[N] = x(N),y[0] = y(0),y[N] = y(N),t[0] = t(0),t[N] = t(N);
#Here is a visualization of the lattice:
Zx := [seq]( x(i),i=0..N);
Zy := [seq]( y(i),i=0..N);
Zt := [seq]( t(i),i=0..N);

Tab:=[seq(seq(seq([Zx[i],Zy[j],Zt[k]],i=1..N+1),j=1..N+1),k=1..N+1)];
pointplot3d(Tab,color = magenta,axes = boxed, symbol = box);

Question: How can I put at each point  (x[i],y[j],t[k]) its name  u[i,j]^k, and  tickmarks in the x and y axis.

I try these lines but I think there is a problem:

textplot([seq([i,-0.1,typeset(u(x+i*h,y,t))],i=-1..1),[0,0.9,typeset(u(x,y+h))],[0,-1.1,typeset(u(x,y-h))]],align={below})],view=[-1.4..1.4,-1.4..1.4],tickmarks=[[-1=x-h,0=x,1=x+h],[-1=y-h,0=y,1=y+h]]

Many thinks.