Advanced Search

vv

14107 Reputation

20 Badges

10 years, 132 days
Contact vv

MaplePrimes Activity


These are answers submitted by vv

You cannot. And you do not need this. Co...

vv 14107
December 03 2018
0 1

You cannot. And you do not need this.
CompressedSparseForm  is supposed to be used for external programs/platforms.

mtaylor...

vv 14107
December 03 2018
2 0
> 

restart;
> 

T:=mtaylor(f(x+u, y+v, z+w), [u,v,w], 4):
> 

T1:=eval(convert(T,diff),[u=fNx/dx, v=(dy-fNy)/dy, w=fNz/dz]):
> 

PDETools[declare](f(x,y,z));  #optional

` f`(x, y, z)*`will now be displayed as`*f

 (1)
> 

T1;

f(x, y, z)+(diff(f(x, y, z), x))*fNx/dx+(diff(f(x, y, z), y))*(dy-fNy)/dy+(diff(f(x, y, z), z))*fNz/dz+(1/2)*(diff(diff(f(x, y, z), x), x))*fNx^2/dx^2+(diff(diff(f(x, y, z), x), y))*fNx*(dy-fNy)/(dx*dy)+(diff(diff(f(x, y, z), x), z))*fNx*fNz/(dx*dz)+(1/2)*(diff(diff(f(x, y, z), y), y))*(dy-fNy)^2/dy^2+(diff(diff(f(x, y, z), y), z))*(dy-fNy)*fNz/(dy*dz)+(1/2)*(diff(diff(f(x, y, z), z), z))*fNz^2/dz^2+(1/2)*(diff(diff(diff(f(x, y, z), x), y), y))*fNx*(dy-fNy)^2/(dx*dy^2)+(diff(diff(diff(f(x, y, z), x), y), z))*fNx*(dy-fNy)*fNz/(dx*dy*dz)+(1/2)*(diff(diff(diff(f(x, y, z), x), z), z))*fNx*fNz^2/(dx*dz^2)+(1/6)*(diff(diff(diff(f(x, y, z), y), y), y))*(dy-fNy)^3/dy^3+(1/2)*(diff(diff(diff(f(x, y, z), y), y), z))*(dy-fNy)^2*fNz/(dy^2*dz)+(1/2)*(diff(diff(diff(f(x, y, z), y), z), z))*(dy-fNy)*fNz^2/(dy*dz^2)+(1/6)*(diff(diff(diff(f(x, y, z), z), z), z))*fNz^3/dz^3+(1/6)*(diff(diff(diff(f(x, y, z), x), x), x))*fNx^3/dx^3+(1/2)*(diff(diff(diff(f(x, y, z), x), x), y))*fNx^2*(dy-fNy)/(dx^2*dy)+(1/2)*(diff(diff(diff(f(x, y, z), x), x), z))*fNx^2*fNz/(dx^2*dz)

 (2)
> 

 

Numeric L...

vv 14107
December 03 2018
2 0

L must be numeric.

restart;
with(IntegrationTools):
L:=5:
simplify(int(f(x), x = 0 .. L*Ts)-Split(int(f(x), x = 0 .. L*Ts), [i*Ts, i = 0 .. L])):lprint(%);
       int(f(x), x = 0 .. 5*Ts)-(int(f(x), x = 0 .. Ts)) -
         (int(f(x), x = Ts .. 2*Ts))-(int(f(x), x = 2*Ts .. 3*Ts))-(int(f(x), x = 3*Ts .. 4*Ts))-(int(f(x), x = 4*Ts .. 5*Ts))
combine(%);

         0

 

Seq...

vv 14107
December 03 2018
0 5

If you just want to parallelize some numerical computations (without int, sum etc, see ?thread-safe) you can use Threads:-Seq.
Example

