## 270 Reputation

8 years, 280 days

## @mehdi jafari    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.

## It's okay. Question solved...

Now, It's okay that is I need. Many thanks.

## @sarra      The seq of al...

The seq of all point is:
x := [seq]( (1/5)*i,i=0..12);  #  x[i] the x-coordinate

y := [seq]( (1/5)*j,j=0..12); # y[j] the y-coordinate

t := [seq]( (1/5)*k,k=0..12);  #  t[k] the t-coordinate

# This table contains all the point must be ploted.

Tab:=seq(seq(seq([x[i],y[j],t[k]],i=1..3),j=1..3),k=1..3);

Now, how can we put these point in the space (x,y,t).

## @mehdi jafari    Thank you for...

Thank you for your remark. But, it's seem for me that your plot is only the point P[i] with coordinate (x[i],y[i],t[i]).  But for me, I want all the point P[i,j,k] with coordinate (x[i],y[j],t[k]).

For example, when you put k=1, you will get all the point  P[i,j,1]  with coordinate (x[i],y[j],t[1]) for all i,j and those point is in the horizontal plane fixed by t=t[1].

I think we will get many horizontal surface caracterized by, t=t[k], and in these plane we put all the point P(i,j,k=fixed).

## Good idea...

Thank you very much, There is a problem of stability. I just change the value of c=1/2, I have not any problem. I get the same figure.

the condition of stability is : c1^2+c2^2<=1. In my old file, I work with c1=1 and c2=1, so that the method is unstable. For example, one van take  c=1/2 , we have not any problem.

This night  I will try to finish the code of stability.

one again, thank you very much.

Please, if I want to plot the error ( between the exact and numerical ), how can I make the difference between the two matrix and plot the error at each time

## @Carl Love I hope it's working.Anal...

Analytic_numerical_s.mw

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

But, my solution is u(x,y,t), when i fixe t from Vect_T:=seq(k*0.1,k=0..10):   and i plot the two solution it's seem for some t, the two curves are the same, but if one take for exemple t=Vect_T[5], i.e working with t=0.4 the two curves are not the same also, take t=t=Vect_T[10], the same conclusion.  Please, run the following code.

Analytic_numerical_s.mw

## @Preben Alsholm    Thank you v...

Thank you very much for the answer.  It is goal.

Many thinks.

Now, I will writte a code to compute the exact solution, then compute the error.

## @mehdi jafari  First special thanks...

First special thanks for your remarks.

Le me recall that our function is a function with three varibles (x,y,t)  so the curve is a surface in R^3 and not in R^2. Can we plot the curve in R^3 and Not in R^2.

u(i,j,k) is an approximation of u(x,y,t) at point (xi,yj,tk) can we plot all the point in R^3.

For example, we can fix two times, (i.e fix k)  and plot the solution for two different times.

There is two differents curves for two time in the last page in the link.

How can I have the same figure.

Thank you very much for your help.

Many thanks.

## @mehdi jafari  Thank you for your a...

Thank you for your answer, I agrre with you,  I forget to put the value of c.  Recall that "u(x,y,t)" is the exact solution of my problem. After discretisation (finite difference):

u[i,j,k] is an approximation  of the analytic  solution  at point (x(i), y(j),t(k)) with

i=1..Fx; j=1..Fy, k=1..Ft.   I want to plot this Numerical  solution.

If I use Matlab, for me I can plot the Numerical solution, for exmaple I can do:

for i= 2:Ft + 1
figure(1);
mesh(y, x, u(:,:,i));
title(‘Numerical solution’);
end

Thanks.

## @Carl Love  I make the change, but ...

I make the change, but now my code run without stop, afther five minutes I Interrupte the computation, maye be a collocation memory needed, or a big error let the code run without stopping.

