Maple 2015 Questions and Posts

These are Posts and Questions associated with the product, Maple 2015

Hello;

I am facing "Error, too many levels of recursion" in loops. kindly guide me. Thanks

 

1D1P.mw


 

restart; printlevel := 3

restart; with(LinearAlgebra); with(linalg); with(CodeGeneration); with(plots); Digits := 30; `ε` := 0.1e-1; k := .5; m := 1; c := 1; q := 1

L[time] := 1:

.200000000000000000000000000000

 

.200000000000000000000000000000

 

.200000000000000000000000000000

(1)

``

ICff := 1.0*exp(-p^2/(2.0))*(1.0+epsilon*cos(1.0*k*x))/sqrt(2.0*Pi):

for i from 0 while i <= N[x] do for j from 0 while j <= N[p] do f[0, i, j] := eval(ICff, [x = i*`&Delta;x`, p = j*`&Delta;p`]) end do end do:

for n from 0 while n <= T do for j from 0 while j <= N[p] do f[n, 0, j] := 0; f[n, N[x], j] := 0 end do end do:

for n from 0 while n <= T do for i from 0 while i <= N[x] do f[n, i, 0] := 0; f[n, i, N[p]] := 0 end do end do:

``

for n1 from 0 while n1 <= T-1 do for i1 while i1 <= N[x] do for j1 while j1 <= N[p] do f[n1+1, i1, j1] := f[n1, i1, j1]-`&Delta;t`*j1*`&Delta;p`*(f[n1+1, i1+1, j1]-f[n1+1, i1-1, j1]+f[n1, i1+1, j1]-f[n1, i1-1, j1])/((4*`&Delta;x`)*(1.0)) end do end do end do

Error, too many levels of recursion

 


 

(Maple 2015)

Why is m1(t) not evaluated correctly when I use a compound condition?

restart:
m2(t) = piecewise(x1(t) < c and v1(t) < 0, 1,  0);
       m2(t) = piecewise(x1(t) < e2 and v1(t) < 0, 1, 0)

m1(t) = piecewise(x1(t) = c and v1(t) > 0, 0,  1);
                           m1(t) = 1

This doesn't happen for simple conditions
 

m2(t) = piecewise(x1(t) < c, 1,  0);
m1(t) = piecewise(x1(t) = c, 0,  1);
               m2(t) = piecewise(x1(t) < c, 1, 0)
               m1(t) = piecewise(x1(t) = c, 0, 1)

TIA

Dear Users!

Hope everyone is fine here. Let me explain my problem first for this consider
diff(Y(xi), xi) = mu*(1-Y(xi)^2)
Then the derivative of a function U=u(Y(xi)) using chain rule (and expression menstiones as red) is given as,
diff(U, xi) = (diff(diff(Y, xi), Y))*U and (diff(diff(Y, xi), Y))*U = mu*(1-Y(xi)^2)*(diff(U, Y))
Similarly the second-order derivaitve of U=u(Y(xi)) using chain rule (and expression menstiones as red) is given as,
((&DifferentialD;)^(2))/(&DifferentialD; xi^(2))U=(&DifferentialD;)/(&DifferentialD; xi)(mu (1-Y^(2)(xi))*(&DifferentialD;)/(&DifferentialD; Y)U)=((&DifferentialD;)/(&DifferentialD; Y)*(&DifferentialD;)/(&DifferentialD; xi)Y)(mu (1-Y^(2)(xi))*(&DifferentialD;)/(&DifferentialD; Y)U)=(&DifferentialD;)/(&DifferentialD; Y)(mu^(2) (1-Y^(2)(xi))^(2)*(&DifferentialD;)/(&DifferentialD; Y)U)=-2 Y(xi) mu^(2) (1-Y^(2)(xi))*(&DifferentialD;)/(&DifferentialD; Y)U+ mu^(2) (1-Y^(2)(xi))^(2)*((&DifferentialD;)^(2))/(&DifferentialD; Y^(2))U;
In the similar way I want to compute the higher-order (like 5th order) derivaitve of U w.r.t. xi using the chain rule  (and expression menstiones as red) explained in above. Kindly help me soolve my problem