f := proc(i) local k; 
if   i=1 then    add(1/evalf((k+sin(k))^(k/(k+1))),k=1..5*10^4)
elif i=2 then    add(1/evalf((k+sin(k))^(k/(k+2))),k=1..6*10^4)
elif i=3 then    add(1/evalf((k+sin(k))^(k/(k+3))),k=1..6*10^4)
elif i=4 then    add(1/evalf((k+sin(k))^(k/(k+4))),k=1..3*10^4)
fi
end:

CodeTools:-Usage([Threads:-Seq(f(i), i=1..4)]);
CodeTools:-Usage([seq(f(i), i=1..4)]);  # without threads, for comparison

 

int...

vv 14107
December 03 2018
2 11
So, you want to approximate

 

> 

J:=Int(exp(I*x)/((1+2*sin(x)^2)*(x+I)), x = -infinity .. infinity);

Int(exp(I*x)/((1+2*sin(x)^2)*(x+I)), x = -infinity .. infinity)

 (1)

(I have taken t=1, which was not given a numerical value)

 

 

The integral is very slowly convergent and oscillates.
It must be transformed in order to be able to approximate. We'll use an integration by parts.

 

> 

F1:=int(cos(x)/((1+2*sin(x)^2)),x);

(1/2)*2^(1/2)*arctan(sin(x)*2^(1/2))

 (2)
> 

F2:=int(sin(x)/((1+2*sin(x)^2)),x);

-(1/6)*6^(1/2)*arctanh((1/3)*cos(x)*6^(1/2))

 (3)
> 

J1:=evalf[15](Int(F1/(x+I)^2, x=-infinity..infinity, 'epsilon'=1e-3, method = _d01amc));

0.-.877869090540869*I

 (4)
> 

J2:=evalf[15](Int(F2/(x+I)^2, x=-infinity..infinity, 'epsilon'=1e-3, method = _d01amc));

.507034848870174+0.*I

 (5)
> 

'J'=evalf[3](J1 + I*J2);

J = 0.-.371*I

 (6)
> 

# For a better accuracy more work is needed.

 

It can be shown that:
J = - 0.37103823863605621848784160669915713405866830221926*I

 

 

 

 

A simple and faster proc...

vv 14107
December 01 2018
3 4 
P:=proc(n::nonnegint,k::nonnegint) 
if n<1 or k<1 or k>n then return NULL fi;
if k=1 then return [n] fi;
if k=n then return [1$n] fi;
seq([u[]+~1], u=[P(n-k,k)]), seq([1,u[]], u=[P(n-1,k-1)])
end:

CodeTools:-Usage(nops([P(100, 5)]));
memory used=323.13MiB, alloc change=76.01MiB, cpu time=2.07s, real time=1.94s, gc time=374.40ms
                             38225

CodeTools:-Usage(nops([P(80, 10)]));
memory used=3.28GiB, alloc change=435.74MiB, cpu time=29.73s, real time=23.37s, gc time=9.70s
                             533975

 

sort...

vv 14107
November 30 2018
4 3

sort~(MM);          

vectors and matrices...

vv 14107
November 27 2018
1 0

1. The exponential of a vector v does not exist in maths. If you want to apply exp elementwise then use exp~(v)  or map(exp, v)
==>  < exp(v[1]), exp(v[2]), ... > .  Note that for a square matrix M there is MatrixExponential(M).

2. To integrate a matrix M (elementwise also)  use  map(int, M, t=a..b)  or map(int, M, [x=a..b, y=c..d,...] ).

3. You cannot have in an integral both  vx  and  vx(x,t)  and integrate wrt  vx.

alternatives...

vv 14107
November 27 2018
1 3

LinearMultivariateSystem has a bug.
Use solve or SolveTools:-SemiAlgebraic  instead.

P.S. If you really need LinearMultivariateSystem  then:

sys1:=subsindets(sys, indexed, e -> cat(op(0,e),__, op(e)));
SolveTools:-Inequality:-LinearMultivariateSystem(sys1,indets(sys1));

 

eval vs subs...

vv 14107
November 24 2018
0 0