restart:
ode := diff(y(x), x) = 2*x-y(x);
f:=(x,y)->2*x-y;
analyticsol := rhs(dsolve({ode, y(0) = 1}));
d
--- y(x) = 2 x - y(x)
dx
(x, y) -> 2 x - y
-2 + 2 x + 3 exp(-x)
#starting and ending values of a<x<b and local error tolerance
a:=0; b:=1; epsilon:=1e-2;
## Procedure
test := proc (f, a, b, epsilon, N)
local n, x, y,k,z,R,h,hmin, hmax,hstar,p;
## Initial condition
x[0] := 0; y[0] := 1; z[0]:=1;
R[0]:=0; h[0]:=evalf(b-a)/N; hmax:=h[0];
hmin:=(1e-4)*hmax;

##calculate the RK steps
p:=2;
for n from 0 to N-1 do
x[n+1]:= x[n]+(n+1)*h[n];
k[1] := f(x[n], y[n]);
k[2] := f(x[n]+h[n], y[n]+h[n]*k[1]);
k[3] := f(x[n]+h[n]/2, y[n]+h[n]/4*(k[1]+k[2]));
## 2-stage runge Kutta.
z[n+1] := z[n]+(h[n]/2)*(k[1]+k[2]);
y[n+1] := y[n]+(h[n]/6)*(k[1]+k[2]+4*k[3]);

while x[n]<b and h[n]>hmin do
n:=n;    #set last point to the end point
if x[n]+h[n]>b then
h[n]:=b-x[n];
end
end do;

## local Error
R[n+1]:=abs(y[n+1]-z[n+1]);
hstar[n+1]:=(epsilon/R[n+1])^(1/p);
if R[n+1]<=epsilon
then
h[n+1]:=h[n]; # we take the same stepsize
x[n] :=x[n]+h[n];
z[n] :=z[n+1];
y[n] :=y[n+1];
end if;
if epsilon <=R[n+1] then h[n+1]= min(0.9*hstar[n+1],hmax); end if;
if x[n]<b  then
printf("Integration did not reach end point requested");
end if;
end do;
#return [seq([x[i], z[i]], i = 0 .. N)];
#return[seq([x[i], y[i]], i = 0 .. N)];
return[seq([x[i], R[i]], i = 0 .. N)];
end proc:

#test(f,0,1,1e-2,4);
test(f,0,1,1e-3,4);
Warning,  computation interrupted

## @Carl Love  I used the same Method ...

I used the same Method here, in Runge Kutta with setpsize varibale, but it's run with stop.  I think maybe must I allocate memory or something else.

This is the code.

restart:
ode := diff(y(x), x) = 2*x-y(x);
f:=(x,y)->2*x-y;
analyticsol := rhs(dsolve({ode, y(0) = 1}));
d
--- y(x) = 2 x - y(x)
dx
(x, y) -> 2 x - y
-2 + 2 x + 3 exp(-x)
#starting and ending values of a<x<b and local error tolerance
a:=0; b:=1; epsilon:=1e-2;
## Procedure
test := proc (f, a, b, epsilon, N)
local n, x, y,k,z,R,h,hmin, hmax,hstar,p;
## Initial condition
x[0] := 0; y[0] := 1; z[0]:=1;
R[0]:=0; h[0]:=evalf(b-a)/N; hmax:=h[0];
hmin:=(1e-4)*hmax;

##calculate the RK steps
p:=2;
for n from 0 to N-1 do
x[n+1]:= x[n]+(n+1)*h[n];
k[1] := f(x[n], y[n]);
k[2] := f(x[n]+h[n], y[n]+h[n]*k[1]);
k[3] := f(x[n]+h[n]/2, y[n]+h[n]/4*(k[1]+k[2]));
## 2-stage runge Kutta.
z[n+1] := z[n]+(h[n]/2)*(k[1]+k[2]);
y[n+1] := y[n]+(h[n]/6)*(k[1]+k[2]+4*k[3]);

while x[n]<b and h[n]>hmin do
n:=n;    #set last point to the end point
if x[n]+h[n]>b then
h[n]:=b-x[n];
end
end do;