I am waiting for positive response.

Dear Users!

Hope you are doing well. I have a funtion give bellow:
beta[1]*exp(x*alpha[1]+y*beta[1]-z*sqrt(-alpha[1]^2-beta[1]^2))/(1+exp(x*alpha[1]+y*beta[1]-z*sqrt(-alpha[1]^2-beta[1]^2)));
For any value of alpha[1] and beta[1] the term highlighted red becomes the imaginary form. I want to separate the real and imaginary parts of this function. Kindly help me in this matter, thanks

I want to compute the series expansion of i3_r wrt (x, y, z) at point (x=y=z=0):

i2   := (x,y) -> -(1/2)*I*(exp(I*x)*(sin(x)/x)-exp(I*y)*(sin(y)/y))/(x-y):
i3_r := -(1/2)*I*(i2(y,z)-i2(y,x))/(z-x);

My first attempt was to compute this mulltiple series expansion this way:

ordre := 3:
sx := convert( series(i3_r, x, ordre), polynom);
sy := convert( series(sx  , y, ordre), polynom);
sz := convert( series(sy  , z, ordre), polynom);

But this gives me sy=sz=0 whatever the expansion order.

I then do this:

sx :=              convert(series(i3_r , x, ordre), polynom):
sy := add(map(u -> convert(series(u    , y, ordre), polynom), [op(expand(sx))])):
sz := add(map(u -> convert(series(u    , z, ordre), polynom), [op(expand(sy))]));

and obtained non zero results for both sy and sz (but are they are correct ?).

Could you explain me what happens and tell me how to find the series expansion of i3_r wrt (x, y, z) ?

TIA

