## 270 Reputation

8 years, 292 days

## @mehdi jafari  Hi Mehdi, sorry, I ...

Hi Mehdi,

sorry, I was tried last days.

## Thanks for all the remarks....

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.

Hi Alsholm.

## @Carl Love  Thank you for remark....

Thank you for remark.

## @Carl Love Hi Carl,Many thanks for ...

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 you...

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 th...

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:

"(&PartialD;)/(&PartialD; t)u(t,x)=d((&PartialD;)^(2))/(&PartialD;

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;
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);
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.

## I have this code using animation in time...

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

animate_curve.mw

## change coeficient to a function....

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.

## code change, it's better. I have simple ...

I modify the code, It's better working now.

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.

## ok for subs...

Fred.mw

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 ok...

Thank you. It's okay.

## @Carl Love    Sorry, I delete ...

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...

Thank for the reference.

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

## @sarra    Thanks for al...

## 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.

 4 5 6 7 8 9 10 Page 6 of 10
﻿