## local Error
R[n+1]:=abs(y[n+1]-z[n+1]);
hstar[n+1]:=(epsilon/R[n+1])^(1/p);
if R[n+1]<=epsilon
then
h[n+1]:=h[n]; # we take the same stepsize
x[n] :=x[n]+h[n];
z[n] :=z[n+1];
y[n] :=y[n+1];
end if;
if epsilon <=R[n+1] then h[n+1]= min(0.9*hstar[n+1],hmax); end if;
if x[n]<b  then
printf("Integration did not reach end point requested");
end if;
end do;
#return [seq([x[i], z[i]], i = 0 .. N)];
#return[seq([x[i], y[i]], i = 0 .. N)];
return[seq([x[i], R[i]], i = 0 .. N)];
end proc:
0
1
0.01
#test(f,0,1,1e-2,4);
test(f,0,1,1e-3,4);

## Mistake, can correct my code,...

Hi. Please  I need to correct my code.

restart:
ode := diff(y(x), x) = 2*x-y(x);
f:=(x,y)->-y;
analyticsol := rhs(dsolve({ode, y(0) = 1}));
d
--- y(x) = 2 x - y(x)
dx
F:(x, y) -> 2*x-y

#starting and ending values of a<x<b and local error tolerance
a:=0; b:=1; epsilon:=1e-2;
## Procedure
RKadaptivestepsize := proc (f, a, b, epsilon, N)
local n, x, y,k,z,R,h,hmin, hmax,hstar,p;
## Initial condition
x[0] := a; y[0] := 1; z[0]:=1;
R[0]:=0; h:=evalf(b-a)/N; hmax:=h;
hmin:=(1e-4)*hmax;

while (x[n])<b) &(h>hmin)
n:=n+1;    %set last point to the end point
if x[n]+h[n]>b
h[n]:=b-x[n];
end
end
##calculate the RK steps
p:=2;
for n from 0 to N-1 do
x[n+1] := x[n]+(n+1)*h[n];
k[1] := f(x[n], y[n]);
k[2] := f(x[n]+h[n], y[n]+h[n]*k[1]);
k[3] := f(x[n]+h[n]/2, y[n]+h[n]/4*(k[1]+k[2]));
## 2-stage runge Kutta.
z[n+1] := z[n]+(h[n]/2)*(k[1]+k[2]);
y[n+1] := y[n]+(h[n]/6)*(k[1]+k[2]+4*k[3]);
## local Error
R[n+1]:=abs(y[n+1]-z[n+1]);
hstar[n+1]:=sqrt(epsilon/R[n+1]);
if R[n+1]<=epsilon
then
h[n+1]:=h[n]; # we take the same stepsize
x[n]:=x[n+1];
z[n]:=z[n+1];
y[n]:=y[n+1];
end if
if R[n+1]:=0 then h[n+1]:= min(0.9*hstar[n+1],hmax);  end if
if (x[n]<b)
display('Integration did not reach end point requested');
end if;
end do;
#return [seq([x[i], z[i]], i = 0 .. N)];
#return[seq([x[i], y[i]], i = 0 .. N)];
return[seq([x[i], R[i]], i = 0 .. N)];
end proc:

Error, `)` unexpected

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

Thank you for your remark. I agree with you.  I make the change.

aaa:=proc(x,y,f,N)
local bb,i;
for i  from 1 to N do

bb[i]:=f(x[i],y);

if bb[i] <=2  then  y:=3*x[i]  else  y:=x[i]^2

end if ;

end do:

end proc:

## @Carl Love  I work in the interval ...

I work in the interval [0,1].

In my code, when I run it, The points x[n+1]=x[n]+n*h, dont gives all the nodes in my interval, the bigest nodes is O.34. Why. But If I remove the condition if ( conditon about R<=epsilon) I get all the nodes. Please see this problem.

 5 6 7 8 9 10 Page 7 of 10
﻿