Dear all

Please, I need you help, to make a simple procedure with two output.  I try this simple code :

test:= proc(a,b)

local variable .....;

for i from 1 to 10 do

expression

end do

retur  fff;

return ggg;

end proc

Then I do  : test(a,b) but I get only ggg and not fff

Any help please

Dear all,

Thank you for your Help.

h: stepsize;

x in [0,x0];

I give all the step of my code, but I think there is a mistake. I wait for your Help.

I would like to compute the error between  Method Huen with step size h and step size 2h using the definition of epsilon given below:

 ## The error written epsilon(x0,h)= sqrt(1/(N+1) * sum_i=0^N  (y_i^{2h}-y_(2i)^h)^2 ), where y_i^(2h) is the approximation of y(i*2*h).

 ## We want : loglog epsilon versus h.

  epsilon:=(x0,h)->sqrt( 1/(N+1)*add( (Heun1(f,x0,i)-Heun2(f,x0,i))^2,i=0..N ) );

  f:=(x,y)=1/(1+cos(y)); 

  ode:=diff(y(x),x)=f(x,y);

ic:=y(0)=1;  h:=x0/(2*N);

## Method Heun with step size 2h

> Heun1 := proc (f, x0,)

local x, y, i, h, k;

y := Array(0 .. N);

x := Array(0 .. N);

h := evalf((1/2)*x0/N);

x[0] := 0;

y[0] := 1;

for i from 0 to N do

x[i+1] := (2*i+2)*h;

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

k[2] := f(x[i]+h, y[i]+h*k[1]);

y[i+1] := y[i]+h*((1/2)*k[1]+(1/2)*k[2]);

end do;

[seq([x[i], y[i]], i = 0 .. N)];

end proc;

### Now Heun with step size h  ( the same h)

> Heun2 := proc (f, x0,)

local x, y, i, h, k;

y := Array(0 .. N);

x := Array(0 .. N);

h := evalf((1/2)*x0/N);

x[0] := 0;

y[0] := 1;

for i from 0 to N do

x[i+1] := (i+1)*h;

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

k[2] := f(x[i]+h, y[i]+h*k[1]);

y[i+1] := y[i]+h*((1/2)*k[1]+(1/2)*k[2]);

end do;

[seq([x[2*i], y[2*i]], i = 0 .. N)];

end proc;

Thanks you for your help.

Dear all

I want to display the value of k   for which the stop test is verified. Thanks for your help. 

epsilon:=(y0,N)->evalf(abs(

 testerror:=array(1..75);
for k from 2 by 2 to 50 do
testerror[k]:=data[RK3][k][2];
if testerror[k]<=0.0001 then break; end if; end do;

Dear all;

Than you for help.

how  many steps are required to achieve a error of 1.e-3 in the numerical value of y(1).

Here The 3 -step procedure  Range Kutta Method.

## Exact  solution  

### We will modifty N ( number of steps to get error =10^(-3). )

## Procedure Range Kutta

> RK3 := proc (f, a, b, y0, N)

local x, y, n, h, k, vectRK3;

y := Array(0 .. N);

x := Array(0 .. N);

h := evalf(b-a)/N;

x[0] := a; y[0] := 1;

for n from 0 to N-1 do

x[n+1] := a+(n+1)*h;

k[1] := f(x[n], y[n]);

k[2] := f(x[n]+(1/2)*h, y[n]+(1/2)*h*k[1]);

k[3] := f(x[n]+h, y[n]+h*(-k[1]+2*k[2]));

y[n+1] := y[n]+(1/6)*h*(k[1]+4*k[2]+k[3])

end do;

[seq([x[n], y[n]], n = 0 .. N)]; y[1];

end proc;

## Now  we compute the error between y(1) and exact  solution for different value of  N

### I have a problem in this part

 errorRk3 := array(1 .. 29);
 for N from  2 to 30 do

errorrRk3[N] := abs(eval(rhs(res), x = 1)-RK3((x,y)->-y,0,1,N));

if errorrRk3[N] =10^{-3} end ;
end  do ;

Dear all,

I need your help

I have a problem with dchange ....  use the dchange command to transform pde to the
                                 x[1], x[2].

> restart;
> with(PDEtools); with(LinearAlgebra);
> pde := diff(u(t), t, t)+2*GAMMA*(diff(u(t), t))+omega^2*u(t) = 0;
           / d  / d      \\           / d      \        2         
           |--- |--- u(t)|| + 2 GAMMA |--- u(t)| + omega  u(t) = 0
           \ dt \ dt     //           \ dt     /                  
> deq1 := diff(u(t), t) = v(t);
                                d             
                               --- u(t) = v(t)
                                dt            
>
> deq2 := subs(deq1, pde);
                 / d      \                       2         
                 |--- v(t)| + 2 GAMMA v(t) + omega  u(t) = 0
                 \ dt     /                                 
>
> dsolve({deq1, deq2}, {u(t), v(t)});
>
> eqns := [rhs(deq1) = lhs(deq1), rhs(deq2) = lhs(deq2)];
       [        d            / d      \                       2     ]
       [v(t) = --- u(t), 0 = |--- v(t)| + 2 GAMMA v(t) + omega  u(t)]
       [        dt           \ dt     /                             ]
> y := [u, v]; b := diff(y(t), t);
                                   [u, v]
                            [ d         d      ]
                            [--- u(t), --- v(t)]
                            [ dt        dt     ]
> A, b := GenerateMatrix(eqns, y(t));
          Matrix(%id = 122038892), Vector[column](%id = 135944696)
 # Return a vector of eigenvalue of A and matrix  whose columns are eigenvectors of A
> gnat := Eigenvectors(A);
> lambda := gnat[1]; Lambda := gnat[2];
                       Vector[column](%id = 135975976)
                           Matrix(%id = 136787860)
> Y := Vector([y]);
                       Vector[column](%id = 123771808)
> tr := solve(GenerateEquations(Lambda, [x[1], x[2]], Y), {u, v});
     /           /                                    (1/2)             
     |      1    |                   /     2        2\                  
    < u = ------ \-x[1] GAMMA - x[1] \GAMMA  - omega /      - x[2] GAMMA
     |         2                                                        
     \    omega                                                         

                               (1/2)\                 \
              /     2        2\     |                 |
       + x[2] \GAMMA  - omega /     /, v = x[1] + x[2] >
                                                      |
                                                      /
>
> dchange(tr, pde, [x[1], x[2]]);