Advanced Search

mmcdara

7891 Reputation

22 Badges

9 years, 56 days
Contact mmcdara

MaplePrimes Activity


These are answers submitted by mmcdara

To obtain the MajorAxis and the MinorAxi...

mmcdara 7891
October 22 2022
0 2

To obtain the MajorAxis and the MinorAxis of an ellipse in canonical (cartesian) form just do this

restart
with(geometry):
_EnvHorizontalName := 'X': _EnvVerticalName := 'Y':
ellipse(ell, (9/4)*X^2+9*Y^2 = 1):
MajorAxis(ell), MinorAxis(ell);

Another way (same package) is 
 

restart with(geometry): _EnvHorizontalName := 'X': _EnvVerticalName := 'Y':
conic(c1, (9/4)*X^2+9*Y^2 = 1,[X, Y]):
detail(c1);
                 /                            
   GeometryDetail|["name of the object", c1], 
                 \                            

     ["form of the object", ellipse2d], ["center", [0, 0]], 

     [        [[  1  (1/2)   ]  [1  (1/2)   ]]]  
     ["foci", [[- - 3     , 0], [- 3     , 0]]], 
     [        [[  3          ]  [3          ]]]  

     [                            4]  
     ["length of the major axis", -], 
     [                            3]  

     [                            2]  
     ["length of the minor axis", -], 
     [                            3]  

     [                                9  2      2    ]\
     ["equation of the ellipse", -1 + - X  + 9 Y  = 0]|
     [                                4              ]/

Better if you are nor sure of the type of the conic

To complete Carl's reply...

mmcdara 7891
October 21 2022
0 0

Firstly what you have written in your quesion is not syntactically correct, neither in Maple nor in mathematics.
Secondly, your "problem" is not well defined (2 equations and 7 formal quantities)

I guess you are a newbye in Maple.
So read the help page (search for solve, this should be a good start) to write syntactically correct relations: this will fix the first point baove.
For the second one it's up to you to correctly formalize what you want to do, for neither MAple, nor any other CAS will be able to give you the answer of your ill defined problem.

> 

restart