I am trying to adapt the code with Direct Search 2 for my needs but do not quite understand the syntaxis :( . I have system with 8 first order but nonlinear odes and 20 parameters but for starters I fix most of them (to estimate time consumption for full calculations etc). The part of the code where it ``breaks'' looks as follows:

#argument
X0:=Vector([24, 36, 48, 60, 72, 84], datatype = float):
#value
Y0:=Vector([13.8e+9, 12.13333333e+9, 10.7e+9, 10.5e+9, 9.8e+9, 10.1e+9], datatype = float):

#so this is the data - time in X0 and corresponding observables (at that time instances) at Y0


par_sol:=dsolve(dsys, numeric, range=0..85,  parameters=[eta_L, L0]);

#parametric numerical solution with eta_L parameter of the model and L0 initial condition to fit

ffit2 := proc(etaLValue, L0Value, tValue) global old_eta, old_L0; local res;
if not [etaLValue, L0Value,tValue]::list(numeric) then
      return 'procname'(args); end if;
if old_eta<>etaLValue and old_L0<>L0Value then
      (old_eta,old_L0) := etaLValue,L0Value;
      par_sol('parameters'=[eta_L=etaLValue,L0=L0Value]);
   end if;
   res:=rhs(par_sol(tValue)[4]); end proc;


# this is just adapted copy of your function

 

par_sol(parameters=[eta_L=0.004, L0=1.1e10]);

par_sol(84)[4];ffit2(0.004, 1.1e10, 84);

#their output coincide OK; same for another trio of parameters

DS2:=CodeTools:-Usage(
        DirectSearch:-DataFit(ffit2,X0,Y0,t,
                              [ eta_L=0.001 .. 0.1, L0=1.5e10 .. 1.7e10
 #                               
                              ] )
);

#and here it breaks:

Warning, initial point [L0Value = .9, tValue = .9] does not satisfy the inequality constraints; trying to find a feasible initial point
Error, (in DirectSearch:-Search) cannot find feasible initial point; specify a new one


But tValue is a time variable - I integrate over it, it is not a parameter to fit - I evaluate ffit2 at it. So basically ffit2 return exactly the model value which need to be fited with data in Y0 at times specified in X0. I thought it is easier than the other problems I saw in examples with Direct Search 2 and still I can't get through syntax neither with it not with NonlinearFit of core Maple...

I would appreciate if you could point me what I am doing wrong. Just in case I also attach full version of the file.01test.mw

Thank you in advance!

Hey everyone!

I am trying to merge two vectors into a matrix. Let us say that both of the vectors has only one line and 4 columns:
A:=[a1,a2,a3];B:=[b1,b2,b3,b4];
Therefore, the matrix that I am trying to get should have this form:
M=Matrix([[a1,b1,0,0,0,0,0,0],[0,0,a2,b2,0,0,0,0],[0,0,0,0,a3,b3,0,0],[0,0,0,0,0,0,a4,b4]]);
Thanks in advance!

restart; N := 100; dx := evalf(2*Pi/N)

100

 

0.6283185308e-1

(1)

f1 := proc (x) options operator, arrow; sin(x) end proc;

proc (x) options operator, arrow; sin(x) end proc

(2)

DiscX := proc (N, dx) local i, xv; xv := Vector(N); for i to N do xv[i] := evalf((i-(1/2)*N-1)*dx) end do; return xv end proc:

Xfun := proc (f1) local i, xa1, xa2; xa1 := Vector(N); xa2 := Vector(N); for i to N do xa1[i] := evalf(subs(x = a[i], f1(x))) end do; return xa1 end proc:

IntNum := proc (N, a, c) local i, xv1; xv1 := Vector(N); for i from 2 to N-1 do xv1[i] := evalf((1/2)*(a[i]-a[i+1])*(c[i]+c[i+1])) end do; return xv1 end proc:

a := DiscX(N, dx); a[1]; a[100]

a := Vector(4, {(1) = ` 1 .. 100 `*Vector[column], (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 

-3.141592654

 

3.078760801

(3)

c := Xfun(f1); c[1]

c := Vector(4, {(1) = ` 1 .. 100 `*Vector[column], (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 

0.4102067615e-9

(4)

k := IntNum(N, a, c)

k := Vector(4, {(1) = ` 1 .. 100 `*Vector[column], (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

(5)

plot([seq([a[i], k[i]], i = 2 .. N)], style = line, title = typeset("Integ_of_sin(x)"), titlefont = [times, bold, 30]);

 

``

Download test.mw

Hey everyone!

I would like to plot the integral of a discrete function. For simplicity, I choose the function sin(x) between –Pi and Pi, which has as an integral -cos(x). I tried to implement that in Maple using the Trapezoidal rule, but the result is simply wrong. Any help would be much appreciated!

I need to declare a whole set of variables as local. The variable names are generates algorithmically using assign. Like so:

seq(seq(assign(cat(S,i,j)=Vector(datatype=float)),i=1..9),j=1..9);

Stand-alone, this works and creates all these Vectors for later use. But this:

local seq(seq(assign(cat(S,i,j)=Vector(datatype=float)),i=1..9),j=1..9);

does not work; I get an "error; '(' unexpected".

I really do not want to type all these by hand... on the other hand, if I do not declare these as local I get 99 warnings about implicit local declaration; not nice.

Is there a way to do this?

Thanks,

M.D.

PS: I do not upload as the one line really is all that is needed. At the lowest level one does not get the implicit-declaration warning, but with "local" it still fails.

hello,my friends.I got a problem in using maple.

I want project a cuboid on a plane, and integrate function on that area.

for example, a cuboid which have a spatial diagonal between (0,0,0) and (1,2,1).
Then it was projected to a plane defined by three points (-1,0,0), (0,-1,0) and (0,0,-1).
Finally, a function x^2+y^2 was integrated on that projected area.

I found the cuboid can be defined with plottools:-cuboid and projected with plottools:-project in maple.
but how to integration next?
I can't find any commands in maple can do the integration with that irregular area.

I need some help.Thank you very much in advance.

Dear Useres!

Hope everyone is fine here! I want to compare the coeficient of exp(k*eta[3]+m*eta[1]+n*eta[2]) for k=0,1,2,3,...,N,n=0,1,2,3,...,N and m=0,1,2,3,...,N for N=10 in the following attached file. But I got some error, please have a look and try to fix it as early as possible. Please take care and thanks

Compare_coeff.mw

I am interested in the inner workings of SignalProcessing:-Convolution. I know I can list it with a higher setting of verboseproc:

with(SignalProcessing);
interface(verboseproc=3); # actually, 2 is enough here...
eval(Convolution);

and get

Obviously the real work happens in IPP:-Convolution, but that seems unknown. How can I list that??

TIA,

Mac Dude

Convolution.mw

I want to estimate numerically the value of P (that is the probability that f exceeds the value 12 when x, y and z are uniformly distributed within the box [-Pi, Pi]3).

restart
f := sin(x)+7*sin(y)^2+0.1*z^4*sin(x);
Omega := [x, y, z] =~ [(-Pi..Pi)$3]:
                                2        4       
               sin(x) + 7 sin(y)  + 0.1 z  sin(x)


h := Heaviside(f-12):
P := Int(h, Omega);


Here is a simple Monte Carlo estimation of P.

f_MC  := x -> sin(x[1])+7*sin(x[2])^2+0.1*x[3]^4*sin(x[1]); 
h_MC  := x -> Heaviside(f_MC(x) - 12);
omega := -Pi, Pi;
                         
P_MC := proc(N)
  local Z := Statistics:-Sample(Uniform(omega), [N, 3]):
  local F := Vector(N, n -> h_MC(Z[n])):
  local K := add(F):
  local P := evalf(add(F)/N):
  local q := Statistics:-Quantile(Normal(0, 1), 0.005, numeric):
  local e := -q*sqrt(P*(1-P)/N):
  printf("A bilateral 99%% confidence interval that Prob(f > 12) is %1.3e +/- %1.2e\n", P, e);
end proc:

P_MC(10^6)
A bilateral 99% confidence interval that Prob(f > 12) is 1.643e-02 +/- 3.27e-04

 

I was hoping to get a value for P using  evalf/Int with some method.
But I didn't succeed with any of the methods for multiple integration:

  • The following command returns 0 for n=1 and 2 and Int(h, Omega) if n >=3
    evalf( Int(h, Omega, 'epsilon=1e-n', 'method=_MonteCarlo') );
    

     

  • And all my attempts with methods from the Cuba library have resulted in an unevaluated Int(h, Omega) integral. 
     

Could you help me to estimate P using  evalf/Int ?
TIA

Hi,
Can anyone tell me what's going wrong with int(Heaviside(f-1/2), x=0..1, y=0..1); ?
 

restart:

f := x*y:

JH := int(Heaviside(f-1/2), x=0..1, y=0..1);  #???
                               1
                               -
                               4

 JP := int(piecewise(x*y>1/2, 1, 0), x=0..1, y=0..1); 
                          1   1      
                          - - - ln(2)
                          2   2      

# integrate "1" over the domain where f > 1/2

S := solve(f=1/2, x):
JS := int(1, x=S..1, y=1/2..1);
                          1   1      
                          - - - ln(2)
                          2   2      

Thanks in advance
 

Hi,

I solve numerically an ODE system which depends on 25 parameters.
I want to know the maximum value of the time at which a specific event is triggered, when these parameters vary (independently the one of the other) within a 25 dimensional hyperbox.

To solve this maximization problem, which I assume is local1, I would like to use NLPSolve

The attached zip (I use Mac OSX) contains:

  • a file that gives a full description of the problem (the original ".pages" file and its export as pdf),
  • an m file wich contains the solution procedure plus a few other variables needed to solve the problem,
  • the mw file which contains:
    • the reading of the m file,
    • the procedure (OBJ) which returns the trigerring time,
    • many attemps to find its maximum value using NLPSolve.

None of my attempts at using NLPSolve gave me the expected answer

Could you help me to fix this?

Thanks in advance

mw_and_m_files.zip

1: Initializing a local method (see the explanation file) with different points gave me different optimizers but all of them led to rather close values of the objective function. Thus the problem is either global instead of local, or convergence might not be achieved for all the initialization points (but keep in mind that I'm not interested in the location of the minimizer(s) but in the maximum value of if the time when the event is triggered).

1 2 3 4 5 6 7 Last Page 3 of 69