This form of eval is an "intelligent" subs. In the second example, a subs was possible. In the first one, subs would produce a nonsense, so the expression was left practically untouched (you can check this using lprint). It is the typesetting system which produced the ...|{x=0} display.

simplex...

vv 14107
November 24 2018
0 2 
NonegComb:=proc(M::Matrix,v::Vector)
local n:=upperbound(M,2), L, t, S;
if upperbound(M,1) <> upperbound(v) then error "Incompatible vector" fi;
L:=convert(add(M[..,i]*t[i],i=1..n)-v,list);
S:=simplex:-minimize(0, [(L=~0)[],  seq(t[i]>=0, i=1..n)]);
if S={} then return false fi;
eval([seq(t[i],i=1..n)],  S);
end:

M:=Matrix([[3, 0, 2, 2, 1, 1], [-1, -1, 0, -1, 0, -1], [0, 3, 0, 1, 1, 2], [3, 0, 1, 2, 0, 1], [-1, -1, -1, -1, -1, -1]]):
v1:= Vector[column](5, [3, 1, 0, 0, -2]):
v2 := M[..,1]+7*M[..,2]+9*M[..,4]:
NonegComb(M,v1);
                             false
NonegComb(M,v2);
                     [0, 13/2, 0, 21/2, 0, 0]

 

seq...

vv 14107
November 24 2018
3 0 
P:=n->seq(seq([a,b,n-a-b],b=max(a, floor(n/2)-a+1)..floor((n-a)/2)), a=1..floor(n/3)):

P(36); nops([%]);
[2,17,17],[3,16,17],[4,15,17],[4,16,16],[5,14,17],[5,15,16],[6,13,17],[6,14,16],[6,15,15],[7,12,17],[7,13,16],[7,14,15],[8,11,17],[8,12,16],[8,13,15],[8,14,14],[9,10,17],[9,11,16],[9,12,15],[9,13,14],[10,10,16],[10,11,15],[10,12,14],[10,13,13],[11,11,14],[11,12,13],[12,12,12]
      27

PolyhedralSets...

vv 14107
November 22 2018
1 1

I have also tried SolveTools:-SemiAlgebraic  and I had to interrupt it.
Probably the decomposition is very large.
The problem can be managed with PolyhedralSets.

restart;
ineqs:= [-a1-a2+a3+a4 <= 2, a1+a2 <= 3, -2*a1 <= -1, -2*a2 <= -1, -2*a3 <= -1, -2*a4 <= -1, -a1+a2+a3-a4 <= 4, -a1+a2-a3+a4 <= 4, a1-a2+a3-a4 <= 4, a1-a2-a3+a4 <= 4]:
with(PolyhedralSets):
ps := PolyhedralSet(ineqs):
vertices:=VerticesAndRays(ps)[1];

vertices := [[1/2, 5/2, 1/2, 5/2], [1/2, 5/2, 3/2, 7/2], [3/2, 3/2, 1/2, 9/2], [1/2, 5/2, 5/2, 1/2], [1/2, 5/2, 7/2, 3/2], [3/2, 3/2, 9/2, 1/2], [1/2, 5/2, 1/2, 1/2], [5/2, 1/2, 3/2, 7/2], [5/2, 1/2, 1/2, 5/2], [5/2, 1/2, 1/2, 1/2], [5/2, 1/2, 7/2, 3/2], [5/2, 1/2, 5/2, 1/2], [1/2, 3/2, 1/2, 7/2], [1/2, 3/2, 7/2, 1/2], [1/2, 1/2, 5/2, 1/2], [1/2, 1/2, 1/2, 5/2], [1/2, 1/2, 1/2, 1/2], [3/2, 1/2, 1/2, 7/2], [3/2, 1/2, 7/2, 1/2]]

So, there are 19 vertices. The solution of your system can be expressed as a convex combination with 19 parameters
(actually 18 if  t1+...+t19=1  is used).
Conclusion: the problem is not as simple as it seems.

 

Workaround...

vv 14107
November 21 2018
1 0

