Advanced Search

acer

32343 Reputation

29 Badges

19 years, 326 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

Contact acer

MaplePrimes Activity


These are answers submitted by acer

dot...

acer 32343
May 31 2017
1 1

What do you see displayed if you replace that * (asterisk) with a . (period)?

something...

acer 32343
May 31 2017
1 1

You might be able to use this kind of approach.

> 

restart

first define input x and output y as vector/list

> 

x := [x1, x2, x3, x4, x5];

[x1, x2, x3, x4, x5]

[y1, y2, y3]

define output values

> 

y1 := exp(-x1/T)*x2+log(x3)*x4+x5^2;

exp(-x1/T)*x2+ln(x3)*x4+x5^2

sin(x1)*x5+x3*x4*x5

c1*x1+c2*x2^2+c3*x3^3*x5+x1*x2*x4*exp(-x1/T)

calculate Jacobian and Hessian matrices

> 

Jac := VectorCalculus:-Jacobian(y, x)

Jac := Matrix(3, 5, {(1, 1) = -exp(-x1/T)*x2/T, (1, 2) = exp(-x1/T), (1, 3) = x4/x3, (1, 4) = ln(x3), (1, 5) = 2*x5, (2, 1) = cos(x1)*x5, (2, 2) = 0, (2, 3) = x4*x5, (2, 4) = x3*x5, (2, 5) = sin(x1)+x3*x4, (3, 1) = c1+x2*x4*exp(-x1/T)-x1*x2*x4*exp(-x1/T)/T, (3, 2) = 2*c2*x2+x1*x4*exp(-x1/T), (3, 3) = 3*c3*x3^2*x5, (3, 4) = x1*x2*exp(-x1/T), (3, 5) = c3*x3^3})

> 

H1 := VectorCalculus:-Hessian(y1, x); H2 := VectorCalculus:-Hessian(y2, x); H3 := VectorCalculus:-Hessian(y3, x)

H1 := Matrix(5, 5, {(1, 1) = exp(-x1/T)*x2/T^2, (1, 2) = -exp(-x1/T)/T, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (2, 1) = -exp(-x1/T)/T, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -x4/x3^2, (3, 4) = 1/x3, (3, 5) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 1/x3, (4, 4) = 0, (4, 5) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 2})