How many syntax error does your pseudo-code contain ?
(You forgot almost everywhere the `*` sign!

After executing each line look where the pursor is set on the input line to see where the error comes from

> 

{[[b/(sin(x-y-z))=a/(sin(z)),],[(a^(2)+b^(2)+2 abcos(x-y))cos(m(y+z)+n)=acos(mx+n)+bcos(my+n),]];

Error, `]` unexpected

  
> 

{[[b/(sin(x-y-z))=a/(sin(z))],[(a^(2)+b^(2)+2 abcos(x-y))cos(m(y+z)+n)=acos(mx+n)+bcos(my+n),]];

Error, missing operator or `;`

  
> 

{[[b/(sin(x-y-z))=a/(sin(z))],[(a^(2)+b^(2)+2*a*b*cos(x-y))cos(m(y+z)+n)=acos(mx+n)+bcos(my+n),]];

Error, missing operator or `;`

  
> 

{[[b/(sin(x-y-z))=a/(sin(z))],[(a^(2)+b^(2)+2*a*b*cos(x-y))*cos(m(y+z)+n)=acos(mx+n)+bcos(my+n),]];

Error, `]` unexpected

  
> 

{[[b/(sin(x-y-z))=a/(sin(z))],[(a^(2)+b^(2)+2*a*b*cos(x-y))*cos(m(y+z)+n)=acos(mx+n)+bcos(my+n)]];

Error, `;` unexpected

  

Here it is!

> 

[[b/(sin(x-y-z))=a/(sin(z))],[(a^(2)+b^(2)+2*a*b*cos(x-y))*cos(m(y+z)+n)=acos(mx+n)+bcos(my+n)]];

[[b/sin(x-y-z) = a/sin(z)], [(a^2+b^2+2*a*b*cos(x-y))*cos(m(y+z)+n) = acos(mx+n)+bcos(my+n)]]

 (1)

Nevertheless I think that acos is a typo mistake and should be a*cos?
Same thing for bcos, my instead of m*y, et caetera

> 

[[b/(sin(x-y-z))=a/(sin(z))],[(a^(2)+b^(2)+2*a*b*cos(x-y))*cos(m*(y+z)+n)=a*cos(m*x+n)+b*cos(m*y+n)]];

[[b/sin(x-y-z) = a/sin(z)], [(a^2+b^2+2*a*b*cos(x-y))*cos(m*(y+z)+n) = a*cos(m*x+n)+b*cos(m*y+n)]]

 (2)

There is no need for the inner square brackets

> 

sys := [b/(sin(x-y-z))=a/(sin(z)), (a^(2)+b^(2)+2*a*b*cos(x-y))*cos(m*(y+z)+n)=a*cos(m*x+n)+b*cos(m*y+n)];

[b/sin(x-y-z) = a/sin(z), (a^2+b^2+2*a*b*cos(x-y))*cos(m*(y+z)+n) = a*cos(m*x+n)+b*cos(m*y+n)]

 (3)

Two equalities and 7 indeterminates a, b, m, n, x, y, z
Here is a solution w.r.t [a, b]

> 

solve(sys, [a, b])

[[a = 0, b = 0], [a = tan(z)*(-cos(x-y)*tan(z)*cos(m*y+n)+sin(x-y)*cos(m*y+n)+tan(z)*cos(m*x+n))/(sin(x-y)^2*cos(m*y+m*z+n)*(tan(z)^2+1)), b = -(tan(z)^2*sin(x-y)^2*cos(m*y+n)+tan(z)^2*cos(m*x+n)*cos(x-y)+2*tan(z)*sin(x-y)*cos(m*y+n)*cos(x-y)-tan(z)^2*cos(m*y+n)-tan(z)*sin(x-y)*cos(m*x+n)-sin(x-y)^2*cos(m*y+n))/(sin(x-y)^2*cos(m*y+m*z+n)*(tan(z)^2+1))]]

 (4)

Harder: the solution  w.r.t [m, n] (no solution found)

> 

solve(sys, [m, n])

Download Big_problems.mw

I am not very familiar with Python and e...

mmcdara 7891
October 21 2022
1 3

I am not very familiar with Python and even less with the HyperNext package.
So I cannot say if Maple 2022's GraphTheory package offers the same features that  HyperNext does.

But I believe the answer to your last question is "yes", at least if you accept to exchange the information between Maple and Python through files.

I use to use this mechanism to delegate to R (a programming language dedicated to statistics) things that Maple cannot do and recover in Maple the results from R.

Here is an example for Python
Python_Hypergraph_from_Maple.mw

> 

restart:
> 

interface(version)

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

 (1)
> 

with(GraphTheory):
> 

e_nodes := seq(e||i, i=1..4):
v_nodes := seq(v||i, i=1..7):
nodes   := [e_nodes, v_nodes]:
G := Graph(nodes, {[e1, v1], [e1, v2], [e1, v3], [e2, v2], [e2, v3], [e3, v3], [e3, v5], [e3, v6], [e4, v4]})

GRAPHLN(directed, unweighted, [e1, e2, e3, e4, v1, v2, v3, v4, v5, v6, v7], Array(%id = 18446744078138018262), `GRAPHLN/table/15`, 0)

 (2)
> 

DrawGraph(G);

 

The Python code below come from this source
 https://stackoverflow.com/questions/9658373/using-hypergraphs-in-pygraph-need-verification-that-example-is-correct

from pygraph.classes.hypergraph import hypergraph

from pygraph.readwrite.dot import write_hypergraph

 

h = hypergraph()

 

h.add_nodes(['v1', 'v2', 'v3', 'v4', 'v5', 'v6', 'v7'])

h.add_hyperedges(['e1', 'e2', 'e3', 'e4'])

 

h.link('v1', 'e1')

h.link('v2', 'e1')

h.link('v3', 'e1')

h.link('v2', 'e2')

h.link('v3', 'e2')

h.link('v3', 'e3')

h.link('v5', 'e3')

h.link('v6', 'e3')

h.link('v4', 'e4')

 

with open('hypergraph.dot', 'w') as f:

    f.write(write_hypergraph(h))

 

> 

edges := IncidentEdges(G, [e_nodes]);
> 

# how to Maple write the code above?

printf("from pygraph.classes.hypergraph import hypergraph\n"):
printf("from pygraph.readwrite.dot import write_hypergraph\n"):
printf("h = hypergraph()\n"):
printf(cat("h.add_nodes([", seq(cat("'", u, "', "), u in v_edges[1..-2]), cat("'", v_edges[-1], "'"), "])\n"));
printf(cat("h.add_hyperedges([", seq(cat("'", u, "', "), u in e_edges[1..-2]), cat("'", e_edges[-1], "'"), "])\n"));
for edge in edges do
  printf(cat("h.link(", "'", edge[2], "', '", edge[1], "')\n"))
end do;
printf("with open('hypergraph.dot', 'w') as f:\n"):
printf("f.write(write_hypergraph(h))\n")
 

from pygraph.classes.hypergraph import hypergraph
from pygraph.readwrite.dot import write_hypergraph

h = hypergraph()
h.add_nodes(['v1', 'v2', 'v3', 'v4', 'v5', 'v6', 'v7'])
h.add_hyperedges(['e1', 'e2', 'e3', 'e4'])
h.link('v1', 'e1')
h.link('v2', 'e1')
h.link('v3', 'e1')
h.link('v2', 'e2')
h.link('v3', 'e2')
h.link('v3', 'e3')
h.link('v5', 'e3')
h.link('v6', 'e3')
h.link('v4', 'e4')
with open('hypergraph.dot', 'w') as f:
f.write(write_hypergraph(h))

  
> 

# Instead of using printf use fprintf to write these lines in some ".py" script file
# for instance test.py
#
# run Python with script file "test.py"
# ssystem("python3 test.py")

In case you would like to recover in Maple the results produced by Python it'is "enough" to read from Maple the output file created by Python.
"enough" means "provided that the output format is compatible with Maple"

Here is your file corrected.The correcti...

mmcdara 7891
October 20 2022
2 1

Here is your file corrected.
The corrections and/or explanations are written in brown (courier font when you will open the file with Maple).

(1) Question1: why the follwing solve output nothing???
In such a cas I use to use infolevel[the_function_I_use] to try and see what could happen.
Here infolevel[solve] := 4; says that 

Main: Entering solver with 15 equations in 1 variable
Transformer:   solving the uncoupled linear subsystem in {Ia1}
Linear: solving 15 linear equations
Polynomial: # of equations is: 2
Transformer:   system has a inconsistent linear subsystem
Main: solving successful - now forming solutions
Main: Exiting solver returning 0 solutions
solve: Warning: no solutions found

The important message is the first one.
Follow some brown commands which should (IMO) give some insights about the structure of your system and the collection of its "unknowns".

(2) From this quick look, and from what I understand when you write "what I expected is Ia1=f(Va1) which Va1 as indepentent,but the above solution have no Va1 ", I suggest you this:

  • The system of 15 equations contains 16 indeterminates.
  • One of them is Ia1.
  • Give Ia1 the special status of a "parameter".
    ( You could have define equList_Y11all := unapply(equList_Y11all, Ia1)) and solve(equList_Y11all(Ia1)) )
  • Define the "true unknowns" by  
    unknowns := remove(has, var_Y11, Ia1);
  • Solve equList_Y11all with respect to unknowns .
  • Let Y11_all the solution : Select from Y11_all the component of the form Va1 = ... .
  • Check that this component depends on Ia1
  • Isolate Ia1 from this component to get the desired (?) expression Ia1=f(Va1, ...)

     

Series Connect Couple Line Question

 

> 

restart

with(LinearAlgebra): 

 

 

1.the first 4 equation

Ia(4*1)=Ya(4*4)*Va(4*1)

 

Ia := Vector(4, [Ia1, Ia2, Ia3, Ia4])

Vector[column]([[Ia1], [Ia2], [Ia3], [Ia4]])

 (1)

Va := [Va1, Va2, Va3, Va4]; #must be a list!

[Va1, Va2, Va3, Va4]

 (2)

Ya := Matrix([[Ya11, Ya12, Ya13, Ya14], [Ya21, Ya22, Ya23, Ya24], [Ya31, Ya32, Ya33, Ya34], [Ya41, Ya42, Ya43, Ya44]])

Matrix([[Ya11, Ya12, Ya13, Ya14], [Ya21, Ya22, Ya23, Ya24], [Ya31, Ya32, Ya33, Ya34], [Ya41, Ya42, Ya43, Ya44]])

 (3)

equList_Ya := GenerateEquations(Ya, Va, Ia);

[Va1*Ya11+Va2*Ya12+Va3*Ya13+Va4*Ya14 = Ia1, Va1*Ya21+Va2*Ya22+Va3*Ya23+Va4*Ya24 = Ia2, Va1*Ya31+Va2*Ya32+Va3*Ya33+Va4*Ya34 = Ia3, Va1*Ya41+Va2*Ya42+Va3*Ya43+Va4*Ya44 = Ia4]

 (4)

equation(4)is the first 4 equation set

 

 

 

2.the 2nd 4 equation

Ib(4*1)=Yb(4*4)*Vb(4*1)

 

Ib := Vector(4, [Ib1, Ib2, Ib3, Ib4])

Vector[column]([[Ib1], [Ib2], [Ib3], [Ib4]])

 (5)

Vb := [Vb1, Vb2, Vb3, Vb4]; #must be a list!

[Vb1, Vb2, Vb3, Vb4]

 (6)

Yb := Matrix([[Yb11, Yb12, Yb13, Yb14], [Yb21, Yb22, Yb23, Yb24], [Yb31, Yb32, Yb33, Yb34], [Yb41, Yb42, Yb43, Yb44]])

Matrix([[Yb11, Yb12, Yb13, Yb14], [Yb21, Yb22, Yb23, Yb24], [Yb31, Yb32, Yb33, Yb34], [Yb41, Yb42, Yb43, Yb44]])

 (7)

equList_Yb := LinearAlgebra:-GenerateEquations(Yb, Vb, Ib);

[Vb1*Yb11+Vb2*Yb12+Vb3*Yb13+Vb4*Yb14 = Ib1, Vb1*Yb21+Vb2*Yb22+Vb3*Yb23+Vb4*Yb24 = Ib2, Vb1*Yb31+Vb2*Yb32+Vb3*Yb33+Vb4*Yb34 = Ib3, Vb1*Yb41+Vb2*Yb42+Vb3*Yb43+Vb4*Yb44 = Ib4]

 (8)

equation(8)is the 2nd 4 equation set

 

 

3.the 3nd 4 equation

equList_ConnectCondition := [Ia3 = -Ib1, Ia4 = -Ib2, Va3 = Vb1, Va4 = Vb2];

[Ia3 = -Ib1, Ia4 = -Ib2, Va3 = Vb1, Va4 = Vb2]

 (9)

 

 

 

unitl now,we have 12 equations

equList_all := {op(equList_Ya), op(equList_Yb), op(equList_ConnectCondition)}

{Ia3 = -Ib1, Ia4 = -Ib2, Va3 = Vb1, Va4 = Vb2, Va1*Ya11+Va2*Ya12+Va3*Ya13+Va4*Ya14 = Ia1, Va1*Ya21+Va2*Ya22+Va3*Ya23+Va4*Ya24 = Ia2, Va1*Ya31+Va2*Ya32+Va3*Ya33+Va4*Ya34 = Ia3, Va1*Ya41+Va2*Ya42+Va3*Ya43+Va4*Ya44 = Ia4, Vb1*Yb11+Vb2*Yb12+Vb3*Yb13+Vb4*Yb14 = Ib1, Vb1*Yb21+Vb2*Yb22+Vb3*Yb23+Vb4*Yb24 = Ib2, Vb1*Yb31+Vb2*Yb32+Vb3*Yb33+Vb4*Yb34 = Ib3, Vb1*Yb41+Vb2*Yb42+Vb3*Yb43+Vb4*Yb44 = Ib4}

 (10)

equList_all_num := nops(equList_all)

12

 (11)

 

what we want to get is Yparameter

(Vector(4, {(1) = Ia1, (2) = Ia2, (3) = Ib3, (4) = Ib4})) = Typesetting[delayDotProduct](Matrix(4, 4, {(1, 1) = Y11, (1, 2) = Y12, (1, 3) = Y13, (1, 4) = Y14, (2, 1) = Y21, (2, 2) = Y22, (2, 3) = Y23, (2, 4) = Y24, (3, 1) = Y31, (3, 2) = Y32, (3, 3) = Y33, (3, 4) = Y34, (4, 1) = Y41, (4, 2) = Y42, (4, 3) = Y43, (4, 4) = Y44}), Vector(4, {(1) = Va1, (2) = Va2, (3) = Vb3, (4) = Vb4}), true)

 

extra equations to get Y11

equList_Y11Condition := [Va2 = 0, Vb3 = 0, Vb4 = 0]

[Va2 = 0, Vb3 = 0, Vb4 = 0]

 (12)

 

all equations for Y11

equList_Y11all := {op(equList_all), op(equList_Y11Condition)}

{Ia3 = -Ib1, Ia4 = -Ib2, Va2 = 0, Va3 = Vb1, Va4 = Vb2, Vb3 = 0, Vb4 = 0, Va1*Ya11+Va2*Ya12+Va3*Ya13+Va4*Ya14 = Ia1, Va1*Ya21+Va2*Ya22+Va3*Ya23+Va4*Ya24 = Ia2, Va1*Ya31+Va2*Ya32+Va3*Ya33+Va4*Ya34 = Ia3, Va1*Ya41+Va2*Ya42+Va3*Ya43+Va4*Ya44 = Ia4, Vb1*Yb11+Vb2*Yb12+Vb3*Yb13+Vb4*Yb14 = Ib1, Vb1*Yb21+Vb2*Yb22+Vb3*Yb23+Vb4*Yb24 = Ib2, Vb1*Yb31+Vb2*Yb32+Vb3*Yb33+Vb4*Yb34 = Ib3, Vb1*Yb41+Vb2*Yb42+Vb3*Yb43+Vb4*Yb44 = Ib4}

 (13)

 

equList_all_num := nops(equList_Y11all)

15

 (14)

 

16 var，15 equation, so I hope to get Ia1=f(Va1) to calculate Y11=Ia1/Va1

var_Y11 := [Va1, Va2, Va3, Va4, Ia1, Ia2, Ia3, Ia4, Vb1, Vb2, Vb3, Vb4, Ib1, Ib2, Ib3, Ib4]

[Va1, Va2, Va3, Va4, Ia1, Ia2, Ia3, Ia4, Vb1, Vb2, Vb3, Vb4, Ib1, Ib2, Ib3, Ib4]

 (15)

 

Question1: why the follwing solve output nothing???

infolevel[solve] := 4:

Main: Entering solver with 15 equations in 1 variable

Transformer:   solving the uncoupled linear subsystem in {Ia1}
Linear: solving 15 linear equations

Polynomial: # of equations is: 2
Transformer:   system has a inconsistent linear subsystem

Main: solving successful - now forming solutions
Main: Exiting solver returning 0 solutions
solve: Warning: no solutions found

  

 

> 

# 15 equations and 47 indeterminates
# note that 3 equations are trivial : Va2=Vb3=Vb4=0

print~(equList_Y11all):
print():
`number of equations` = numelems(equList_Y11all);
indeterminates := indets(equList_Y11all, name);
`number of indeterminates` = numelems(indeterminates);
 

Ia3 = -Ib1

  

Ia4 = -Ib2

  

Va2 = 0

  

Va3 = Vb1

  

Va4 = Vb2

  

Vb3 = 0

  

Vb4 = 0

  

Va1*Ya11+Va2*Ya12+Va3*Ya13+Va4*Ya14 = Ia1

  

Va1*Ya21+Va2*Ya22+Va3*Ya23+Va4*Ya24 = Ia2

  

Va1*Ya31+Va2*Ya32+Va3*Ya33+Va4*Ya34 = Ia3

  

Va1*Ya41+Va2*Ya42+Va3*Ya43+Va4*Ya44 = Ia4

  

Vb1*Yb11+Vb2*Yb12+Vb3*Yb13+Vb4*Yb14 = Ib1

  

Vb1*Yb21+Vb2*Yb22+Vb3*Yb23+Vb4*Yb24 = Ib2

  

Vb1*Yb31+Vb2*Yb32+Vb3*Yb33+Vb4*Yb34 = Ib3

  

Vb1*Yb41+Vb2*Yb42+Vb3*Yb43+Vb4*Yb44 = Ib4

  

 

`number of equations` = 15

  

{Ia1, Ia2, Ia3, Ia4, Ib1, Ib2, Ib3, Ib4, Va1, Va2, Va3, Va4, Vb1, Vb2, Vb3, Vb4, Ya11, Ya12, Ya13, Ya14, Ya21, Ya22, Ya23, Ya24, Ya31, Ya32, Ya33, Ya34, Ya41, Ya42, Ya43, Ya44, Yb11, Yb12, Yb13, Yb14, Yb21, Yb22, Yb23, Yb24, Yb31, Yb32, Yb33, Yb34, Yb41, Yb42, Yb43, Yb44}

  

`number of indeterminates` = 48

 (16)
> 

member(Ia1, indeterminates)

true

 (17)
> 

# in your command solve the 2nd argument is intended to represent the variables you want to solve the system for
# one problem is that Ia1 is not an indeterminate of the system equList_Y11all)

if member(Ia1, indeterminates) then
  print("Ia1 is an indeterminate of the system"):
else
  print("Ia1 is not an indeterminate of the system"):
end if:
 

print~(equList_Y11all):

"Ia1 is an indeterminate of the system"

  

Ia3 = -Ib1

  

Ia4 = -Ib2

  

Va2 = 0

  

Va3 = Vb1

  

Va4 = Vb2

  

Vb3 = 0

  

Vb4 = 0

  

Va1*Ya11+Va2*Ya12+Va3*Ya13+Va4*Ya14 = Ia1

  

Va1*Ya21+Va2*Ya22+Va3*Ya23+Va4*Ya24 = Ia2

  

Va1*Ya31+Va2*Ya32+Va3*Ya33+Va4*Ya34 = Ia3

  

Va1*Ya41+Va2*Ya42+Va3*Ya43+Va4*Ya44 = Ia4

  

Vb1*Yb11+Vb2*Yb12+Vb3*Yb13+Vb4*Yb14 = Ib1

  

Vb1*Yb21+Vb2*Yb22+Vb3*Yb23+Vb4*Yb24 = Ib2

  

Vb1*Yb31+Vb2*Yb32+Vb3*Yb33+Vb4*Yb34 = Ib3

  

Vb1*Yb41+Vb2*Yb42+Vb3*Yb43+Vb4*Yb44 = Ib4

 (18)
> 

# the 4th message seems to say that solve operates on 12 equations only: remember 
# that 3 are trivial (see above).

unknowns := remove(has, var_Y11, Ia1);
Y11_all := solve(equList_Y11all, unknowns, symbolic = true):

[Va1, Va2, Va3, Va4, Ia2, Ia3, Ia4, Vb1, Vb2, Vb3, Vb4, Ib1, Ib2, Ib3, Ib4]

  

Main: Entering solver with 15 equations in 15 variables
Main: attempting to solve as a linear system
Linear: solving 15 linear equations
Polynomial: # of equations is: 12
Polynomial: best equation / unknown Ib1 Ia3 1
Polynomial: # of equations is: 11
Polynomial: best equation / unknown Ib2 Ia4 1
Polynomial: # of equations is: 10
Polynomial: best equation / unknown -Vb1 Va3 1
Polynomial: # of equations is: 9
Polynomial: best equation / unknown -Vb2 Va4 1
Polynomial: # of equations is: 8
Polynomial: best equation / unknown Vb1*Yb11+Vb2*Yb12 Ib1 -1
Polynomial: # of equations is: 7
Polynomial: best equation / unknown Vb1*Yb21+Vb2*Yb22 Ib2 -1
Polynomial: # of equations is: 6
Polynomial: best equation / unknown Vb1*Yb31+Vb2*Yb32 Ib3 -1
Polynomial: # of equations is: 5
Polynomial: best equation / unknown Vb1*Yb41+Vb2*Yb42 Ib4 -1
Polynomial: # of equations is: 4
Polynomial: best equation / unknown Vb1*Ya13+Vb2*Ya14-Ia1 Va1 Ya11
Polynomial: # of equations is: 3
Polynomial: best equation / unknown (Ya11*Ya23-Ya13*Ya21)*Vb1+(Ya11*Ya24-Ya14*Ya21)*Vb2+Ia1*Ya21 Ia2 -Ya11
Polynomial: # of equations is: 2
Polynomial: best equation / unknown (Ya11*Ya34+Ya11*Yb12-Ya14*Ya31)*Vb2+Ia1*Ya31 Vb1 Ya11*Ya33+Ya11*Yb11-Ya13*Ya31
Polynomial: # of equations is: 1
Polynomial: best equation / unknown -Ia1*Ya31*Ya43-Ia1*Ya31*Yb21+Ia1*Ya33*Ya41+Ia1*Ya41*Yb11 Vb2 Ya11*Ya33*Ya44+Ya11*Ya33*Yb22-Ya11*Ya34*Ya43-Ya11*Ya34*Yb21-Ya11*Ya43*Yb12+Ya11*Ya44*Yb11+Ya11*Yb11*Yb22-Ya11*Yb12*Yb21-Ya13*Ya31*Ya44-Ya13*Ya31*Yb22+Ya13*Ya34*Ya41+Ya13*Ya41*Yb12+Ya14*Ya31*Ya43+Ya14*Ya31*Yb21-Ya14*Ya33*Ya41-Ya14*Ya41*Yb11
Polynomial: backsubstitution at: 12
Polynomial: backsubstitution at: 11
Polynomial: backsubstitution at: 10
Polynomial: backsubstitution at: 9
Polynomial: backsubstitution at: 8
Polynomial: backsubstitution at: 7
Polynomial: backsubstitution at: 6
Polynomial: backsubstitution at: 5
Polynomial: backsubstitution at: 4
Polynomial: backsubstitution at: 3
Polynomial: backsubstitution at: 2
Polynomial: backsubstitution at: 1
Main: Linear solver successful. Exiting solver returning 1 solution

  

 

> 

# let's check the solution

whattype(Y11_all);
numelems(Y11_all);  # consistent with the last displayed message by the solve command

# so let's redefine Y11_all this way:

Y11_all := Y11_all[1]:
whattype(Y11_all);
numelems(Y11_all);

list

  

1

  

list

  

15

 (19)
> 

# what unknown do depend on Ia1?

for k from 1 to numelems(Y11_all) do
  if  indets(Y11_all[k], name) intersect {Ia1} <> {} then
    printf("%a depends on Ia1\n", lhs(Y11_all[k]));
  else
    printf("%a does not depend on Ia1\n", lhs(Y11_all[k]));
  end if:
end do;

Va1 depends on Ia1

Va2 does not depend on Ia1
Va3 depends on Ia1
Va4 depends on Ia1
Ia2 depends on Ia1
Ia3 depends on Ia1
Ia4 depends on Ia1
Vb1 depends on Ia1
Vb2 depends on Ia1
Vb3 does not depend on Ia1
Vb4 does not depend on Ia1
Ib1 depends on Ia1
Ib2 depends on Ia1
Ib3 depends on Ia1
Ib4 depends on Ia1

  
> 

# select Va1 from Y11_all

Va1_expr := select(has, Y11_all, Va1)[]

Va1 = Ia1*(Ya33*Ya44+Ya33*Yb22-Ya34*Ya43-Ya34*Yb21-Ya43*Yb12+Ya44*Yb11+Yb11*Yb22-Yb12*Yb21)/(Ya11*Ya33*Ya44+Ya11*Ya33*Yb22-Ya11*Ya34*Ya43-Ya11*Ya34*Yb21-Ya11*Ya43*Yb12+Ya11*Ya44*Yb11+Ya11*Yb11*Yb22-Ya11*Yb12*Yb21-Ya13*Ya31*Ya44-Ya13*Ya31*Yb22+Ya13*Ya34*Ya41+Ya13*Ya41*Yb12+Ya14*Ya31*Ya43+Ya14*Ya31*Yb21-Ya14*Ya33*Ya41-Ya14*Ya41*Yb11)

 (20)
> 

# isolate Ia1

isolate(Va1_expr, Ia1);

Ia1 = Va1*(Ya11*Ya33*Ya44+Ya11*Ya33*Yb22-Ya11*Ya34*Ya43-Ya11*Ya34*Yb21-Ya11*Ya43*Yb12+Ya11*Ya44*Yb11+Ya11*Yb11*Yb22-Ya11*Yb12*Yb21-Ya13*Ya31*Ya44-Ya13*Ya31*Yb22+Ya13*Ya34*Ya41+Ya13*Ya41*Yb12+Ya14*Ya31*Ya43+Ya14*Ya31*Yb21-Ya14*Ya33*Ya41-Ya14*Ya41*Yb11)/(Ya33*Ya44+Ya33*Yb22-Ya34*Ya43-Ya34*Yb21-Ya43*Yb12+Ya44*Yb11+Yb11*Yb22-Yb12*Yb21)

 (21)
> 

 

NULL

Download CoupleLineSeriesConnect_mmcdara.mw

A way......

mmcdara 7891
October 15 2022
2 0


... waiting for probably more astute proposals:

(be careful, not all H matrices define an ellipse)
ellipse_mmcdara.mw

> 

restart:
> 

H := Matrix(2$2, [3, 1, 1, 1]);
eigv := LinearAlgebra:-Eigenvectors(H)[2];

H := Matrix(2, 2, {(1, 1) = 3, (1, 2) = 1, (2, 1) = 1, (2, 2) = 1})

  

eigv := Matrix(2, 2, {(1, 1) = 1/(2^(1/2)-1), (1, 2) = 1/(-1-2^(1/2)), (2, 1) = 1, (2, 2) = 1})

 (1)
> 

c   := <1, 2>:
xy  := <x, y>:
d   := xy-c;
ell := expand(d^+ . H . d)=1

d := Vector(2, {(1) = x-1, (2) = y-2})

  

3*x^2+2*x*y+y^2-10*x-6*y+11 = 1

 (2)

solve(ell, x)[1]:
select(has, [op(%)], y^2)[]:
y_bounds := solve(%^2, y):

solve(ell, y)[1]:
select(has, [op(%)], x^2)[]:
x_bounds := solve(%^2, x):

ranges := x=x_bounds[1]..x_bounds[2], y=y_bounds[1]..y_bounds[2];

plots:-display(
  plot([convert(c, list)], style=point, symbol=circle, symbolsize=20, color=blue)
  , plots:-implicitplot(ell, ranges, grid=[50, 50])
  , plots:-implicitplot(eigv[..,1].d=0, ranges, linestyle=3, color=blue)
  , plots:-implicitplot(eigv[..,2].d=0, ranges, linestyle=3, color=blue)
  , view=[0..2, 0..4]
  , scaling=constrained
);
 

  
 
> 

 

 

(-2)^x is nor necessarily real....

mmcdara 7891
October 15 2022
2 0

Here is a quick example for 

x := [-2, -3/2, -1, -1/2, 0, 1/2, 1, 3/2, 2]
(-2)^~x;
[1/4, (1/4)*sqrt(-2), -1/2, -(1/2)*sqrt(-2), 1, sqrt(-2), -2, -2*sqrt(-2), 4]

A workaround might be to use plots:-complexplot (sorry I can only insert the link not the content of the file)

Maple Worksheet - Error
Failed to load the worksheet /maplenet/convert/explanation.mw

Download explanation.mw

Here is a proposal where the OperatorFo...

mmcdara 7891
October 14 2022
2 3


Here is a proposal where the OperatorForm is used for NLPSolve.
One interest here is to trace the evolution of the optimization process (I use MAPLE 2015 and it doesn't seem that a trace mode is avliable through some option ???)
(To trace the evolution uncomment the yellow highlighted command)
 

> 

restart:with(LinearAlgebra):
> 

N:=3:
> 

m:=Vector([ 1 , log(x+b3) , b2/(x+b3) ]):
> 

A:=m.m^+:
> 

for i to N do
m||i:=eval(A,[x=x||i]);
od:
> 

M:=add(w||i*m||i,i=1..N-1)+(1-add(w||i,i=1..N-1))*m||N:
> 

MM:=( LinearAlgebra:-Trace(MatrixInverse(M)) ):
> 

IF1 := proc(__w1, __w2, __x1, __x2, __x3)
  local J:
  J := evalf(
         Int(
           eval(MM, [w1=__w1, w2=__w2, x1=__x1, x2=__x2, x3=__x3])
           , [b2=1..2,b3=1..2]
           , method = _d01ajc
           , epsilon=0.001
         )
       ):
  # print([_passed], J);
  return J:
end proc:
> 

# check

IF1(0.6, .1, 8, 7, 5);

2204832.780

 (1)
> 

s:= Optimization:-NLPSolve(
  IF1
  , 0..1, 0..1, 7..10, 5..10, 1..10
  , initialpoint=[0.6, .1, 8, 7, 5]
  , maximize=false
)

s := [482.354032179177977, Vector(5, {(1) = 1.0, (2) = 1.0, (3) = 10.0, (4) = 5.0, (5) = 1.0}, datatype = float[8])]

 (2)
> 

 

 

Download LinearLog-A-Bayesian_1_mmcdara.mw

 

The reason is...

mmcdara 7891
October 14 2022
1 2


Functions 1 and 4 are of the form 

+/-sqrt(-2+2*cos(z))

and are purely imaginary for any valu of z.

It's up to you to know how to fix this :-)



 

Yes, it's possible (edited answer)...

mmcdara 7891
October 13 2022
2 3

s this possible?

  • formally : NO YES
  • numerically YES

Here is my first reply where solve gave no results. I keep the file for it contains a plot you might find useful.

Half_Width_mmcdara.mw

And here is the file where the two points are computed formally:

> 

restart
> 

f := (x, x__0, S, C) -> -C*(exp(S*(x - x__0)^2*(-1/2)) - 1)*Heaviside(x - x__0)

proc (x, x__0, S, C) options operator, arrow; -C*(exp(-(1/2)*S*(x-x__0)^2)-1)*Heaviside(x-x__0) end proc

 (1)
> 

X0   := 0:
fDOG := unapply(f(x, X0, S__a, 1) - f(x, X0, S__d, 1), (x, S__a, S__d))

proc (x, S__a, S__d) options operator, arrow; -(exp(-(1/2)*S__a*x^2)-1)*Heaviside(x)+(exp(-(1/2)*S__d*x^2)-1)*Heaviside(x) end proc

 (2)
> 

DfDOG := diff(fDOG(x, S__a, S__d), x) assuming x > X0

S__a*x*exp(-(1/2)*S__a*x^2)-S__d*x*exp(-(1/2)*S__d*x^2)

 (3)
> 

solve({DfDOG=0, x > X0}, x) assuming S__a > S__d;
tpeak := eval(x, %[1])

[{x = 2^(1/2)*((S__a-S__d)*ln(S__a/S__d))^(1/2)/(S__a-S__d)}, x = x]

  

2^(1/2)*((S__a-S__d)*ln(S__a/S__d))^(1/2)/(S__a-S__d)

 (4)
> 

hpeak := fDOG(tpeak, S__a, S__d) assuming S__a > S__d, S__d > 0

-exp(-S__a*ln(S__a/S__d)/(S__a-S__d))+exp(-S__d*ln(S__a/S__d)/(S__a-S__d))

 (5)
> 

half_hpeak := hpeak/2;

-(1/2)*exp(-S__a*ln(S__a/S__d)/(S__a-S__d))+(1/2)*exp(-S__d*ln(S__a/S__d)/(S__a-S__d))

 (6)
> 

eq := fDOG(x, S__a, S__d) = half_hpeak assuming x > X0;

MyChange := ln(S__a/S__d) = (u/sqrt(2))^2*(S__a-S__d); # u is > 0
eval(eq, MyChange);
eval(%, [S__a=2, S__d=1]);
 

-exp(-(1/2)*S__a*x^2)+exp(-(1/2)*S__d*x^2) = -(1/2)*exp(-S__a*ln(S__a/S__d)/(S__a-S__d))+(1/2)*exp(-S__d*ln(S__a/S__d)/(S__a-S__d))

  

ln(S__a/S__d) = (1/2)*u^2*(S__a-S__d)

  

-exp(-(1/2)*S__a*x^2)+exp(-(1/2)*S__d*x^2) = -(1/2)*exp(-(1/2)*S__a*u^2)+(1/2)*exp(-(1/2)*S__d*u^2)

  

-exp(-x^2)+exp(-(1/2)*x^2) = -(1/2)*exp(-u^2)+(1/2)*exp(-(1/2)*u^2)

 (7)
> 

xofu := [solve(%, x)];
 

[(-2*ln(1/2+(1/2)*(1+2*(exp(-(1/2)*u^2))^2-2*exp(-(1/2)*u^2))^(1/2)))^(1/2), -(-2*ln(1/2+(1/2)*(1+2*(exp(-(1/2)*u^2))^2-2*exp(-(1/2)*u^2))^(1/2)))^(1/2), (-2*ln(1/2-(1/2)*(1+2*(exp(-(1/2)*u^2))^2-2*exp(-(1/2)*u^2))^(1/2)))^(1/2), -(-2*ln(1/2-(1/2)*(1+2*(exp(-(1/2)*u^2))^2-2*exp(-(1/2)*u^2))^(1/2)))^(1/2)]

  

2^(1/2)*((S__a-S__d)*ln(S__a/S__d))^(1/2)/(S__a-S__d), -2^(1/2)*((S__a-S__d)*ln(S__a/S__d))^(1/2)/(S__a-S__d)

 (8)
> 

solve(MyChange, u);
x_half_peak := eval(xofu, u=%[1]);

2^(1/2)*((S__a-S__d)*ln(S__a/S__d))^(1/2)/(S__a-S__d), -2^(1/2)*((S__a-S__d)*ln(S__a/S__d))^(1/2)/(S__a-S__d)

  

[(-2*ln(1/2+(1/2)*(1+2*(exp(-ln(S__a/S__d)/(S__a-S__d)))^2-2*exp(-ln(S__a/S__d)/(S__a-S__d)))^(1/2)))^(1/2), -(-2*ln(1/2+(1/2)*(1+2*(exp(-ln(S__a/S__d)/(S__a-S__d)))^2-2*exp(-ln(S__a/S__d)/(S__a-S__d)))^(1/2)))^(1/2), (-2*ln(1/2-(1/2)*(1+2*(exp(-ln(S__a/S__d)/(S__a-S__d)))^2-2*exp(-ln(S__a/S__d)/(S__a-S__d)))^(1/2)))^(1/2), -(-2*ln(1/2-(1/2)*(1+2*(exp(-ln(S__a/S__d)/(S__a-S__d)))^2-2*exp(-ln(S__a/S__d)/(S__a-S__d)))^(1/2)))^(1/2)]

 (9)
> 

{abs~(x_half_peak)[]}

{2^(1/2)*abs(ln(1/2-(1/2)*(1+2*(exp(-ln(S__a/S__d)/(S__a-S__d)))^2-2*exp(-ln(S__a/S__d)/(S__a-S__d)))^(1/2)))^(1/2), 2^(1/2)*abs(ln(1/2+(1/2)*(1+2*(exp(-ln(S__a/S__d)/(S__a-S__d)))^2-2*exp(-ln(S__a/S__d)/(S__a-S__d)))^(1/2)))^(1/2)}

 (10)
> 

sa := 2:
sd := 1:

eval(x_half_peak, [S__a=sa, S__d=sd]):
select(is@`>`, %, 0);
evalf(%);

[(-2*ln(1/2+(1/4)*2^(1/2)))^(1/2), (-2*ln(1/2-(1/4)*2^(1/2)))^(1/2)]

  

[.5627560464, 1.960150176]

 (11)
> 

# Applying "select" directly on x_half_peak doesn't work.
# I believe the following is correct but this must be checked

`correct?` := {abs~(x_half_peak)[]};
eval(`correct?`, [S__a=sa, S__d=sd]):
evalf(%);

{2^(1/2)*abs(ln(1/2-(1/2)*(1+2*(exp(-ln(S__a/S__d)/(S__a-S__d)))^2-2*exp(-ln(S__a/S__d)/(S__a-S__d)))^(1/2)))^(1/2), 2^(1/2)*abs(ln(1/2+(1/2)*(1+2*(exp(-ln(S__a/S__d)/(S__a-S__d)))^2-2*exp(-ln(S__a/S__d)/(S__a-S__d)))^(1/2)))^(1/2)}

  

{.5627560463, 1.960150176}

 (12)
> 

 

Download Formal__Half_Width_mmcdara.mw

This has been made possible after introducing a change of variables

u = sqrt(2)*sqrt((S__a-S__d)*ln(S__a/S__d))/(S__a-S__d)

Yet another way...

mmcdara 7891
October 09 2022
1 0


 

> 

restart:
> 

H_lines := seq(CURVES([[0, i], [1, i]], LINESTYLE(2)), i in [seq](0..1, 0.1)):
V_lines := seq(CURVES([[i, 0], [i, 1]], LINESTYLE(2)), i in [seq](0..1, 0.1)):
H_vec   := plottools:-arrow([0, 0], [0.1, 0], .005, .02, .4, color=black):
V_vec   := plottools:-arrow([0, 0], [0, 0.1], .005, .02, .4, color=black):
> 

p := PLOT(H_lines, V_lines, H_vec, V_vec):
plots:-display(p, axes=none);

 
> 

T := (a, b) -> plottools:-transform((x, y) -> [x+a*y, y+b*x]):

plots:-display(T(2, 3)(p), axes=none)

 
> 

 


 

Download One_method.mw

Agree:...

mmcdara 7891
October 08 2022
1 1

for some unknown reason Maple (2015... what is the version you use?) is wrong when compared ro reference Statistical Softwares

Result obtained with R:

shapiro.test(S)

	Shapiro-Wilk normality test

data:  S
W = 0.98259, p-value = 0.3481

for fun: SWWT.mw

For your information:
I had already mentionned years ago some problemes with Statistics[ChiSquareGoodnessOfFitTest]

As you seem familiar with dedicated statistical softwares, I can't help but recommend you to use them to verify the results provided by Maple, before using this latter with full confidence.

It is highly likely that the formal inte...

mmcdara 7891
October 01 2022
0 1

It is highly likely that the formal integration cannot be done.

Ir remains you the numerical integration, which requires setting the values of p, r and "i" (the imaginary unit?).


 

> 

restart:
> 

local exp;  # In Maple "exp" stands for "exponential",
            # thus you must say explicitely that "exp" is a local name

Is "i" the imaginary unit?
If it is so you must use "I"instead

> 

exp := -sin(alpha)*i*r*(-sin(alpha)*cos(phi)*cos(theta)+sin(theta)*cos(alpha))/(4*sqrt(-2*sin(theta)*sin(alpha)*cos(phi)*p*r-2*cos(alpha)*cos(theta)*p*r+p^2+r^2)*Pi*(-2*sin(theta)*sin(alpha)*cos(phi)*p*r-2*cos(alpha)*cos(theta)*p*r+p^2+r^2)*(-2+sqrt(2))*Pi)

-(1/4)*sin(alpha)*i*r*(-sin(alpha)*cos(phi)*cos(theta)+sin(theta)*cos(alpha))/((-2*sin(theta)*sin(alpha)*cos(phi)*p*r-2*cos(alpha)*cos(theta)*p*r+p^2+r^2)^(3/2)*Pi^2*(-2+2^(1/2)))

 (1)

You do not need any assumptions for they are implicitely defined by the integration bounds.
The conditions about theta can or cannot be usefull, we will see this later

> 

Int(exp*p^2*sin(alpha), p = 0 .. 1, alpha = 0 .. (1/4)*Pi, phi = 0 .. 2*Pi)

Int(-(1/4)*sin(alpha)^2*i*r*(-sin(alpha)*cos(phi)*cos(theta)+sin(theta)*cos(alpha))*p^2/((-2*sin(theta)*sin(alpha)*cos(phi)*p*r-2*cos(alpha)*cos(theta)*p*r+p^2+r^2)^(3/2)*Pi^2*(-2+2^(1/2))), p = 0 .. 1, alpha = 0 .. (1/4)*Pi, phi = 0 .. 2*Pi)

 (2)

Let us to compute this "simple integral

> 

I_p := Int(exp*p^2*sin(alpha), p = 0 .. 1)

Int(-(1/4)*sin(alpha)^2*i*r*(-sin(alpha)*cos(phi)*cos(theta)+sin(theta)*cos(alpha))*p^2/((-2*sin(theta)*sin(alpha)*cos(phi)*p*r-2*cos(alpha)*cos(theta)*p*r+p^2+r^2)^(3/2)*Pi^2*(-2+2^(1/2))), p = 0 .. 1)

 (3)
> 

# The integrand is of the form

b := add(B[i]*p^i, i=0..2):
c := add(C[i]*p^i, i=0..2):

J_p := A*p^2 / ( sqrt(b)*c )

A*p^2/((p^2*B[2]+p*B[1]+B[0])^(1/2)*(p^2*C[2]+p*C[1]+C[0]))

 (4)
> 

# The formal integration alreaqy fayls on this even simpler expression

infolevel[int] := 4:
p^2 / sqrt(b);
int(%, p=0..1)

p^2/(p^2*B[2]+p*B[1]+B[0])^(1/2)

  

Definite Integration:   Integrating expression on p=0..1
Definite Integration:   Using the integrators [distribution, piecewise, series, o, polynomial, ln, lookup, cook, ratpoly, elliptic, elliptictrig, meijergspecial, improper, asymptotic, ftoc, ftocms, meijerg, contour]
Definite Integration:   Trying method distribution.
Definite Integration:   Trying method piecewise.
Definite Integration:   Trying method series.
Definite Integration:   Trying method o.
Definite Integration:   Trying method polynomial.
Definite Integration:   Trying method ln.
Definite Integration:   Trying method lookup.
LookUp Integrator:   unable to find the specified integral in the table
Definite Integration:   Trying method cook.
Definite Integration:   Trying method ratpoly.
Definite Integration:   Trying method elliptic.
int/elliptic: trying elliptic integration
Definite Integration:   Trying method elliptictrig.
Definite Integration:   Trying method meijergspecial.
Definite Integration:   Trying method improper.
Definite Integration:   Trying method asymptotic.
Definite Integration:   Trying method ftoc.

Warning,  computation interrupted

  

Attempt to do numerical integration

> 

infolevel[int] := 0:
> 

I_theta := proc(__theta, __r, __i,  meth)
  eval(exp*p^2*sin(alpha), [theta=__theta, r=__r, i=1]):
  __i*evalf(Int(%, p=0..1, alpha=0..(1/4)*Pi, phi=0..2*Pi, method = meth))
end proc:
> 

I_theta(0, 1, I, _d01ajc);     # fails

I*(Int(Int(Int(-0.4324151988e-1*sin(alpha)^3*cos(phi)*p^2/(1.-2.*cos(alpha)*p+p^2)^(3/2), p = 0. .. 1.), alpha = 0. .. .7853981635), phi = 0. .. 6.283185308))

 (5)
> 

I_theta(0, 1, I, _CubaVegas);  # interrupted

Warning,  computation interrupted

  
> 

I_theta(0, 1, I, _MonteCarlo); # fails
 

I*(Int(Int(Int(-0.4324151988e-1*sin(alpha)^3*cos(phi)*p^2/(1.-2.*cos(alpha)*p+p^2)^(3/2), p = 0. .. 1.), alpha = 0. .. .7853981635), phi = 0. .. 6.283185308))

 (6)

Monte-Carlo integration can be done "by hand":

> 

with(Statistics):
S_N     := 10^3:   # Change this value to increase the precision
                   # A good practice is to repeat the integration with different input
                   # samples of the same size and assess the integration error.
                   # From a statistical point of view it can be better to do K Monte-Carlo
                   # integrations with S_N points that a single one with K*S_N points.
                   # The Bootstrap method (see help pages) can be use then to refine
                   # the estimation of the integral.

S_p     := Sample(Uniform(0, 1)   , S_N):
S_alpha := Sample(Uniform(0, Pi/4), S_N):
S_phi   := Sample(Uniform(0, 2*Pi), S_N):

S_I := proc(__theta, __i, __r, M)
  local f:
  f := unapply(eval(exp*p^2*sin(alpha), [theta=__theta, r=__r, i=1]), [p, alpha, phi]):
  __i * add(f(S_p[n], S_alpha[n], S_phi[n]), n=1..M) / M:
  return evalf(simplify(%));
end proc:
  
> 

# Example

StringTools:-FormatTime("%H:%M:%S");
DocumentTools:-Tabulate(
  [
    plot3d('Re(S_I(P, R, I, 100))', P=0..1, R=0..1, labels=["p", "r", "Int"], grid=[20, 20], title="Real part"),
    plot3d('Im(S_I(P, R, I, 100))', P=0..1, R=0..1, labels=["p", "r", "Int"], grid=[20, 20], title="Imaginary part")
  ],
  width=60
):
StringTools:-FormatTime("%H%M%S");

"10:00:11"

  

"100014"

 (7)
> 

`I'm not in a hurry` := false:

if `I'm not in a hurry` then
  DocumentTools:-Tabulate(
    [
      plot3d('Re(S_I(P, R, I, S__N))', P=0..1, R=0..1, labels=["p", "r", "Int"], title="Real part"),
      plot3d('Im(S_I(P, R, I, S__N))', P=0..1, R=0..1, labels=["p", "r", "Int"], title="Imaginary part")
    ],
    width=60
  ):
end if:
> 

# Sketch of the intergation error
# (just to give you an idea of the variability of the Monte Carlo intagration)



S_I2 := proc(__theta, __i, __r, M1, M2)
  local f:
  f := unapply(eval(exp*p^2*sin(alpha), [theta=__theta, r=__r, i=1]), [p, alpha, phi]):
  __i * add(f(S_p[n], S_alpha[n], S_phi[n]), n=M1..M2) / (M2-M1+1):
  return evalf(simplify(%));
end proc:

for k from 1 to 10 do
  printf("Integration over block %2d = %+1.3e  %+1.3e * I\n", k, (Re, Im)(S_I2(0, 1, I, 1+(k-1)*100, k*100)));
end do;

Integration over block  1 = -7.288e-05  -1.275e-04 * I
Integration over block  2 = -4.107e-05  -7.174e-05 * I
Integration over block  3 = -1.142e-06  +8.639e-06 * I
Integration over block  4 = -1.725e-06  +1.745e-05 * I
Integration over block  5 = +1.180e-04  +7.142e-05 * I
Integration over block  6 = +8.563e-05  +2.254e-04 * I
Integration over block  7 = -1.454e-06  -1.150e-04 * I
Integration over block  8 = +6.503e-05  +1.921e-04 * I
Integration over block  9 = -7.424e-05  -2.054e-04 * I
Integration over block 10 = +2.240e-05  -3.212e-05 * I

  
> 

V := Vector(S_N, n -> evalf(eval(exp*p^2*sin(alpha), [p=S_p[n], alpha=S_alpha[n], phi=S_phi[n], theta=0, r=1]))):
Mean(eval(V, i=1));
Boot := Bootstrap('Mean', eval(V, i=1), replications=1000, output=['value', 'standarderror']);

HFloat(-1.9365481621754508e-4)

  

[HFloat(-1.824804060830804e-4), 0.352362351056818162e-3]

 (8)
> 

# The Bootstrap estimation of the integral is val = -0.000186623004800676
# and its error is err = 0.000350218565372046231
#
# A 95% bilateral confidence interval is

Boot[1]+Boot[2]*Quantile(Normal(0, 1), 0.025) .. Boot[1]+Boot[2]*Quantile(Normal(0, 1), 0.975)

HFloat(-8.730979236620875e-4) .. HFloat(5.081371114959267e-4)

 (9)
> 

 


 

Download Int_mmcdara.mw

 

A possibility >...

mmcdara 7891
September 28 2022
1 4

A possibility
(it seems that all dthe A[i] but A[2] are null, does this seem right to you?)

> 

restart
> 

with(LinearAlgebra):
> 

with(PDEtools):
> 

Setup(mathematicalnotation = true);

Setup(mathematicalnotation = true)

 (1)
> 

assume(x::real);
> 

assume(t::real);
> 

assume(xi::real);
> 

assume(tau::real);
> 

interface(showassumed = 0);

0

 (2)
> 

alias(u = u(x, t), q = q(x, t));

u, q

 (3)
> 

Eq := diff(u, `$`(x, 2))-u^2-a*u-b*x-c*x^2 = 0;

diff(diff(u, x), x)-u^2-a*u-b*x-c*x^2 = 0

 (4)
> 

va := {u = A[0]/(x-x[0])^p};

{u = A[0]/(x-x[0])^p}

  

{diff(diff(u, x), x) = A[0]*(x-x[0])^(-p-2)*p*(p+1)}

  

A[0]*(x-x[0])^(-p-2)*p*(p+1)

 (5)
> 

Eq1 := simplify(eval(Eq, va), size);

A[0]*p^2/((x-x[0])^p*(x-x[0])^2)+A[0]*p/((x-x[0])^p*(x-x[0])^2)-A[0]^2/((x-x[0])^p)^2-a*A[0]/(x-x[0])^p-b*x-c*x^2 = 0

 (6)
> 

va1 := {x = x[0]+eta, x-x[0] = eta};

{x = x[0]+eta, x-x[0] = eta}

 (7)
> 

Eq2 := simplify(eta^(p+2)*simplify(eval(Eq1, va1)));

A[0]*p^2+A[0]*p-eta^(-p+2)*A[0]^2-eta^2*a*A[0]-eta^(p+3)*b-eta^(p+2)*b*x[0]-eta^(p+4)*c-2*eta^(p+3)*c*x[0]-eta^(p+2)*c*x[0]^2 = 0

 (8)
> 

``
> 

simplify(eval(Eq2, [p = 2]), size);

-c*eta^6-2*c*eta^5*x[0]-c*eta^4*x[0]^2-b*eta^5-b*eta^4*x[0]-a*eta^2*A[0]-A[0]^2+6*A[0] = 0

 (9)
> 

> 

u = expand(sum(A[m]*eta^m, m = 0 .. infinity))/eta^2;

u = (sum(A[m]*eta^m, m = 0 .. infinity))/eta^2

 (10)
> 

> 

Eq3 := subs(u = expand(sum(A[m]*eta^m, m = 0 .. infinity))/eta^2, Eq);

diff(diff((sum(A[m]*eta^m, m = 0 .. infinity))/eta^2, x), x)-(sum(A[m]*eta^m, m = 0 .. infinity))^2/eta^4-a*(sum(A[m]*eta^m, m = 0 .. infinity))/eta^2-b*x-c*x^2 = 0

 (11)
> 

GenSys := proc(N)
   local expr         := collect(numer(normal(eval(lhs(Eq3), infinity=N))), eta);
   local coefficients := [coeffs(expr, eta, 'k')];
   local unknowns     := [select(has, indets(coefficients, name), A)[]];
   local NbC          := numelems(coefficients);
  
   coefficients := select(has, coefficients=~[0$NbC], A):
   NbC          := numelems(coefficients);

   return [coefficients, unknowns]
end proc:

deg := 1:
g   := GenSys(deg);
sol := solve(g[1], g[2]);

[[-A[0]^2 = 0, -a*A[1] = 0, -a*A[0]-A[1]^2 = 0, -2*A[0]*A[1] = 0], [A[0], A[1]]]

  

[[A[0] = 0, A[1] = 0]]

 (12)
> 

deg := 2:
g   := GenSys(deg);
sol := solve(g[1], g[2]);

[[-A[0]^2 = 0, -c*x^2-a*A[2]-b*x-A[2]^2 = 0, -a*A[1]-2*A[1]*A[2] = 0, -a*A[0]-2*A[0]*A[2]-A[1]^2 = 0, -2*A[0]*A[1] = 0], [A[0], A[1], A[2]]]

  

[[A[0] = 0, A[1] = 0, A[2] = RootOf(c*x^2+_Z^2+_Z*a+b*x)]]

 (13)
> 

vals := allvalues(eval(A[2], sol[])):
map(u -> [remove(has, sol[], A[2])[], A[2]=u], [vals]):
print~(%):

[A[0] = 0, A[1] = 0, A[2] = -(1/2)*a+(1/2)*(-4*c*x^2+a^2-4*b*x)^(1/2)]

  

[A[0] = 0, A[1] = 0, A[2] = -(1/2)*a-(1/2)*(-4*c*x^2+a^2-4*b*x)^(1/2)]

 (14)
> 

 

NULL

 

Download PA_mmcdara.mw

Improvement (automatic treatment of equations of the form A[m]=RootOf(...))

deg := 10: # for instance
g   := GenSys(deg):
sol := solve(g[1], g[2]):

NonZero   := map(u -> if rhs(u) <> 0 then u end if, sol[]);

if NonZero <> [] then 
  NonZero_u := lhs~(NonZero):
  NonZero_v := [allvalues~(rhs~(NonZero))]:
  map(
    i -> 
      seq(
        [
          remove(has, sol[], NonZero_u)[], 
          NonZero_u[i]=NonZero_v[i][j]
        ], 
        j=1..numelems(NonZero_v[i])
      ),
      [$1..numelems(NonZero_u)]
  ):
  print~(%):
else
  print(sol[])
end if:

 

Looking at the integrand it appears it i...

mmcdara 7891
September 25 2022
2 0


Looking at the integrand it appears it is almost constant wrt y and oscillating wrt x.
Using evalf(Int(integrand, x=..., y=..., method=...)) seems to be challenging for no integration method seems perfecly suited to handle this situation.
Here is a proposal (all depends if the precision you need)

(PS: I use CurveFitting:-Spline only because I'm a lazzy person)

> 

restart:
> 

afa:=0.3:
vh:=3.5:
u:=3.12:
mu:=5.5:
gama:=-4*10^(-29)*(1-cos(6*afa))*(1-1*10^(-8)*I):
d1:=1.78*10^(-9):
d2:=48.22*10^(-9):
HBAR:=1.05457266*10^(-34):
ME:=9.1093897*10^(-31):
ELEC:=1.60217733*10^(-19):
Kh:=2.95*10^10:
kc:=sqrt(2*ME*ELEC/HBAR^2):
k:=kc*sqrt(mu):
k0:=sqrt(k^2-k^2*sin(x)^2):
kh:=sqrt(k^2-(Kh-k*sin(x)*cos(y))^2 - k^2 * sin(x)^2 * sin(y)^2):
khg:=sqrt(k^2-(2*Kh*sin(afa/2)*sin(afa/2)-k*sin(x)*cos(y))^2-(2*Kh*sin(afa/2)*cos(afa/2)+k*sin(x)*sin(y))^2):
kg1:=sqrt(k^2-(Kh*cos(Pi/3-afa)-k*sin(x)*cos(y))^2-(Kh*sin(Pi/3-afa)+k*sin(x)*sin(y))^2):
kg2:=sqrt(k^2-(Kh*cos(afa)-k*sin(x)*cos(y))^2-(k*sin(x)*sin(y)-Kh*sin(afa))^2):
k0pl:=sqrt(k^2-k^2*sin(x)^2+kc^2*vh-kc^2*u):
k0mi:=sqrt(k^2-k^2*sin(x)^2-kc^2*vh-kc^2*u):
khpl:=sqrt(k^2-(Kh-k*sin(x)*cos(y))^2-k^2*sin(x)^2*sin(y)^2+kc^2*vh-kc^2*u):
khmi:=sqrt(k^2-(Kh-k*sin(x)*cos(y))^2-k^2*sin(x)^2*sin(y)^2-kc^2*vh-kc^2*u):
k0plpl:=sqrt(k^2-k^2*sin(x)^2+2*kc^2*vh):
k0mimi:=sqrt(k^2-k^2*sin(x)^2-2*kc^2*vh):
khplpl:=sqrt(k^2-(Kh-k*sin(x)*cos(y))^2-k^2*sin(x)^2*sin(y)^2+2*kc^2*vh):
khmimi:=sqrt(k^2-(Kh-k*sin(x)*cos(y))^2-k^2*sin(x)^2*sin(y)^2-2*kc^2*vh):
khgplpl:=sqrt(k^2-(2*Kh*sin(afa/2)*sin(afa/2)-k*sin(x)*cos(y))^2-(2*Kh*sin(afa/2)*cos(afa/2)+k*sin(x)*sin(y))^2+2*kc^2*vh):
khgmimi:=sqrt(k^2-(2*Kh*sin(afa/2)*sin(afa/2)-k*sin(x)*cos(y))^2-(2*Kh*sin(afa/2)*cos(afa/2)+k*sin(x)*sin(y))^2-2*kc^2*vh):
kg1plpl:=sqrt(k^2-(Kh*cos(Pi/3-afa)-k*sin(x)*cos(y))^2-(Kh*sin(Pi/3-afa)+k*sin(x)*sin(y))^2+2*kc^2*vh):
kg1mimi:=sqrt(k^2-(Kh*cos(Pi/3-afa)-k*sin(x)*cos(y))^2-(Kh*sin(Pi/3-afa)+k*sin(x)*sin(y))^2-2*kc^2*vh):
kg2plpl:=sqrt(k^2-(Kh*cos(afa)-k*sin(x)*cos(y))^2-(k*sin(x)*sin(y)-Kh*sin(afa))^2+2*kc^2*vh):
kg2mimi:=sqrt(k^2-(Kh*cos(afa)-k*sin(x)*cos(y))^2-(k*sin(x)*sin(y)-Kh*sin(afa))^2-2*kc^2*vh):
A1:=1/(1+I*ME*gama/(HBAR^2*k0pl))*exp(I*k0pl*d1)/2:
B1:=1/(1+I*ME*gama/(HBAR^2*k0pl))*exp(I*k0pl*d1)/2:
A2:=1/(1+I*ME*gama/(HBAR^2*khpl))*exp(I*khpl*d1)/2:
B2:=1/(1+I*ME*gama/(HBAR^2*khpl))*exp(I*khpl*d1)/2:
A3:=1/(1+I*ME*gama/(HBAR^2*k0mi))*exp(I*k0mi*d1)/2:
B3:=1/(1+I*ME*gama/(HBAR^2*k0mi))*exp(I*k0mi*d1)/2:
A4:=1/(1+I*ME*gama/(HBAR^2*khmi))*exp(I*khmi*d1)/2:
B4:=1/(1+I*ME*gama/(HBAR^2*khmi))*exp(I*khmi*d1)/2:
> 

T1:=1/4*Re(abs(A1)^2*k0plpl*exp(I*(k0plpl-conjugate(k0plpl))*d2)+abs(A1)^2*kg1plpl*exp(I*(kg1plpl-conjugate(kg1plpl))*d2)+abs(A2)^2*khplpl*exp(I*(khplpl-conjugate(khplpl))*d2)+abs(A2)^2*khgplpl*exp(I*(khgplpl-conjugate(khgplpl))*d2)+abs(B3)^2*k0mimi*exp(I*(k0mimi-conjugate(k0mimi))*d2)+abs(B3)^2*kg1mimi*exp(I*(kg1mimi-conjugate(kg1mimi))*d2)+abs(B4)^2*khmimi*exp(I*(khmimi-conjugate(khmimi))*d2)+abs(B4)^2*khgmimi*exp(I*(khgmimi-conjugate(khgmimi))*d2)+abs(B1+A3)^2*k0*exp(I*(k0-conjugate(k0))*d2)+abs(B1-A3)^2*kg1*exp(I*(kg1-conjugate(kg1))*d2)+abs(A4-B2)^2*kh*exp(I*(kh-conjugate(kh))*d2)+abs(A4+B2)^2*khg*exp(I*(khg-conjugate(khg))*d2)+conjugate(A1)*B3*k0mimi*exp(I*(k0mimi-conjugate(k0plpl))*d2)+A1*conjugate(B3)*k0plpl*exp(I*(k0plpl-conjugate(k0mimi))*d2)+conjugate(A1)*(B1+A3)*k0*exp(I*(k0-conjugate(k0plpl))*d2)+A1*conjugate(B1+A3)*k0plpl*exp(I*(k0plpl-conjugate(k0))*d2)+conjugate(B3)*(B1+A3)*k0*exp(I*(k0-conjugate(k0mimi))*d2)+B3*conjugate(B1+A3)*k0mimi*exp(I*(k0mimi-conjugate(k0))*d2)-conjugate(A1)*B3*kg1mimi*exp(I*(kg1mimi-conjugate(kg1plpl))*d2)-A1*conjugate(B3)*kg1plpl*exp(I*(kg1plpl-conjugate(kg1mimi))*d2)+conjugate(A1)*(A3-B1)*kg1*exp(I*(kg1-conjugate(kg1plpl))*d2)+A1*conjugate(A3-B1)*kg1plpl*exp(I*(kg1plpl-conjugate(kg1))*d2)+conjugate(B3)*(B1-A3)*kg1*exp(I*(kg1-conjugate(kg1mimi))*d2)+B3*conjugate(B1-A3)*kg1mimi*exp(I*(kg1mimi-conjugate(kg1))*d2)-conjugate(A2)*B4*khmimi*exp(I*(khmimi-conjugate(khplpl))*d2)-A2*conjugate(B4)*khplpl*exp(I*(khplpl-conjugate(khmimi))*d2)+conjugate(A2)*(B2-A4)*kh*exp(I*(kh-conjugate(khplpl))*d2)+A2*conjugate(B2-A4)*khplpl*exp(I*(khplpl-conjugate(kh))*d2)+conjugate(B4)*(A4-B2)*kh*exp(I*(kh-conjugate(khmimi))*d2)+B4*conjugate(A4-B2)*khmimi*exp(I*(khmimi-conjugate(kh))*d2)+conjugate(A2)*B4*khgmimi*exp(I*(khgmimi-conjugate(khgplpl))*d2)+A2*conjugate(B4)*khgplpl*exp(I*(khgplpl-conjugate(khgmimi))*d2)-conjugate(A2)*(A4+B2)*khg*exp(I*(khg-conjugate(khgplpl))*d2)-A2*conjugate(A4+B2)*khgplpl*exp(I*(khgplpl-conjugate(khg))*d2)-conjugate(B4)*(A4+B2)*khg*exp(I*(khg-conjugate(khgmimi))*d2)-B4*conjugate(A4+B2)*khgmimi*exp(I*(khgmimi-conjugate(khg))*d2)):
> 


T2:=1/4*Re(abs(A1)^2*k0plpl*exp(I*(k0plpl-conjugate(k0plpl))*d2)+abs(A1)^2*kg2plpl*exp(I*(kg2plpl-conjugate(kg2plpl))*d2)+abs(A2)^2*khplpl*exp(I*(khplpl-conjugate(khplpl))*d2)+abs(A2)^2*khgplpl*exp(I*(khgplpl-conjugate(khgplpl))*d2)+abs(B3)^2*k0mimi*exp(I*(k0mimi-conjugate(k0mimi))*d2)+abs(B3)^2*kg2mimi*exp(I*(kg2mimi-conjugate(kg2mimi))*d2)+abs(B4)^2*khmimi*exp(I*(khmimi-conjugate(khmimi))*d2)+abs(B4)^2*khgmimi*exp(I*(khgmimi-conjugate(khgmimi))*d2)+abs(B1+A3)^2*k0*exp(I*(k0-conjugate(k0))*d2)+abs(B1-A3)^2*kg2*exp(I*(kg2-conjugate(kg2))*d2)+abs(A4-B2)^2*kh*exp(I*(kh-conjugate(kh))*d2)+abs(A4+B2)^2*khg*exp(I*(khg-conjugate(khg))*d2)+conjugate(A1)*B3*k0mimi*exp(I*(k0mimi-conjugate(k0plpl))*d2)+A1*conjugate(B3)*k0plpl*exp(I*(k0plpl-conjugate(k0mimi))*d2)+conjugate(A1)*(B1+A3)*k0*exp(I*(k0-conjugate(k0plpl))*d2)+A1*conjugate(B1+A3)*k0plpl*exp(I*(k0plpl-conjugate(k0))*d2)+conjugate(B3)*(B1+A3)*k0*exp(I*(k0-conjugate(k0mimi))*d2)+B3*conjugate(B1+A3)*k0mimi*exp(I*(k0mimi-conjugate(k0))*d2)-conjugate(A1)*B3*kg2mimi*exp(I*(kg2mimi-conjugate(kg2plpl))*d2)-A1*conjugate(B3)*kg2plpl*exp(I*(kg2plpl-conjugate(kg2mimi))*d2)+conjugate(A1)*(A3-B1)*kg2*exp(I*(kg2-conjugate(kg2plpl))*d2)+A1*conjugate(A3-B1)*kg2plpl*exp(I*(kg2plpl-conjugate(kg2))*d2)+conjugate(B3)*(B1-A3)*kg2*exp(I*(kg2-conjugate(kg2mimi))*d2)+B3*conjugate(B1-A3)*kg2mimi*exp(I*(kg2mimi-conjugate(kg2))*d2)-conjugate(A2)*B4*khmimi*exp(I*(khmimi-conjugate(khplpl))*d2)-A2*conjugate(B4)*khplpl*exp(I*(khplpl-conjugate(khmimi))*d2)+conjugate(A2)*(B2-A4)*kh*exp(I*(kh-conjugate(khplpl))*d2)+A2*conjugate(B2-A4)*khplpl*exp(I*(khplpl-conjugate(kh))*d2)+conjugate(B4)*(A4-B2)*kh*exp(I*(kh-conjugate(khmimi))*d2)+B4*conjugate(A4-B2)*khmimi*exp(I*(khmimi-conjugate(kh))*d2)+conjugate(A2)*B4*khgmimi*exp(I*(khgmimi-conjugate(khgplpl))*d2)+A2*conjugate(B4)*khgplpl*exp(I*(khgplpl-conjugate(khgmimi))*d2)-conjugate(A2)*(A4+B2)*khg*exp(I*(khg-conjugate(khgplpl))*d2)-A2*conjugate(A4+B2)*khgplpl*exp(I*(khgplpl-conjugate(khg))*d2)-conjugate(B4)*(A4+B2)*khg*exp(I*(khg-conjugate(khgmimi))*d2)-B4*conjugate(A4+B2)*khgmimi*exp(I*(khgmimi-conjugate(khg))*d2)):
> 

#evalf(Int(k*sin(x)*T1,x=0..Pi/2,y=-Pi/6..-Pi/6+afa))
> 

indets(T1, name)

{x, y}

 (1)
> 

T1:=1/4*(abs(A1)^2*k0plpl*exp(I*(k0plpl-conjugate(k0plpl))*d2)+abs(A1)^2*kg1plpl*exp(I*(kg1plpl-conjugate(kg1plpl))*d2)+abs(A2)^2*khplpl*exp(I*(khplpl-conjugate(khplpl))*d2)+abs(A2)^2*khgplpl*exp(I*(khgplpl-conjugate(khgplpl))*d2)+abs(B3)^2*k0mimi*exp(I*(k0mimi-conjugate(k0mimi))*d2)+abs(B3)^2*kg1mimi*exp(I*(kg1mimi-conjugate(kg1mimi))*d2)+abs(B4)^2*khmimi*exp(I*(khmimi-conjugate(khmimi))*d2)+abs(B4)^2*khgmimi*exp(I*(khgmimi-conjugate(khgmimi))*d2)+abs(B1+A3)^2*k0*exp(I*(k0-conjugate(k0))*d2)+abs(B1-A3)^2*kg1*exp(I*(kg1-conjugate(kg1))*d2)+abs(A4-B2)^2*kh*exp(I*(kh-conjugate(kh))*d2)+abs(A4+B2)^2*khg*exp(I*(khg-conjugate(khg))*d2)+conjugate(A1)*B3*k0mimi*exp(I*(k0mimi-conjugate(k0plpl))*d2)+A1*conjugate(B3)*k0plpl*exp(I*(k0plpl-conjugate(k0mimi))*d2)+conjugate(A1)*(B1+A3)*k0*exp(I*(k0-conjugate(k0plpl))*d2)+A1*conjugate(B1+A3)*k0plpl*exp(I*(k0plpl-conjugate(k0))*d2)+conjugate(B3)*(B1+A3)*k0*exp(I*(k0-conjugate(k0mimi))*d2)+B3*conjugate(B1+A3)*k0mimi*exp(I*(k0mimi-conjugate(k0))*d2)-conjugate(A1)*B3*kg1mimi*exp(I*(kg1mimi-conjugate(kg1plpl))*d2)-A1*conjugate(B3)*kg1plpl*exp(I*(kg1plpl-conjugate(kg1mimi))*d2)+conjugate(A1)*(A3-B1)*kg1*exp(I*(kg1-conjugate(kg1plpl))*d2)+A1*conjugate(A3-B1)*kg1plpl*exp(I*(kg1plpl-conjugate(kg1))*d2)+conjugate(B3)*(B1-A3)*kg1*exp(I*(kg1-conjugate(kg1mimi))*d2)+B3*conjugate(B1-A3)*kg1mimi*exp(I*(kg1mimi-conjugate(kg1))*d2)-conjugate(A2)*B4*khmimi*exp(I*(khmimi-conjugate(khplpl))*d2)-A2*conjugate(B4)*khplpl*exp(I*(khplpl-conjugate(khmimi))*d2)+conjugate(A2)*(B2-A4)*kh*exp(I*(kh-conjugate(khplpl))*d2)+A2*conjugate(B2-A4)*khplpl*exp(I*(khplpl-conjugate(kh))*d2)+conjugate(B4)*(A4-B2)*kh*exp(I*(kh-conjugate(khmimi))*d2)+B4*conjugate(A4-B2)*khmimi*exp(I*(khmimi-conjugate(kh))*d2)+conjugate(A2)*B4*khgmimi*exp(I*(khgmimi-conjugate(khgplpl))*d2)+A2*conjugate(B4)*khgplpl*exp(I*(khgplpl-conjugate(khgmimi))*d2)-conjugate(A2)*(A4+B2)*khg*exp(I*(khg-conjugate(khgplpl))*d2)-A2*conjugate(A4+B2)*khgplpl*exp(I*(khgplpl-conjugate(khg))*d2)-conjugate(B4)*(A4+B2)*khg*exp(I*(khg-conjugate(khgmimi))*d2)-B4*conjugate(A4+B2)*khgmimi*exp(I*(khgmimi-conjugate(khg))*d2)):
> 

indets(k*sin(x)*T1, name)

{x, y}

 (2)
> 

# evalf(Int(k*sin(x)*T1,x=0..Pi/2,y=-Pi/6..-Pi/6+afa))

Warning,  computation interrupted

  
> 

# the integrand seems almost invariant regarding y:

plot3d(Re(k*sin(x)*T1),x=0..Pi/2,y=-Pi/6..-Pi/6+afa, labels=[x, y, z])

 

 

method = _d01akc

--

for finite interval of integration, oscillating integrands; uses adaptive Gauss 30-point and Kronrod 61-point rules.
 

 

 

> 

# as the integrand is x-oscillating and y-constant (or almost), a suggestion is to do this

i1 := evalf(Int(eval(k*sin(x)*Re(T1), y=-Pi/6. ), x=0..Pi/2, method=_d01akc));
i2 := evalf(Int(eval(k*sin(x)*Re(T1), y=-Pi/12.), x=0..Pi/2, method=_d01akc));
i3 := evalf(Int(eval(k*sin(x)*Re(T1), y=     0.), x=0..Pi/2, method=_d01akc));
i4 := evalf(Int(eval(k*sin(x)*Re(T1), y=+Pi/12.), x=0..Pi/2, method=_d01akc));
i5 := evalf(Int(eval(k*sin(x)*Re(T1), y=+Pi/6. ), x=0..Pi/2, method=_d01akc));

0.1179159365e20

  

0.1171354144e20

  

0.1171354144e20

  

0.1171354144e20

  

0.1171354144e20

 (3)
> 

# here is an evaluation of the integral

CurveFitting:-Spline(Pi*~[seq(-1/6..1/6, 1/12)], [seq(i||k, k=1..5)], t, degree=1):
approx_1 := int(%, t=-Pi/6..Pi/6);

0.1227660893e20

 (4)
> 

# improvements: T1 depends on y only in a thin layer -Pi/6 .. -Pi/6.*(0.94)
evalf(Int(eval(k*sin(x)*Re(T1), y=-Pi/6.*(0.94)), x=0..Pi/2, method=_d01akc));

0.1171354144e20

 (5)
> 

# thus this sevond approximation of the integral

CurveFitting:-Spline([-Pi/6, -Pi/6*0.94, Pi/6], [i1, i2, i2], t, degree=1):
approx_2 := int(%, t=-Pi/6..Pi/6);

0.1226761795e20

 (6)
> 

 


 

Download text_program_mmcdara.mw

For what it worths.......

mmcdara 7891
September 23 2022
1 9

 

As dualxisplot can't do the job (as Carl said it is not intended to this), maybe yu will be interested in that:

KindOf.mw

First 32 33 34 35 36 37 38 Last Page 34 of 65