The problem is more general. E.g.    int(f(alpha),alpha=0..1) * L;   latex(%);
Workaround:

restart;
MyLI_:=eval(`latex/int`):
`latex/int` := proc() MyLI_(args)," " end:

 

SemiAlgebraic...

vv 14107
November 20 2018
1 2

 

Equations containing roots and parameters are difficult.
In general Maple will give generic results, which could be valid for some values of the parameters.
I said (only) could, because (mainly in the real case), this is not guaranteed. For your equation even the "obvious" solution alpha=1/2 is not general because for l=2 the denominator is 0.

 

To solve completely such equations in Maple, the best way is to convert them into polynomial equalities and inequalities and then use SolveTools:-SemiAlgebraic.

Let's do it for the numerator only and ignoring alpha=1/2.

 

> 

restart;
> 

f:=(6*alpha^4*l^2-7*alpha^3*l^2+6*alpha^3*l+2*alpha^2*l^2-6*alpha^2*l+alpha*l+3*alpha-2*sqrt(alpha^3*l^2*(alpha*l-l+1)*(9*alpha^4*l-13*alpha^3*l+9*alpha^3+6*alpha^2*l-12*alpha^2-alpha*l+6*alpha-1))-1);

6*alpha^4*l^2-7*alpha^3*l^2+6*alpha^3*l+2*alpha^2*l^2-6*alpha^2*l+alpha*l+3*alpha-2*(alpha^3*l^2*(alpha*l-l+1)*(9*alpha^4*l-13*alpha^3*l+9*alpha^3+6*alpha^2*l-12*alpha^2-alpha*l+6*alpha-1))^(1/2)-1

 (1)
> 

f1:=select(type,f,polynom);
f2:=f1-f;

6*alpha^4*l^2-7*alpha^3*l^2+6*alpha^3*l+2*alpha^2*l^2-6*alpha^2*l+alpha*l+3*alpha-1

  

2*(alpha^3*l^2*(alpha*l-l+1)*(9*alpha^4*l-13*alpha^3*l+9*alpha^3+6*alpha^2*l-12*alpha^2-alpha*l+6*alpha-1))^(1/2)

 (2)
> 

sys:=simplify([f1^2-f2^2, f1>=0])[];

(alpha*l+1)*(4*alpha^2*l-3*alpha*l+1)*(alpha^2*l-3*alpha+1)^2, 0 <= -1+6*alpha^4*l^2+(-7*l^2+6*l)*alpha^3+(2*l^2-6*l)*alpha^2+(l+3)*alpha

 (3)
> 

SolveTools:-SemiAlgebraic([sys, l>=0], parameters=[l]);

piecewise(l < 0, [], l = 0, [[alpha = 1/3]], l < 16/9, [[alpha = -(-3+sqrt(-4*l+9))/(2*l)], [alpha = (3+sqrt(-4*l+9))/(2*l)]], l = 16/9, [[alpha = 27/32-3*sqrt(17)*(1/32)], [alpha = 27/32+3*sqrt(17)*(1/32)]], l < 2, [[alpha = -(-3+sqrt(-4*l+9))/(2*l)], [alpha = (3+sqrt(-4*l+9))/(2*l)]], l = 2, [[alpha = 1], [alpha = 1/2]], l < 9/4, [[alpha = -(-3+sqrt(-4*l+9))/(2*l)], [alpha = (3+sqrt(-4*l+9))/(2*l)]], l = 9/4, [[alpha = 2/3]], l < 9/2, [], l = 9/2, [[alpha = 2/3]], 9/2 < l, [[alpha = (3*l+sqrt(9*l^2-16*l))/(8*l)]])

 (4)

I have added l>=0  for easier visualization. You can omit l>=0  or solve the system

> 

SolveTools:-SemiAlgebraic([sys, l<0], parameters=[l]);

but the solution is very large in this case and it is not listed here.

 

 

 

 

First 61 62 63 64 65 66 67 Last Page 63 of 121