H2 := Matrix(5, 5, {(1, 1) = -sin(x1)*x5, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = cos(x1), (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = x5, (3, 5) = x4, (4, 1) = 0, (4, 2) = 0, (4, 3) = x5, (4, 4) = 0, (4, 5) = x3, (5, 1) = cos(x1), (5, 2) = 0, (5, 3) = x4, (5, 4) = x3, (5, 5) = 0})

H3 := Matrix(5, 5, {(1, 1) = -2*x2*x4*exp(-x1/T)/T+x1*x2*x4*exp(-x1/T)/T^2, (1, 2) = x4*exp(-x1/T)-x1*x4*exp(-x1/T)/T, (1, 3) = 0, (1, 4) = exp(-x1/T)*x2-x1*x2*exp(-x1/T)/T, (1, 5) = 0, (2, 1) = x4*exp(-x1/T)-x1*x4*exp(-x1/T)/T, (2, 2) = 2*c2, (2, 3) = 0, (2, 4) = x1*exp(-x1/T), (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 6*c3*x3*x5, (3, 4) = 0, (3, 5) = 3*c3*x3^2, (4, 1) = exp(-x1/T)*x2-x1*x2*exp(-x1/T)/T, (4, 2) = x1*exp(-x1/T), (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 3*c3*x3^2, (5, 4) = 0, (5, 5) = 0})

> 

CodeGeneration:-C(Jac, optimize);

t1 = 0.1e1 / T;

t2 = x1 * t1;
t3 = exp(-t2);
t4 = (int) (x3 * x3);
cg[0][0] = -t1 * t3 * x2;
cg[0][1] = t3;
cg[0][2] = 0.1e1 / x3 * (double) x4;
cg[0][3] = log(x3);
cg[0][4] = 2 * x5;
cg[1][0] = cos(x1) * (double) x5;
cg[1][1] = 0;
cg[1][2] = x4 * x5;
cg[1][3] = x3 * (double) x5;
cg[1][4] = sin(x1) + x3 * (double) x4;
cg[2][0] = -x2 * (double) x4 * t3 * (-0.1e1 + t2) + c1;
cg[2][1] = x1 * (double) x4 * t3 + 0.2e1 * c2 * x2;
cg[2][2] = 3 * c3 * t4 * x5;
cg[2][3] = x1 * x2 * t3;
cg[2][4] = (double) c3 * x3 * (double) t4;

> 

CodeGeneration:-C(codegen[makeproc]([
                     seq(seq('Jac'[i,j]=Jac[i,j],j=1..5),i=1..3)
                     , 'NULL'],[]), optimize);

#include <math.h>

void cg0 (void)
{
  double t1;
  double t12;
  double t20;
  double t3;
  double t7;
  double t8;
  t1 = 0.1e1 / T;
  t3 = exp(-x1 * t1);
  Jac[0][0] = -t1 * t3 * x2;
  Jac[0][1] = t3;
  Jac[0][2] = 0.1e1 / x3 * x4;
  Jac[0][3] = log(x3);
  Jac[0][4] = (double) (2 * x5);
  t7 = cos(x1);
  Jac[1][0] = t7 * (double) x5;
  Jac[1][1] = 0.0e0;
  Jac[1][2] = x4 * (double) x5;
  Jac[1][3] = x3 * (double) x5;
  t8 = sin(x1);
  Jac[1][4] = x3 * x4 + t8;
  t12 = x1 * x2;
  Jac[2][0] = -t1 * t12 * x4 * Jac[0][1] + x2 * x4 * Jac[0][1] + c1;
  Jac[2][1] = x1 * x4 * Jac[0][1] + 0.2e1 * c2 * x2;
  t20 = x3 * x3;
  Jac[2][2] = 0.3e1 * c3 * t20 * (double) x5;
  Jac[2][3] = t12 * Jac[0][1];
  Jac[2][4] = c3 * t20 * x3;
  ;
}

> 

CodeGeneration:-C(codegen[makeproc]([
                     seq(seq('Jac'[i,j]=Jac[i,j],j=1..5),i=1..3)
                     , seq(seq('H1'[i,j]=H1[i,j],j=1..5),i=1..5)
                     , seq(seq('H2'[i,j]=H2[i,j],j=1..5),i=1..5)
                     , seq(seq('H3'[i,j]=H3[i,j],j=1..5),i=1..5)
                     ],[]));

#include <math.h>

double cg1 (void)
{
  Jac[0][0] = -0.1e1 / T * exp(-x1 / T) * x2;
  Jac[0][1] = exp(-x1 / T);
  Jac[0][2] = 0.1e1 / x3 * x4;
  Jac[0][3] = log(x3);
  Jac[0][4] = (double) (2 * x5);
  Jac[1][0] = cos(x1) * (double) x5;
  Jac[1][1] = 0.0e0;
  Jac[1][2] = x4 * (double) x5;
  Jac[1][3] = x3 * (double) x5;
  Jac[1][4] = sin(x1) + x3 * x4;
  Jac[2][0] = c1 + x2 * x4 * exp(-x1 / T) - x1 * x2 * x4 / T * exp(-x1 / T);
  Jac[2][1] = 0.2e1 * c2 * x2 + x1 * x4 * exp(-x1 / T);
  Jac[2][2] = 0.3e1 * c3 * x3 * x3 * (double) x5;
  Jac[2][3] = x1 * x2 * exp(-x1 / T);
  Jac[2][4] = c3 * pow(x3, 0.3e1);
  H1[0][0] = pow(T, -0.2e1) * exp(-x1 / T) * x2;
  H1[0][1] = -0.1e1 / T * exp(-x1 / T);
  H1[0][2] = 0.0e0;
  H1[0][3] = 0.0e0;
  H1[0][4] = 0.0e0;
  H1[1][0] = -0.1e1 / T * exp(-x1 / T);
  H1[1][1] = 0.0e0;
  H1[1][2] = 0.0e0;
  H1[1][3] = 0.0e0;
  H1[1][4] = 0.0e0;
  H1[2][0] = 0.0e0;
  H1[2][1] = 0.0e0;
  H1[2][2] = -pow(x3, -0.2e1) * x4;
  H1[2][3] = 0.1e1 / x3;
  H1[2][4] = 0.0e0;
  H1[3][0] = 0.0e0;
  H1[3][1] = 0.0e0;
  H1[3][2] = 0.1e1 / x3;
  H1[3][3] = 0.0e0;
  H1[3][4] = 0.0e0;
  H1[4][0] = 0.0e0;
  H1[4][1] = 0.0e0;
  H1[4][2] = 0.0e0;
  H1[4][3] = 0.0e0;
  H1[4][4] = 0.2e1;
  H2[0][0] = -sin(x1) * (double) x5;
  H2[0][1] = 0.0e0;
  H2[0][2] = 0.0e0;
  H2[0][3] = 0.0e0;
  H2[0][4] = cos(x1);
  H2[1][0] = 0.0e0;
  H2[1][1] = 0.0e0;
  H2[1][2] = 0.0e0;
  H2[1][3] = 0.0e0;
  H2[1][4] = 0.0e0;
  H2[2][0] = 0.0e0;
  H2[2][1] = 0.0e0;
  H2[2][2] = 0.0e0;
  H2[2][3] = (double) x5;
  H2[2][4] = x4;
  H2[3][0] = 0.0e0;
  H2[3][1] = 0.0e0;
  H2[3][2] = (double) x5;
  H2[3][3] = 0.0e0;
  H2[3][4] = x3;
  H2[4][0] = cos(x1);
  H2[4][1] = 0.0e0;
  H2[4][2] = x4;
  H2[4][3] = x3;
  H2[4][4] = 0.0e0;
  H3[0][0] = -0.2e1 * x2 * x4 / T * exp(-x1 / T) + x1 * x2 * x4 * pow(T, -0.2e1) * exp(-x1 / T);
  H3[0][1] = x4 * exp(-x1 / T) - x1 * x4 / T * exp(-x1 / T);
  H3[0][2] = 0.0e0;
  H3[0][3] = exp(-x1 / T) * x2 - x1 * x2 / T * exp(-x1 / T);
  H3[0][4] = 0.0e0;
  H3[1][0] = x4 * exp(-x1 / T) - x1 * x4 / T * exp(-x1 / T);
  H3[1][1] = 0.2e1 * c2;
  H3[1][2] = 0.0e0;
  H3[1][3] = x1 * exp(-x1 / T);
  H3[1][4] = 0.0e0;
  H3[2][0] = 0.0e0;
  H3[2][1] = 0.0e0;
  H3[2][2] = 0.6e1 * c3 * x3 * (double) x5;
  H3[2][3] = 0.0e0;
  H3[2][4] = 0.3e1 * c3 * x3 * x3;
  H3[3][0] = exp(-x1 / T) * x2 - x1 * x2 / T * exp(-x1 / T);
  H3[3][1] = x1 * exp(-x1 / T);
  H3[3][2] = 0.0e0;
  H3[3][3] = 0.0e0;
  H3[3][4] = 0.0e0;
  H3[4][0] = 0.0e0;
  H3[4][1] = 0.0e0;
  H3[4][2] = 0.3e1 * c3 * x3 * x3;
  H3[4][3] = 0.0e0;
  H3[4][4] = 0.0e0;
  return(H3[4][4]);
}

> 

CodeGeneration:-C(codegen[makeproc]([
                     seq(seq('Jac'[i,j]=Jac[i,j],j=1..5),i=1..3)
                     , seq(seq('H1'[i,j]=H1[i,j],j=1..5),i=1..5)
                     , seq(seq('H2'[i,j]=H2[i,j],j=1..5),i=1..5)
                     , seq(seq('H3'[i,j]=H3[i,j],j=1..5),i=1..5)
                     ,'NULL'],[]), optimize);

#include <math.h>

void cg2 (void)
{
  double t1;
  double t10;
  double t12;
  double t18;
  double t20;
  double t21;
  double t24;
  double t25;
  double t27;
  double t3;
  double t6;
  double t7;
  double t8;
  t1 = 0.1e1 / T;
  t3 = exp(-x1 * t1);
  Jac[0][0] = -t1 * t3 * x2;
  Jac[0][1] = t3;
  t6 = 0.1e1 / x3;
  Jac[0][2] = t6 * x4;
  Jac[0][3] = log(x3);
  Jac[0][4] = (double) (2 * x5);
  t7 = cos(x1);
  Jac[1][0] = t7 * (double) x5;
  Jac[1][1] = 0.0e0;
  Jac[1][2] = x4 * (double) x5;
  Jac[1][3] = x3 * (double) x5;
  t8 = sin(x1);
  Jac[1][4] = x3 * x4 + t8;
  t10 = x2 * x4;
  t12 = x1 * x2;
  Jac[2][0] = -t1 * t12 * x4 * Jac[0][1] + t10 * Jac[0][1] + c1;
  t18 = x1 * x4;
  Jac[2][1] = 0.2e1 * c2 * x2 + t18 * Jac[0][1];
  t20 = x3 * x3;
  t21 = c3 * t20;
  Jac[2][2] = 0.3e1 * t21 * (double) x5;
  Jac[2][3] = t12 * Jac[0][1];
  Jac[2][4] = c3 * t20 * x3;
  t24 = T * T;
  t25 = 0.1e1 / t24;
  H1[0][0] = t25 * Jac[0][1] * x2;
  t27 = t1 * Jac[0][1];
  H1[0][1] = -t27;
  H1[0][2] = 0.0e0;
  H1[0][3] = 0.0e0;
  H1[0][4] = 0.0e0;
  H1[1][0] = H1[0][1];
  H1[1][1] = 0.0e0;
  H1[1][2] = 0.0e0;
  H1[1][3] = 0.0e0;
  H1[1][4] = 0.0e0;
  H1[2][0] = 0.0e0;
  H1[2][1] = 0.0e0;
  H1[2][2] = -0.1e1 / t20 * x4;
  H1[2][3] = t6;
  H1[2][4] = 0.0e0;
  H1[3][0] = 0.0e0;
  H1[3][1] = 0.0e0;
  H1[3][2] = H1[2][3];
  H1[3][3] = 0.0e0;
  H1[3][4] = 0.0e0;
  H1[4][0] = 0.0e0;
  H1[4][1] = 0.0e0;
  H1[4][2] = 0.0e0;
  H1[4][3] = 0.0e0;
  H1[4][4] = 0.2e1;
  H2[0][0] = -t8 * (double) x5;
  H2[0][1] = 0.0e0;
  H2[0][2] = 0.0e0;
  H2[0][3] = 0.0e0;
  H2[0][4] = t7;
  H2[1][0] = 0.0e0;
  H2[1][1] = 0.0e0;
  H2[1][2] = 0.0e0;
  H2[1][3] = 0.0e0;
  H2[1][4] = 0.0e0;
  H2[2][0] = 0.0e0;
  H2[2][1] = 0.0e0;
  H2[2][2] = 0.0e0;
  H2[2][3] = (double) x5;
  H2[2][4] = x4;
  H2[3][0] = 0.0e0;
  H2[3][1] = 0.0e0;
  H2[3][2] = (double) x5;
  H2[3][3] = 0.0e0;
  H2[3][4] = x3;
  H2[4][0] = H2[0][4];
  H2[4][1] = 0.0e0;
  H2[4][2] = x4;
  H2[4][3] = x3;
  H2[4][4] = 0.0e0;
  H3[0][0] = t12 * t25 * x4 * Jac[0][1] - 0.2e1 * t10 * t27;
  H3[0][1] = -t18 * t27 + x4 * Jac[0][1];
  H3[0][2] = 0.0e0;
  H3[0][3] = -t12 * t27 + Jac[0][1] * x2;
  H3[0][4] = 0.0e0;
  H3[1][0] = H3[0][1];
  H3[1][1] = 0.2e1 * c2;
  H3[1][2] = 0.0e0;
  H3[1][3] = x1 * Jac[0][1];
  H3[1][4] = 0.0e0;
  H3[2][0] = 0.0e0;
  H3[2][1] = 0.0e0;
  H3[2][2] = 0.6e1 * c3 * x3 * (double) x5;
  H3[2][3] = 0.0e0;
  H3[2][4] = 0.3e1 * t21;
  H3[3][0] = H3[0][3];
  H3[3][1] = H3[1][3];
  H3[3][2] = 0.0e0;
  H3[3][3] = 0.0e0;
  H3[3][4] = 0.0e0;
  H3[4][0] = 0.0e0;
  H3[4][1] = 0.0e0;
  H3[4][2] = H3[2][4];
  H3[4][3] = 0.0e0;
  H3[4][4] = 0.0e0;
  ;
}

> 

 

Download Jacobian_Hessian_CG_1.mw

printing...

acer 32343
May 30 2017
0 7

This seems harder that I would have guessed at first, for a few esoteric reasons.

I'm supposing that what you want is a specific order of the derivatives in the pretty-printing of the differential equation, and that you don't actually want the expression's internal representation to be altered necessarily. (The latter would affect how the op command pulls out the terms.)

While I suppose that your D^4 in your Question is just a convenient shorthand (to illustrate your point quickly), it's not completely clear to me whether you mean a call to diff or Diff or D.

So here is a procedure that can print differential equations in a sorted manner.

> 

restart;
> 

`print/pre`:=proc(expr::{scalar,scalar=scalar})
  local T,TT;
  uses PDEtools;
  if type(expr,`=`) then return map(procname,expr); end if;
  if type(expr,`+`) then
    T:=sort([op(expr)],(a,b)->difforder(a)>difforder(b));
    if interface(':-typesetting')=':-standard' then
      TT:=subs([:-diff=`tools/gensym`(:-diff),
                :-Diff=`tools/gensym`(:-Diff),
                :-D=`tools/gensym`(:-D)],T);
      `+`(op(TT));
    else
      `+`(op(subs([:-diff=:-%diff,
                   :-Diff=:-%Diff,
                   :-D=`tools/gensym`(:-D)],T)));
    end if
  else
    expr;
  end if;
end proc:
> 

ee1:=diff(f(x),x)+diff(f(x),$(x,4))+diff(f(x),$(x,6))+7=f(t):
ee2:=Diff(f(x),x)+Diff(f(x),$(x,4))+Diff(f(x),$(x,6))+7=f(t):
ee3:=D(f)(x)+(D@@4)(f)(x)+(D@@6)(f)(x)+7=f(t):
ff1:=diff(x(t),t)+diff(x(t),$(t,4))+diff(x(t),$(t,6))+7=f(t):
> 

interface(typesetting=standard):
> 

ee1;

diff(f(x), x)+diff(diff(diff(diff(f(x), x), x), x), x)+diff(diff(diff(diff(diff(diff(f(x), x), x), x), x), x), x)+7 = f(t)

> 

pre(ee1);

pre(diff(f(x), x)+diff(diff(diff(diff(f(x), x), x), x), x)+diff(diff(diff(diff(diff(diff(f(x), x), x), x), x), x), x)+7 = f(t))

> 

ee2;

Diff(f(x), x)+Diff(f(x), x, x, x, x)+Diff(f(x), x, x, x, x, x, x)+7 = f(t)

> 

pre(ee2);

pre(Diff(f(x), x)+Diff(f(x), x, x, x, x)+Diff(f(x), x, x, x, x, x, x)+7 = f(t))

> 

ee3;

(D(f))(x)+((D@@4)(f))(x)+((D@@6)(f))(x)+7 = f(t)

> 

pre(ee3);

pre((D(f))(x)+((D@@4)(f))(x)+((D@@6)(f))(x)+7 = f(t))

> 

interface(typesetting=extended):
> 

ee1;

diff(f(x), x)+diff(diff(diff(diff(f(x), x), x), x), x)+diff(diff(diff(diff(diff(diff(f(x), x), x), x), x), x), x)+7 = f(t)

> 

pre(ee1);

pre(diff(f(x), x)+diff(diff(diff(diff(f(x), x), x), x), x)+diff(diff(diff(diff(diff(diff(f(x), x), x), x), x), x), x)+7 = f(t))

> 

ff1;

diff(x(t), t)+diff(diff(diff(diff(x(t), t), t), t), t)+diff(diff(diff(diff(diff(diff(x(t), t), t), t), t), t), t)+7 = f(t)

> 

pre(ff1);

pre(diff(x(t), t)+diff(diff(diff(diff(x(t), t), t), t), t)+diff(diff(diff(diff(diff(diff(x(t), t), t), t), t), t), t)+7 = f(t))

> 

ee2;

Diff(f(x), x)+Diff(f(x), x, x, x, x)+Diff(f(x), x, x, x, x, x, x)+7 = f(t)

> 

pre(ee2);

pre(Diff(f(x), x)+Diff(f(x), x, x, x, x)+Diff(f(x), x, x, x, x, x, x)+7 = f(t))

> 

ee3;

(D(f))(x)+((D@@4)(f))(x)+((D@@6)(f))(x)+7 = f(t)

> 

pre(ee3);

pre((D(f))(x)+((D@@4)(f))(x)+((D@@6)(f))(x)+7 = f(t))

> 

 

Download sortdiff.mw

The above also seems to work in Maple 2017.0, with the default typesetting level of extended, with the prime and dot derivative notations for diff toggled on. Eg, after say,

Typesetting:-Settings(typesetprime=true):
Typesetting:-Settings(typesetdot=true):

something...

acer 32343
May 29 2017
2 0

We've seen requests along these lines before. Here's one idea (from here).

> 

restart;
> 

twostep:=proc(expr::uneval)
  uses Typesetting;
  (Typeset(expr) =
   subsindets(expr,
     And(name,
         satisfies(t->type(eval(t),
           {constant,
            '`*`'({constant,
                   'specfunc'(anything,
                              Units:-Unit)})}))),
     z->Typeset(EV(z))))
   = combine(eval(expr),':-units');
end proc:
> 

v:=2.256*10^8*Unit(m/s);

225600000.0*Units:-Unit(('m')/('s'))

> 

f:=4.74*10^14*Unit(1/s);

0.4740000000e15*Units:-Unit(1/('s'))

> 

lambda:=twostep(v/f);

(Typesetting:-mfrac(Typesetting:-mi("v"), Typesetting:-mi("f")) = Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mn("2.256000000"), Typesetting:-mo("&InvisibleTimes;"), Typesetting:-msup(Typesetting:-mn("10"), Typesetting:-mn("8"), Typesetting:-msemantics = "^")), Typesetting:-mo("&InvisibleTimes;"), Typesetting:-mfenced(Typesetting:-mfrac(Typesetting:-mi("m"), Typesetting:-mi("s")), open = "&lobrk;", close = "&robrk;", Typesetting:-msemantics = "Unit"))/Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mn("4.740000000"), Typesetting:-mo("&InvisibleTimes;"), Typesetting:-msup(Typesetting:-mn("10"), Typesetting:-mn("14"), Typesetting:-msemantics = "^")), Typesetting:-mo("&InvisibleTimes;"), Typesetting:-mfenced(Typesetting:-mfrac(Typesetting:-mn("1"), Typesetting:-mi("s")), open = "&lobrk;", close = "&robrk;", Typesetting:-msemantics = "Unit"))) = 0.4759493671e-6*Units:-Unit('m')

