Advanced Search

vv

13837 Reputation

20 Badges

9 years, 319 days
Contact vv

MaplePrimes Activity


These are answers submitted by vv

int...

vv 13837
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 13837
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 13837
November 30 2018
4 3

sort~(MM);          

vectors and matrices...

vv 13837
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 13837
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 13837
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 13837
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 13837
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 13837
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 13837
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 13837
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.

 

 

 

 

A faster solution...

vv 13837
November 20 2018
2 1

Here is a faster solution.
Note that in Acer's solution a,b  must be a bit larger not to miss solutions.

> 

restart;
> 

P:=isolve(u^2+v^2 = 225):
> 

Q:=isolve(u^2+v^2 = 125):
> 

NS:=0:
for p in [P] do for q in [Q] do
  NS:=NS+1;
  s:=eval([u,v],p); t:=eval([u,v],q);
  x-a=s[1],y-b=s[2],x+1=t[1],y+3=t[2];
  S[NS]:= eval([a,b,x,y],solve({%}))
od od:
S:=sort(convert(S,list)):NS;

192

 (1)
> 

nr:=0:NT:=0:
for i to NS-1 do
if S[i][[1,2]] = S[i+1][[1,2]] then eq:=true; nr:=nr+1 else eq:=false fi;
if eq=false then
   if nr=1 then NT:=NT+1;  lprint(NT=S[i-1],S[i])fi;
   nr:=0
fi;
od:

1 = [-21, -13, -12, -1], [-21, -13, -6, -13]
2 = [-21, 7, -12, -5], [-21, 7, -6, 7]
3 = [-15, -5, -6, 7], [-15, -5, -3, -14]
4 = [-15, -1, -6, -13], [-15, -1, -3, 8]
5 = [-11, -23, -11, -8], [-11, -23, 1, -14]
6 = [-11, -13, -11, 2], [-11, -13, 4, -13]
7 = [-11, 7, -11, -8], [-11, 7, 4, 7]
8 = [-11, 17, -11, 2], [-11, 17, 1, 8]
9 = [-8, -4, 1, 8], [-8, -4, 4, -13]
10 = [-8, -2, 1, -14], [-8, -2, 4, 7]
11 = [-6, -8, -6, 7], [-6, -8, 9, -8]
12 = [-6, 2, -6, -13], [-6, 2, 9, 2]
13 = [-5, -5, 4, 7], [-5, -5, 10, -5]
14 = [-5, -1, 4, -13], [-5, -1, 10, -1]
15 = [-3, -17, -12, -5], [-3, -17, 9, -8]
16 = [-3, -7, -3, 8], [-3, -7, 9, 2]
17 = [-3, 1, -3, -14], [-3, 1, 9, -8]
18 = [-3, 11, -12, -1], [-3, 11, 9, 2]
19 = [-2, -10, -11, 2], [-2, -10, 10, -1]
20 = [-2, 4, -11, -8], [-2, 4, 10, -5]
21 = [0, -10, -12, -1], [0, -10, 9, 2]
22 = [0, 4, -12, -5], [0, 4, 9, -8]
23 = [1, -17, -11, -8], [1, -17, 10, -5]
24 = [1, -7, -11, 2], [1, -7, 1, 8]
25 = [1, 1, -11, -8], [1, 1, 1, -14]
26 = [1, 11, -11, 2], [1, 11, 10, -1]
27 = [3, -5, -12, -5], [3, -5, -6, 7]
28 = [3, -1, -12, -1], [3, -1, -6, -13]
29 = [4, -8, -11, -8], [4, -8, 4, 7]
30 = [4, 2, -11, 2], [4, 2, 4, -13]
31 = [6, -4, -6, -13], [6, -4, -3, 8]
32 = [6, -2, -6, 7], [6, -2, -3, -14]
33 = [9, -23, -3, -14], [9, -23, 9, -8]
34 = [9, -13, -6, -13], [9, -13, 9, 2]
35 = [9, 7, -6, 7], [9, 7, 9, -8]
36 = [9, 17, -3, 8], [9, 17, 9, 2]
37 = [13, -5, 1, -14], [13, -5, 4, 7]
38 = [13, -1, 1, 8], [13, -1, 4, -13]
39 = [19, -13, 4, -13], [19, -13, 10, -1]
40 = [19, 7, 4, 7], [19, 7, 10, -5]

  
> 

 


 

Download aisolve.mw

Inverse...

vv 13837
November 19 2018
0 1

The inverse function is elementary

W:=solve(Y=vs, w):
Y1:=fsolve(W=0):
Y2:=fsolve(W=10):
plot([W,Y, Y=Y2 .. Y1]);

parfrac...

vv 13837
November 18 2018
1 0

You don't have to compute this way the partial fractions. Simply use:

convert(Denom1, parfrac, s);


 

Basis...

vv 13837
November 18 2018
1 1

You must take the leading monomials of the Groebner basis for J, not of J itself.

First 60 61 62 63 64 65 66 Last Page 62 of 120