> 

 


Download 2step.mw

normal...

acer 32343
May 27 2017
0 1

For your examples, how about just,

not ( evalb( normal( expr ) = expr ) )

And of course you could replace that `normal` with something stronger if necessary. Like `simplify `.

typesetting=extended in Maple 2016...

acer 32343
May 26 2017
0 6

In Maple 2016 (or earlier, mostly) a number of special functions are typeset tersely like that, if interface(typesetting) has been set to extended. But, while that setting can be set by the user, that does not cover all the pre-formed Help pages.

For example, using Maple 2016,

> 

interface(typesetting=standard):
> 

LerchPhi(3,4,5);

LerchPhi(3, 4, 5)

> 

interface(typesetting=extended):
> 

LerchPhi(3,4,5);

LerchPhi(3, 4, 5)

> 

 

LerchPhi_2016.mw

In Maple 2017.0, such terse (and obscure, IMHO) pretty-printing of special functions has been turned off, even with the new default of interface(typesetting) as extended . It can, however, be reinstated even on a function-by-function basis, using commands from the Typesetting package. For example,

> 

interface(typesetting=standard):
> 

LerchPhi(3,4,5);

LerchPhi(3, 4, 5)

> 

interface(typesetting=extended):
> 

LerchPhi(3,4,5);

LerchPhi(3, 4, 5)

> 

Typesetting:-QueryTypesetRule("LerchPhi");

{"LerchPhi" = false}

> 

Typesetting:-EnableTypesetRule("LerchPhi"):
> 

LerchPhi(3,4,5);

LerchPhi(3, 4, 5)

> 

Typesetting:-DisableTypesetRule("LerchPhi"):
> 

LerchPhi(3,4,5);

LerchPhi(3, 4, 5)

> 

 

LerchPhi_2017.mw

read from file...

acer 32343
May 25 2017
2 0

If your Maple procedure has been read from a file, then showsource can show you what that source looked like, including comments.

In addition, if the procedure is saved to a .mla Library archive, then even after a restart the source may still be shown using that command.

Note that the argument to showsource is the name of a Maple procedure, not some source file name (string). The name of the source file and any names which get assigned procedures (within that source file) are not necessarily the same or even similar. I just happened to call my source text file f.mpl , which contains the source that assigns a procedure to name f.

> 

restart;
> 

kernelopts(keepdebuginfo=true):
> 

read "C:/TEMP/f.mpl";

proc (x) local res; res := x^2; return res end proc

> 

showsource(f);


--- C:/TEMP/f.mpl -------------------------------------------------------------
1     proc(x)
2       local res;
4,1     # This is my procedure.
5       res:=x^2;
6,2     # This is another comment.
7       return res;

> 

showstat(f);


f := proc(x)
local res;
   1   res := x^2;
   2   return res
end proc

> 

s:="C:/TEMP/nm.mla";

"C:/TEMP/nm.mla"

> 

LibraryTools:-Create(s);
> 

LibraryTools:-Save(f, s);
> 

restart;
> 

kernelopts(keepdebuginfo=true):
> 

s:="C:/TEMP/nm.mla";

"C:/TEMP/nm.mla"

> 

libname:=s,libname:
> 

f(3);

9

> 

showsource(f);


--- C:/TEMP/f.mpl -------------------------------------------------------------
1     proc(x)
2       local res;
4,1     # This is my procedure.
5       res:=x^2;
6,2     # This is another comment.
7       return res;

> 

showstat(f);


f := proc(x)
local res;
   1   res := x^2;
   2   return res
end proc

> 

 

 

Download showsource.mw

one way...

acer 32343
May 25 2017
2 2


 

> 

restart;
> 

col:=proc(ee, c::string)
  uses Typesetting;
  subsindets(Typeset(ee),
             'specfunc'(anything,{mi,mo,mn,ms,mfenced}),
             u->op(0,u)(op(u),':-mathcolor'=c));
end proc:
> 

captext1 := col( "A plot of: ", "green" ):
capother := col( sin(x), "red" ):
captext2 := col( " from ", "green" ),
            col( -Pi, "green" ),
            col( " to ", "green"),
            col( Pi, "green" ):
> 

plot(sin(x), x=-Pi..Pi, caption=typeset(captext1, capother, captext2));

> 

 

 

Download colored_caption.mw

operator form, or alternate function cal...

acer 32343
May 16 2017
1 1 
evalf(Int(Fg, [0..1, 0..2.2], epsilon=5e-10));
Error, (in evalf/int) Cannot obtain the requested accuracy

evalf(Int(Fg, [0..1, 0..2.2], epsilon=1e-9));
                           8.140784249

evalf(Int(Fg, [0..1, 0..2.2], epsilon=1e-6));

                           8.140784161
Or, construct your procedure to return unevaluated when its arguments are not both of type numeric. (That is more robust than using uneval quotes.) 
newFg := proc(x0,y0)
  if not (x0::numeric and y0::numeric) then
    return 'procname'(args);
  end if;
  if (x0>=0)and(x0<=3) and (y0<=x0+2)
   and (y0>=x0-1) and (y0>=0) and (y0 <=3) then
    return y0*(3-y0)*x0*(3-x0)*(x0+2-y0)*(y0-x0+1);
  else
    return 0;
  end if:
end proc:

evalf(Int(newFg(x,y), [x=0..1, y=0..2.2]));

                                8.140784249
Or, do not use the nested form. This seems to be more friendly to the uneval-quote protection (against premature evaluation). 
evalf(Int(Int('Fg'(x,y), x=0..1), y=0..2.2, epsilon=1e-9));
  
                              8.140784249

local...

acer 32343
May 15 2017
0 8

Remove the semicolon that appears before the keyword `local`.

proc(z0, u0, F0);

I mean that semicolon.

 

mint...

acer 32343
May 13 2017
2 0

If your code is in plaintext then you can run the standalong shell tool mint against it. It has been kept much more up-to-date than has the maplemint command, with regards to newer syntax and Maple language developments.

If you don't want to write out your procedure's source code to a text file (and if you haven't developed it as such, or if running shell scripts is awkward on your platform), and if the rather weak maplemint command cannot handle your procedure, then you can try the method in this old Post. That uses a temporary file to get access to the external mint utility.

one way...

acer 32343
May 10 2017
1 1

It isn't clear to me whether you also want to obtain an explicit result, and whether you want it exact or as a float approximation.

expr:=1/(u+v)^2+4*u*v-1:

ans:=[solve({expr=0, u>0, v>0}, real, explicit)]:
ans:=simplify(combine(evalc(simplify(ans)))):

lprint(ans);

[{u = (2*cos(2/9*Pi)+5/4+3^(1/2)*sin(2/9*Pi))
      *(2*(19-6*cos(4/9*Pi)-12*3^(1/2)*sin(2/9*Pi)
      -6*cos(2/9*Pi))^(1/2)+6*cos(4/9*Pi)
      +3*3^(1/2)*sin(2/9*Pi)-3*cos(2/9*Pi)-2)^(1/2),
  v = 1/4*(2*(19-6*cos(4/9*Pi)-12*3^(1/2)*sin(2/9*Pi)
      -6*cos(2/9*Pi))^(1/2)+6*cos(4/9*Pi)
      +3*3^(1/2)*sin(2/9*Pi)-3*cos(2/9*Pi)-2)^(1/2)}]

evalf(ans);                                      
               [{u = 1.594539597, v = 0.1023339994}]

simplify(eval(expr, ans[1]));                    

                                0

evalf(Sum(...))...

acer 32343
May 03 2017
0 16

It seems to me that if Digits is "high enough" then one may achieve enough numerical robustness to get reasonably close to 10 digits of accuracy with fewer terms in the nested sums. But also, using the accelerated numerical summation techniques offered by evalf(Sum(...)) computes such results quite bit quicker than does add, for these examples.

On my (fast) 64bit Linux machine, using Maple 2016.2, I can compute J(3,Pi/6) to about 9 digits in 16sec using evalf(Sum(...)) while it takes add about 54sec. For the expression a*b*(-A*J(3,Pi/6)+B*J(6,Pi/6)) it computes with evalf(Sum(...) about 7 digits in about a quarter of the time than does add. Your mileage may vary. Please double-check, as we can all make mistakes.

> 

restart;
> 

r := 2.8749: a := 0.7747: b := 0.3812: A := 17.4: B := 29000: R := 5.4813: Z := 2:
> 

JSum := proc(n, phi, NN) options operator, arrow;
8*Pi^(3/2)*r*R*(Sum((2*r*R)^(2*i)*pochhammer((1/2)*n, i)*pochhammer((1/2)*n+1/2, i)*(Sum((-1)^j*cos(phi)^(2*j)*(Sum((2*r*cos(phi))^(2.*l)*pochhammer(n+2*i, 2*l)*hypergeom([2*j+2*l+1, .5], [2*j+2*l+1.5], -1)*(.5*Beta(l+.5, n+2*i+l-.5)-sin(arctan(-Z/sqrt(R^2+r^2)))^(2*l+1)*hypergeom([-n-2*i-l+1.5, l+.5], [l+1.5], sin(arctan(-Z/sqrt(R^2+r^2)))^2)/(2*l+1))/(factorial(2*l)*pochhammer(2*j+2*l+1, .5)*(R^2+r^2)^(n+2*i+l-.5)), l = 0 .. NN))/(factorial(i-j)*factorial(j)), j = 0 .. i))/factorial(i), i = 0 .. NN)) end proc:
> 

Digits:=10: NN:=70:
''Digits''=Digits, ''JSum''(3,Pi/6,NN) = CodeTools:-Usage( evalf( JSum(3,Pi/6,NN) ) ): evalf[10](%);

memory used=2.74GiB, alloc change=36.00MiB, cpu time=15.26s, real time=15.27s, gc time=636.00ms

Digits = 10., JSum(3, (1/6)*Pi, 70) = 15.47952908

> 

restart;
> 

r := 2.8749: a := 0.7747: b := 0.3812: A := 17.4: B := 29000: R := 5.4813: Z := 2:
> 

JSum := proc(n, phi, NN) options operator, arrow;
8*Pi^(3/2)*r*R*(Sum((2*r*R)^(2*i)*pochhammer((1/2)*n, i)*pochhammer((1/2)*n+1/2, i)*(Sum((-1)^j*cos(phi)^(2*j)*(Sum((2*r*cos(phi))^(2.*l)*pochhammer(n+2*i, 2*l)*hypergeom([2*j+2*l+1, .5], [2*j+2*l+1.5], -1)*(.5*Beta(l+.5, n+2*i+l-.5)-sin(arctan(-Z/sqrt(R^2+r^2)))^(2*l+1)*hypergeom([-n-2*i-l+1.5, l+.5], [l+1.5], sin(arctan(-Z/sqrt(R^2+r^2)))^2)/(2*l+1))/(factorial(2*l)*pochhammer(2*j+2*l+1, .5)*(R^2+r^2)^(n+2*i+l-.5)), l = 0 .. NN))/(factorial(i-j)*factorial(j)), j = 0 .. i))/factorial(i), i = 0 .. NN)) end proc:
> 

Digits:=30: NN:=70:
''Digits''=Digits, ''JSum''(3,Pi/6,NN) = CodeTools:-Usage( evalf( JSum(3,Pi/6,NN) ) ): evalf[10](%);

memory used=3.49GiB, alloc change=36.00MiB, cpu time=20.12s, real time=20.14s, gc time=800.00ms

Digits = 30., JSum(3, (1/6)*Pi, 70) = .3995579519

> 

restart;
> 

r := 2.8749: a := 0.7747: b := 0.3812: A := 17.4: B := 29000: R := 5.4813: Z := 2:
> 

Jadd := proc(n, phi, NN) options operator, arrow;
8*Pi^(3/2)*r*R*(add((2*r*R)^(2*i)*pochhammer((1/2)*n, i)*pochhammer((1/2)*n+1/2, i)*(add((-1)^j*cos(phi)^(2*j)*(add((2*r*cos(phi))^(2.*l)*pochhammer(n+2*i, 2*l)*hypergeom([2*j+2*l+1, .5], [2*j+2*l+1.5], -1)*(.5*Beta(l+.5, n+2*i+l-.5)-sin(arctan(-Z/sqrt(R^2+r^2)))^(2*l+1)*hypergeom([-n-2*i-l+1.5, l+.5], [l+1.5], sin(arctan(-Z/sqrt(R^2+r^2)))^2)/(2*l+1))/(factorial(2*l)*pochhammer(2*j+2*l+1, .5)*(R^2+r^2)^(n+2*i+l-.5)), l = 0 .. NN))/(factorial(i-j)*factorial(j)), j = 0 .. i))/factorial(i), i = 0 .. NN)) end proc:
> 

Digits:=10: NN:=70:
''Digits''=Digits, ''Jadd''(3,Pi/6,NN) = CodeTools:-Usage( evalf( Jadd(3,Pi/6,NN) ) ): evalf[10](%);

memory used=3.65GiB, alloc change=32.00MiB, cpu time=20.35s, real time=20.39s, gc time=824.00ms

Digits = 10., Jadd(3, (1/6)*Pi, 70) = -155738.6548

> 

restart;
> 

r := 2.8749: a := 0.7747: b := 0.3812: A := 17.4: B := 29000: R := 5.4813: Z := 2:
> 

Jadd := proc(n, phi, NN) options operator, arrow;
8*Pi^(3/2)*r*R*(add((2*r*R)^(2*i)*pochhammer((1/2)*n, i)*pochhammer((1/2)*n+1/2, i)*(add((-1)^j*cos(phi)^(2*j)*(add((2*r*cos(phi))^(2.*l)*pochhammer(n+2*i, 2*l)*hypergeom([2*j+2*l+1, .5], [2*j+2*l+1.5], -1)*(.5*Beta(l+.5, n+2*i+l-.5)-sin(arctan(-Z/sqrt(R^2+r^2)))^(2*l+1)*hypergeom([-n-2*i-l+1.5, l+.5], [l+1.5], sin(arctan(-Z/sqrt(R^2+r^2)))^2)/(2*l+1))/(factorial(2*l)*pochhammer(2*j+2*l+1, .5)*(R^2+r^2)^(n+2*i+l-.5)), l = 0 .. NN))/(factorial(i-j)*factorial(j)), j = 0 .. i))/factorial(i), i = 0 .. NN)) end proc:
> 

Digits:=30: NN:=70:
''Digits''=Digits, ''Jadd''(3,Pi/6,NN) = CodeTools:-Usage( evalf( Jadd(3,Pi/6,NN) ) ): evalf[10](%);

memory used=5.22GiB, alloc change=36.00MiB, cpu time=28.27s, real time=28.30s, gc time=1.27s

Digits = 30., Jadd(3, (1/6)*Pi, 70) = .3995579519

 

The above are all in separate sessions, to ensure the performance measurements are fair, and absolutely no results
are cached, etc.

 

We can see from the above that using evalf(Sum(....)) gets a result of comparable accuracy to add when Digits=30
and the upper bounds of the summation indices is just 70. But the accelerated floating-point summation method
done by evalf(Sum(...)) is faster than the full explicit summation done by add(...) here.

Let's try and show next that, for at least the case of J(3, Pi/6) , the value of NN=70 suffices when done with Digits=30.
That is to say, more terms, using even greater working precision, is not necessary for at least that example.

 

> 

restart;
> 

r := 2.8749: a := 0.7747: b := 0.3812: A := 17.4: B := 29000: R := 5.4813: Z := 2:
> 

JSum := proc(n, phi, NN) options operator, arrow;
8*Pi^(3/2)*r*R*(Sum((2*r*R)^(2*i)*pochhammer((1/2)*n, i)*pochhammer((1/2)*n+1/2, i)*(Sum((-1)^j*cos(phi)^(2*j)*(Sum((2*r*cos(phi))^(2.*l)*pochhammer(n+2*i, 2*l)*hypergeom([2*j+2*l+1, .5], [2*j+2*l+1.5], -1)*(.5*Beta(l+.5, n+2*i+l-.5)-sin(arctan(-Z/sqrt(R^2+r^2)))^(2*l+1)*hypergeom([-n-2*i-l+1.5, l+.5], [l+1.5], sin(arctan(-Z/sqrt(R^2+r^2)))^2)/(2*l+1))/(factorial(2*l)*pochhammer(2*j+2*l+1, .5)*(R^2+r^2)^(n+2*i+l-.5)), l = 0 .. NN))/(factorial(i-j)*factorial(j)), j = 0 .. i))/factorial(i), i = 0 .. NN)) end proc:
> 

Digits:=100: NN:=100:
''Digits''=Digits, ''JSum''(3,Pi/6,NN) = CodeTools:-Usage( evalf( JSum(3,Pi/6,NN) ) ): evalf[10](%);

memory used=15.43GiB, alloc change=114.57MiB, cpu time=82.20s, real time=80.09s, gc time=4.79s

Digits = 100., JSum(3, (1/6)*Pi, 100) = .3995579519

 

So the result with Digits=30 and NN=100 agrees to 10 digits with that obtained more quickly
with only Digits=30 and NN=70.


You can of course experiment to see whether those faster settings are also numerically robust for
other values of arguments n and phi.

It seems that we can also "simplify" the expression to be evaluated, with some additional performance
gain and only an apparent slight drop in accuracy (for this example at least).

 

> 

restart;
> 

r := 2.8749: a := 0.7747: b := 0.3812: A := 17.4: B := 29000: R := 5.4813: Z := 2:
> 

JSum := proc(n, phi, NN) options operator, arrow;
8*Pi^(3/2)*r*R*(Sum((2*r*R)^(2*i)*pochhammer((1/2)*n, i)*pochhammer((1/2)*n+1/2, i)*(Sum((-1)^j*cos(phi)^(2*j)*(Sum((2*r*cos(phi))^(2*l)*pochhammer(n+2*i, 2*l)*hypergeom([2*j+2*l+1, 1/2], [2*j+2*l+3/2], -1)*(1/2*Beta(l+1/2, n+2*i+l-1/2)-sin(arctan(-Z/sqrt(R^2+r^2)))^(2*l+1)*hypergeom([-n-2*i-l+3/2, l+1/2], [l+3/2], sin(arctan(-Z/sqrt(R^2+r^2)))^2)/(2*l+1))/(factorial(2*l)*pochhammer(2*j+2*l+1, 1/2)*(R^2+r^2)^(n+2*i+l-1/2)), l = 0 .. NN))/(factorial(i-j)*factorial(j)), j = 0 .. i))/factorial(i), i = 0 .. NN)) end proc:
> 

raw:=JSum(ii, Pi/6, NNN):
> 

Digits:=30:
new:=simplify(simplify(combine(convert(combine(raw),StandardFunctions))),size) assuming i::nonnegint, j::nonnegint, l::nonnegint, j<=i, ii::posint, NNN::nonnegint:
> 

NN:=70:
''Digits''=Digits, ''JSum''(3,Pi/6,NN) = CodeTools:-Usage( evalf( subs([ii=3,NNN=NN],new) ) ): evalf[10](%);

memory used=3.06GiB, alloc change=4.00MiB, cpu time=16.58s, real time=16.03s, gc time=1.18s

Digits = 30., JSum(3, (1/6)*Pi, 70) = .3995579521

> 

CodeTools:-Usage( evalf( a * b * ( -A * subs([ii=3,NNN=NN,W=Pi/6],new)
                                   + B * subs([ii=6,NNN=NN,W=Pi/6],new) ) ) ): evalf[10](%);

memory used=3.13GiB, alloc change=0 bytes, cpu time=16.92s, real time=16.38s, gc time=1.21s

.6541253939

> 

Digits:=50:
new:=simplify(simplify(combine(convert(combine(raw),StandardFunctions))),size) assuming i::nonnegint, j::nonnegint, l::nonnegint, j<=i, ii::posint, NNN::nonnegint:
> 

NN:=100:
''Digits''=Digits, ''JSum''(3,Pi/6,NN) = CodeTools:-Usage( evalf( subs([ii=3,NNN=NN],new) ) ): evalf[10](%);

memory used=10.66GiB, alloc change=0 bytes, cpu time=55.85s, real time=53.69s, gc time=4.84s

Digits = 50., JSum(3, (1/6)*Pi, 100) = .3995579521

> 

CodeTools:-Usage( evalf( a * b * ( -A * subs([ii=3,NNN=NN,W=Pi/6],new)
                                   + B * subs([ii=6,NNN=NN,W=Pi/6],new) ) ) ): evalf[10](%);

memory used=10.63GiB, alloc change=40.00MiB, cpu time=55.66s, real time=53.61s, gc time=4.57s

.6541254115

> 

restart;

r := 2.8749; a := 0.7747; b := 0.3812; A := 17.4; B := 29000; R := 5.4813; Z := 2;

J := proc (n, phi) options operator, arrow; 8*Pi^(3/2)*r*R*(add((2*r*R)^(2*i)*pochhammer((1/2)*n, i)*pochhammer((1/2)*n+1/2, i)*(add((-1)^j*cos(phi)^(2*j)*(add((2*r*cos(phi))^(2.*l)*pochhammer(n+2*i, 2*l)*hypergeom([2*j+2*l+1, .5], [2*j+2*l+1.5], -1)*(.5*Beta(l+.5, n+2*i+l-.5)-sin(arctan(-Z/sqrt(R^2+r^2)))^(2*l+1)*hypergeom([-n-2*i-l+1.5, l+.5], [l+1.5], sin(arctan(-Z/sqrt(R^2+r^2)))^2)/(2*l+1))/(factorial(2*l)*pochhammer(2*j+2*l+1, .5)*(R^2+r^2)^(n+2*i+l-.5)), l = 0 .. 70))/(factorial(i-j)*factorial(j)), j = 0 .. i))/factorial(i), i = 0 .. 70)) end proc:

Digits:=30:

CodeTools:-Usage( evalf(a*b*(-A*J(3, (1/6)*Pi)+B*J(6, (1/6)*Pi))) ): evalf[10](%);

2.8749

.7747

.3812

17.4

29000

5.4813

2

memory used=10.44GiB, alloc change=36.00MiB, cpu time=57.16s, real time=57.21s, gc time=2.74s

.6541253936

> 

 

Download sumproc.mw

method=alternate...

acer 32343
May 03 2017
1 2

A little poking around in the code reveals an undocumented method="alternate" option, which makes Maple 2016's plots:-contourplot command utilize plots:-implicitplot internally (via the internal procedure Plot:-ContourPlot).

restart;
kernelopts(version);

   Maple 2016.2, X86 64 LINUX, Jan 13 2017, Build ID 1194701

plots:-contourplot( 1/(x^2+y^2), x=-2 .. 2, y=-2 .. 2, method="alternate" );

plots:-contourplot( 1/(x^2+y^2), x=-2 .. 2, y=-2 .. 2, method="alternate",
                    contours=[seq(i,i=0.5..4,0.5)] );

Download alternate_contourplot.mw

It would be even nicer if that internal procedure Plot:-ContourPlot could separate out additionally passed options specific to plots:-implicitplot (so that it didn't pass all of argument _rest to plots:-display, but allowed options like say gridrefine, etc, to be passed along to plots:-implicitplot). Or, for consistency with plots:-inequal it might be made to pass along with an optionsimplicit option.

evalf[n]...

acer 32343
May 03 2017
2 5

evalf[5]( evalf( expression ) )

will round the result to 5 digits. But the inner call to evalf will still use the current Digits working precision (which you may want done).

First 196 197 198 199 200 201 202 Last Page 198 of 336