acer

32333 Reputation

29 Badges

19 years, 323 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

I don't see why it's necessary to call expand here.

This produces a result close to your supplied target expression.

By the way, it looks like you gave a target expression with a wrong sign on one L1/n12 subterm.

restart;

eqn := v1(t) - R1*(-i3(t)/n13 - i2(t)/n12) - L1*diff(-i3(t)/n13
       - i2(t)/n12, t) = n12*(v2(t) - R2*i2(t) - L2*diff(i2(t), t));

v1(t)-R1*(-i3(t)/n13-i2(t)/n12)-L1*(-(diff(i3(t), t))/n13-(diff(i2(t), t))/n12) = n12*(v2(t)-R2*i2(t)-L2*(diff(i2(t), t)))

((L,R)->L=-R)(selectremove(has,collect((rhs-lhs)(eqn),diff),diff));

(-n12*L2-L1/n12)*(diff(i2(t), t))-L1*(diff(i3(t), t))/n13 = -n12*(-R2*i2(t)+v2(t))+v1(t)-R1*(-i3(t)/n13-i2(t)/n12)

(proc(L,R)
   local s; s:=sign(L);
   collect(L*s,diff)=-R*s;
 end proc)(selectremove(has,collect((rhs-lhs)(eqn),diff),diff));

(n12*L2+L1/n12)*(diff(i2(t), t))+L1*(diff(i3(t), t))/n13 = n12*(-R2*i2(t)+v2(t))-v1(t)+R1*(-i3(t)/n13-i2(t)/n12)

simplify((lhs-rhs)(%-eqn));

0

Download collect_diff_ac.mw

You haven't told us what numeric value your "to" has.

So here are a couple of ways to do it purely as numeric quadrature, quickly, with the domain being specified by the parameter of a procedure.

restart;

f1:=proc(c) options remember, system;
  if not c::numeric then return 'procname'(_passed); end if;
  evalf(Int(exp((5+I*6)*x)*sin(1+erf(x)),x=-c..c, method=_Dexp,_rest));
end proc:

CodeTools:-Usage( f1(1.2) );

memory used=4.06MiB, alloc change=32.00MiB, cpu time=76.00ms, real time=76.00ms, gc time=7.24ms

49.09997652+1.408399186*I

CodeTools:-Usage( plot([Re(f1(t,epsilon=1e-3)),
                        Im(f1(t,epsilon=1e-3))],
                       t=0..1, gridlines=false) );

memory used=239.36MiB, alloc change=38.00MiB, cpu time=2.53s, real time=2.41s, gc time=352.02ms

restart;

#
# This is faster.
#

Ref:=evalc(Re(exp((5+I*6)*x)*sin(1+erf(x)))):

Imf:=evalc(Im(exp((5+I*6)*x)*sin(1+erf(x)))):

f2:=proc(c) options remember, system;
  if not c::numeric then return 'procname'(_passed); end if;
  evalf(Int(Ref,x=-c..c, method=_d01ajc,_rest)
        +I*Int(Imf,x=-c..c, method=_d01ajc,_rest));
end proc:

CodeTools:-Usage( f2(1.2) );

memory used=0.51MiB, alloc change=0 bytes, cpu time=9.00ms, real time=9.00ms, gc time=0ns

49.09997652+1.408399186*I

CodeTools:-Usage( plot([Re(f2(t,epsilon=1e-3)),
                        Im(f2(t,epsilon=1e-3))],
                       t=0..1, gridlines=false) );

memory used=13.54MiB, alloc change=0 bytes, cpu time=129.00ms, real time=130.00ms, gc time=0ns

 

Download evalf_int_example.mw

Don't use linalg, which is deprecated. You can extract any row of a Matrix with simple indexing.

You don't need to call print in order to display the scalar value assigned to a name.

I have used the parameter nr of aFshear to denote which row to use.

restart:

L:=Vector[row]([''RA'',''MA'',''thetaA'',''yA'',''RB'',''MB'',''thetaB'',''yB'']):

R:=Matrix(2,8,[1,2,3,4,5,6,7,8,11,12,13,14,15,16,17,18]);

Matrix(2, 8, {(1, 1) = 1, (1, 2) = 2, (1, 3) = 3, (1, 4) = 4, (1, 5) = 5, (1, 6) = 6, (1, 7) = 7, (1, 8) = 8, (2, 1) = 11, (2, 2) = 12, (2, 3) = 13, (2, 4) = 14, (2, 5) = 15, (2, 6) = 16, (2, 7) = 17, (2, 8) = 18})

aFshear:=proc(nr,R,L)
           zip(assign,L,R[nr,..]);
           NULL:
         end proc:

 

aFshear(2,R,L);

MA,thetaB;

12, 17

aFshear(1,R,L);

MA,thetaB;

2, 7

aFshear(2,R,L);

MA,thetaB;

12, 17

 

Download zip_assign.mw

You have not yet given numeric values to Nturns and Nsets.

The unapply command can make an operator from a previously computed expression (or a name to which such has been assigned).

The way you were trying to define M muddled the concepts of the formal parameters z and r of the operator with the names z and r in the expression previously computed and assigned to Bzsum. This is a common misunderstanding. In contrast, your simpler example (z,r)->z*r provides the body of the operator literally.

restart;

isol:=10;
Zres:=3:
Rres:=2:
l:=1;
d:=0.25;
mu[0]:=4*Pi*(10^(-7)):
R:=d/2:
ires:=(isol*Nturns*Nsets)/Zres:
kk:=(4*R*r)/(((R+r)^2)+((z-A)^2)):
K:=int(1/(sqrt((1-(t^2))*(1-(kk)*(t^2)))),t=0..1):
E:=int(sqrt(1-(kk)*(t^2))/(sqrt(1-(t^2))),t=0..1):
Bz:=((mu[0]*ires)/(2*Pi))*(1/(sqrt(((R+r)^2)
    +((z-A)^2))))*(K+E*((R^2)-(r^2)-((z-A)^2))/(((R-r)^2)+((z-A)^2))):
BzL:=eval~(Bz,A=~[seq((n*l)/(Zres),n=((-Zres)/2)..((Zres)/2))]):
Bzsum:=add(BzL[n],n=1..(Zres)):

10

1

.25

M:=unapply(Bzsum,[z,r]):

M1:=Array(1..Zres,i->(i*l)/(2*Zres));
M2:= Array(1..Rres,i->(i*d)/Rres);

Vector[row](3, {(1) = 1/6, (2) = 1/3, (3) = 1/2})

Array(%id = 18446884106455853766)

A=Array(1..numelems(M1),1..numelems(M2),(i,j)->M(M1[i],M2[j]));

A = (Matrix(3, 2, {(1, 1) = 0.6274841782e-6*Nturns*Nsets+(Float(undefined)-0.*I)*Nturns*Nsets, (1, 2) = -0.1152929500e-5*Nturns*Nsets, (2, 1) = 0.2069611699e-5*Nturns*Nsets, (2, 2) = 0.2547600213e-6*Nturns*Nsets, (3, 1) = 0.6580489429e-6*Nturns*Nsets, (3, 2) = 0.3184628994e-6*Nturns*Nsets}))

 

Download unapply_again.mw

As a programming excercise perhaps it's easier to start off with examples that you know have only rational roots, and then afterwards move on to more general randomly generated examples (whose roots may not all be rational).

In other words, start off by figuring out how to generate and test the set of candidates using examples for which you know a non-empty result should be found.

So, below I first deal with the case that all roots are rational. Then I move on to the case that possibly only one root is rational, to check the same technique. That seems like a reasonable beginning for you, before you try deal with the general case of randonly generated coefficients.

If you attack the general case up front then you may get muddled about whether no rational roots exist versus whether your code cannot find any properly.

restart;

randomize():

n := rand(1..6):
d := 2*rand(0..1)-1:
r := ()->d()*n():

p := expand( (r()*x-r()) * (r()*x-r()) * (r()*x-r()) );

72*x^3-24*x^2-22*x+4

Na0 := NumberTheory:-Divisors(coeff(p,x,0));

{1, 2, 4}

Nan := NumberTheory:-Divisors(coeff(p,x,degree(p,x)));

{1, 2, 3, 4, 6, 8, 9, 12, 18, 24, 36, 72}

S := {seq(seq(a/b,a=Na0),b=Nan)}:
S := S union map(`*`,S,-1):
nops(S);

36

select(r->eval(p,x=r)=0, S);

{-1/2, 1/6, 2/3}

solve(p);

1/6, -1/2, 2/3

restart;

randomize():

n := rand(1..6):
d := 2*rand(0..1)-1:
r := ()->d()*n():

p2 := expand( (r()*x-r()) * randpoly(x,degree=2,coeffs=rand(1..9),dense) );

-10*x^3-17*x^2-8*x-1

N2a0 := NumberTheory:-Divisors(coeff(p2,x,0));

{1}

N2an := NumberTheory:-Divisors(coeff(p2,x,degree(p2,x)));

{1, 2, 5, 10}

S2 := {seq(seq(a/b,a=N2a0),b=N2an)};
S2 := S2 union map(`*`,S2,-1):
nops(S2);

{1, 1/2, 1/5, 1/10}

8

ans := select(r->eval(p2,x=r)=0, S2);

{-1, -1/2, -1/5}

normal(p2/(x-ans[1]));
solve(%);

-10*x^2-7*x-1

-1/2, -1/5

 

Download ratroots.mw

I am guessing that your set may contain terms like a+b and that you would want that to be treated like the case of a constant or solitary parameter.

restart;

func := S -> not andmap(has,S,{x,y}):

func( {a*x^3-1, b*y+x, -1, c} );

                true

func( {a*x^3-1, b*y+x, c} );

                true

func( {a*x^3-1, b*y+x, -1} );

                 true

func( {a*x^3-1, b*y+x, b+a} );

                 true

func( {a*x^3-1, b*y+x} );

                 false

If you want to test only for constants or stand-alone parameters then it could be done by making the test check whether the entry s is constant or member(s,[a,b,c]) . Thus this treats the polynomial a+b differently from above. (I'm guessing that this is not what you want.)

params := [a,b,c]:

func3 := S -> ormap(s->s::constant
                       or member(s,params),
                    [op(S)]):

func( {a*x^3-1, b*y+x, -1, c} );

                true

func( {a*x^3-1, b*y+x, c} );

                true

func( {a*x^3-1, b*y+x, -1} );

                 true

func( {a*x^3-1, b*y+x} );

                 false

 

@EB1000 As of Maple 2018 the Library-generated commands appear only in the right side-panel (a.k.a. context-panel) and not in the popup.

The last item in the popup you've shown is "Open Context Panel for more...". You can click on that, or open the right side panel using the double side-chevron icon at right of the regular GUI menubars.

See the help page with Topic updates,Maple2018,ContextPanel

Did you obtain this Maple package from here? If so then you really ought to give proper attribution.

Do you have to have PHCpack installed somewhere on your machine?

Before you can run an example from within Maple do not you have to call the comand phc:-setPHCloc from the Maple package (which you reference in your Question here)? Perhaps you are supposed to call that with the location of PHCpack on you machine, perhaps at the start of each Maple session in which you intend on using this phc package.

I downloaded a precompiled 64bit Linux binary executable of PHCpack from here. I unpacked that to obtain an executable binary named phc . I put that in the same subdirectory as I use below to run the Maple script that creates the Library archive (and to which I save this Maple package).

The I ran the following in Maple 2017.2,

restart;

currentdir(cat(kernelopts(homedir),"/mapleprimes/phc")):

# Create an empty repository.
if FileTools:-Exists("phc.mla") then
  print("deleting old archive");
  FileTools:-Remove("phc.mla"):
end if;

       "deleting old archive"

LibraryTools:-Create("phc.mla"):
LibraryTools:-ShowContents( "phc.mla" );

                      []

read("phc.module"):
LibraryTools:-Save('phc', "phc.mla");
LibraryTools:-ShowContents( "phc.mla" );

    [["phc.m", [2019, 7, 2, 12, 59, 46], 41984, 1736]]

restart;
libname := cat(kernelopts(homedir),"/mapleprimes/phc"), libname:

with(phc);
    [cascade, computeResiduals, decompose, deflationStep,
     drawPaths, embed, eqnbyeqn, factor, filter,
     intersectWitnessSets, makeBertiniInput,
     makeBertiniStart, makeSolution, makeStartSystem,
     makeSystem, makeWitnessSet,
     printSolutions, printSystem, readSolutionsBertini,
     refine, removeBertiniFiles, setPHCloc,
     solutionsAppendToFile, solve, solveBertini,
     subsVariables, systemToFile, track, trackBertini]

setPHCloc(cat(kernelopts(homedir),"/mapleprimes/phc")):

T := makeSystem([x, y], [], [x^2+y^2-1, x^3+y^3-1]):

sols := solve(T):

printSolutions(T, sols);
(1) [x = -.67793e-31-.61630e-32*I, y = 1.0+.15407e-31*I]
(2) [x = -1.0+.70711*I, y = -1.0-.70711*I]
(3) [x = -.15407e-31-.11556e-31*I, y = 1.0+.13710e-31*I]
(4) [x = 1.0-.53357e-35*I, y = .79194e-30-.27271e-30*I]
(5) [x = -1.0-.70711*I, y = -1.0+.70711*I]
(6) [x = 1.0-.38519e-32*I, y = .30815e-32-.24652e-31*I]

 

The error message is coming from your attempts to invert Matrix V. But V is singular and noninvertible.

By the way you are attempting Matrix-Matrix multiplication, and not Vector-Matrix multiplication, since if the inverse of V were obtained then it would be a Matrix. V is a Matrix. There are no Vectors in your example, and VectorMatrixMultiply is not the right command.

You could also use the `.` command (just the dot), for multiplying Vectors and Matrices.

But your central problem is that V is noninvertible. There's a minor possibility that you could utilize a pseudo-inverse of V. I include it on the off chance. It may well not be what you want or expect.

restart

with(LinearAlgebra):

V := Matrix(8, 8, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = (1-theta[h])*alpha*beta*`ε`*PI/mu, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (4, 8) = alpha*beta*PI/mu, (5, 1) = 0, (5, 2) = 0, (5, 3) = gamma[o], (5, 4) = gamma[1], (5, 5) = 0, (5, 6) = 0, (5, 7) = 0, (5, 8) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 0, (6, 7) = 0, (6, 8) = 0, (7, 1) = 0, (7, 2) = 0, (7, 3) = 0, (7, 4) = 0, (7, 5) = 0, (7, 6) = 0, (7, 7) = 0, (7, 8) = 0, (8, 1) = 0, (8, 2) = 0, (8, 3) = alpha*c*lambda[b]/mu[b], (8, 4) = alpha*c*lambda[b]/mu[b], (8, 5) = 0, (8, 6) = 0, (8, 7) = 0, (8, 8) = 0}); F := Matrix(8, 8, {(1, 1) = mu, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = sigma, (1, 7) = 0, (1, 8) = 0, (2, 1) = 0, (2, 2) = mu, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = gamma[o]+mu, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = gamma[1]+delta+mu, (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (4, 8) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = mu, (5, 6) = 0, (5, 7) = 0, (5, 8) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = mu+sigma, (6, 7) = 0, (6, 8) = 0, (7, 1) = 0, (7, 2) = 0, (7, 3) = 0, (7, 4) = 0, (7, 5) = 0, (7, 6) = 0, (7, 7) = mu[b], (7, 8) = 0, (8, 1) = 0, (8, 2) = 0, (8, 3) = 0, (8, 4) = 0, (8, 5) = 0, (8, 6) = 0, (8, 7) = 0, (8, 8) = mu[b]})

V := Matrix(8, 8, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = (1-theta[h])*alpha*beta*`ε`*PI/mu, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (4, 8) = alpha*beta*PI/mu, (5, 1) = 0, (5, 2) = 0, (5, 3) = gamma[o], (5, 4) = gamma[1], (5, 5) = 0, (5, 6) = 0, (5, 7) = 0, (5, 8) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 0, (6, 7) = 0, (6, 8) = 0, (7, 1) = 0, (7, 2) = 0, (7, 3) = 0, (7, 4) = 0, (7, 5) = 0, (7, 6) = 0, (7, 7) = 0, (7, 8) = 0, (8, 1) = 0, (8, 2) = 0, (8, 3) = alpha*c*lambda[b]/mu[b], (8, 4) = alpha*c*lambda[b]/mu[b], (8, 5) = 0, (8, 6) = 0, (8, 7) = 0, (8, 8) = 0})

F := Matrix(8, 8, {(1, 1) = mu, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = sigma, (1, 7) = 0, (1, 8) = 0, (2, 1) = 0, (2, 2) = mu, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = gamma[o]+mu, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = gamma[1]+delta+mu, (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (4, 8) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = mu, (5, 6) = 0, (5, 7) = 0, (5, 8) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = mu+sigma, (6, 7) = 0, (6, 8) = 0, (7, 1) = 0, (7, 2) = 0, (7, 3) = 0, (7, 4) = 0, (7, 5) = 0, (7, 6) = 0, (7, 7) = mu[b], (7, 8) = 0, (8, 1) = 0, (8, 2) = 0, (8, 3) = 0, (8, 4) = 0, (8, 5) = 0, (8, 6) = 0, (8, 7) = 0, (8, 8) = mu[b]})

1/V

Error, (in rtable/Power) singular matrix

VectorMatrixMultiply(F, 1/V)

Error, (in rtable/Power) singular matrix

Vinv := MatrixInverse(V, method = pseudo)

Vinv := Matrix(8, 8, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = -1/(-gamma[o]+gamma[1]), (3, 6) = 0, (3, 7) = 0, (3, 8) = gamma[1]*mu[b]/(alpha*c*lambda[b]*(-gamma[o]+gamma[1])), (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 1/(-gamma[o]+gamma[1]), (4, 6) = 0, (4, 7) = 0, (4, 8) = -gamma[o]*mu[b]/(alpha*c*lambda[b]*(-gamma[o]+gamma[1])), (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 0, (5, 6) = 0, (5, 7) = 0, (5, 8) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 0, (6, 7) = 0, (6, 8) = 0, (7, 1) = 0, (7, 2) = 0, (7, 3) = 0, (7, 4) = 0, (7, 5) = 0, (7, 6) = 0, (7, 7) = 0, (7, 8) = 0, (8, 1) = 0, (8, 2) = -`ε`*mu*(-1+theta[h])/(alpha*beta*PI*(`ε`^2*theta[h]^2-2*`ε`^2*theta[h]+`ε`^2+1)), (8, 3) = 0, (8, 4) = mu/(alpha*beta*PI*(`ε`^2*theta[h]^2-2*`ε`^2*theta[h]+`ε`^2+1)), (8, 5) = 0, (8, 6) = 0, (8, 7) = 0, (8, 8) = 0})

`~`[normal]((V.Vinv)^%T-V.Vinv); `~`[normal](V.Vinv.V-V)

Matrix(8, 8, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (4, 8) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 0, (5, 6) = 0, (5, 7) = 0, (5, 8) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 0, (6, 7) = 0, (6, 8) = 0, (7, 1) = 0, (7, 2) = 0, (7, 3) = 0, (7, 4) = 0, (7, 5) = 0, (7, 6) = 0, (7, 7) = 0, (7, 8) = 0, (8, 1) = 0, (8, 2) = 0, (8, 3) = 0, (8, 4) = 0, (8, 5) = 0, (8, 6) = 0, (8, 7) = 0, (8, 8) = 0})

Matrix(8, 8, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (4, 8) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 0, (5, 6) = 0, (5, 7) = 0, (5, 8) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 0, (6, 7) = 0, (6, 8) = 0, (7, 1) = 0, (7, 2) = 0, (7, 3) = 0, (7, 4) = 0, (7, 5) = 0, (7, 6) = 0, (7, 7) = 0, (7, 8) = 0, (8, 1) = 0, (8, 2) = 0, (8, 3) = 0, (8, 4) = 0, (8, 5) = 0, (8, 6) = 0, (8, 7) = 0, (8, 8) = 0})

F.Vinv

Matrix(8, 8, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = -(gamma[o]+mu)/(-gamma[o]+gamma[1]), (3, 6) = 0, (3, 7) = 0, (3, 8) = (gamma[o]+mu)*gamma[1]*mu[b]/(alpha*c*lambda[b]*(-gamma[o]+gamma[1])), (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = (gamma[1]+delta+mu)/(-gamma[o]+gamma[1]), (4, 6) = 0, (4, 7) = 0, (4, 8) = -(gamma[1]+delta+mu)*gamma[o]*mu[b]/(alpha*c*lambda[b]*(-gamma[o]+gamma[1])), (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 0, (5, 6) = 0, (5, 7) = 0, (5, 8) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 0, (6, 7) = 0, (6, 8) = 0, (7, 1) = 0, (7, 2) = 0, (7, 3) = 0, (7, 4) = 0, (7, 5) = 0, (7, 6) = 0, (7, 7) = 0, (7, 8) = 0, (8, 1) = 0, (8, 2) = -mu[b]*`ε`*mu*(-1+theta[h])/(alpha*beta*PI*(`ε`^2*theta[h]^2-2*`ε`^2*theta[h]+`ε`^2+1)), (8, 3) = 0, (8, 4) = mu[b]*mu/(alpha*beta*PI*(`ε`^2*theta[h]^2-2*`ε`^2*theta[h]+`ε`^2+1)), (8, 5) = 0, (8, 6) = 0, (8, 7) = 0, (8, 8) = 0})

 

Download matinv.mw

Here is another way to handle this example,

restart;

kernelopts(version);

   Maple 2019.1, X86 64 LINUX, May 21 2019, Build ID 1399874

fprime := x-> x^6-(3/2)*x^5+2*x^4+(5/2)*x^3-7*x^2+2:
f := unapply(simplify(int(fprime(x), x)), x):
g := unapply(expand(f(x^2+2*x)), x):

sols := [solve({x>-infinity, x<infinity,
                 diff(g(x),x)=0, diff(g(x),x,x)<>0},x)]:

evalf(sols);

      [{x = -1.}, {x = 0.2834077704}, {x = -2.283407770}, 
       {x = 0.4977473282}, {x = -2.497747328}, {x = -0.3051050611}, 
       {x = -1.694894939}]

sort(map[2](eval,x,evalf(sols)));

    [-2.497747328, -2.283407770, -1.694894939, -1., -0.3051050611, 
      0.2834077704, 0.4977473282]

I was able to compile the C file and link the shared library, in a Linux terminal.

Using the commandline interface for Maple 2019.0 in that Linux terminal I was get it to run from within Maple after a call to define_external.

I added a few lines in the C source to printf the entries, before the computational portion. And I threw in a printf("\n") before the return(0).

Here is what my terminal session looked like:

restart:                                                                                     

loc:=cat(kernelopts(homedir),"/gauss.so"):                                                   

mygauss:=define_external('main',RETURN::integer[4],LIB=loc):                                 

mygauss();                                                                                   

Enter the order of matrix: 
Enter the elements of augmented matrix row-wise:

A[1][1] : A[1][2] : 
A[1][1] = 2.000000

A[1][2] = 3.000000

The solution is: 

x1=1.500000	
                          0

As Carl has mentioned, using scanf in C to input the entries is unusual. It's awkward, and possibly not acting just as intended (prematurely accepting and receiving input values, say). It behaved a little weird for me -- I'm not sure that I'd easily be able to enter a 2x2 example with float entries. You'd be far better off specifying Matrix arguments in the define_external, both to pass in the data, and to act in-place to hold the result.

But the first issue is whether the define_external will succeed for you.

Are you 100% sure that the LIB location that you're using actually points to the correct location of the shared library? Is your directory actually called "/home/user/Desktop"? Where is the shared library file? Is it somewhere like, say, "/home/radaar/gauss.so" ?

Is this acceptable?

I wrote it to try and handle only the kind of indexing in your example, eg, C[1][i,m] etc. No doubt something recursive could be concocted to handle arbitrary depth of indexing and arbitrary degree of indexing (per depth).

restart;

ee := sum(sum(C[1][i, m]*C[1][j, k]*v[0][k, m], k = 1 .. N), m = 1 .. M)-A__0*H__1*(sum(sum(C[1][i, m]*C[1][j, k]*v[0][k, m], k = 1 .. N), m = 1 .. M))-A__1*F*(sum(C[2][i, k]*psi[theta][k, j], k = 1 .. N))+A__1*K__4*(sum(C[1][j, k]*psi[theta][k, i], k = 1 .. N))/r__i+A__0*K__4*(sum(C[1][j, k]*v[0][k, i], k = 1 .. N))/(2*r__i)

sum(sum(C[1][i, m]*C[1][j, k]*v[0][k, m], k = 1 .. N), m = 1 .. M)-A__0*H__1*(sum(sum(C[1][i, m]*C[1][j, k]*v[0][k, m], k = 1 .. N), m = 1 .. M))-A__1*F*(sum(C[2][i, k]*psi[theta][k, j], k = 1 .. N))+A__1*K__4*(sum(C[1][j, k]*psi[theta][k, i], k = 1 .. N))/r__i+(1/2)*A__0*K__4*(sum(C[1][j, k]*v[0][k, i], k = 1 .. N))/r__i

new:=subsindets(ee,
                And(indexed,satisfies(u->op(0,u)::indexed)),
                u->nprintf(`%a__%a__%s`,
                           op([0,0],u),op([0,..],u),
                           cat(op(1,u),
                               seq(cat(`,`,p),p=[op(2..,u)]))));

sum(sum(`C__1__i,m`*`C__1__j,k`*`v__0__k,m`, k = 1 .. N), m = 1 .. M)-A__0*H__1*(sum(sum(`C__1__i,m`*`C__1__j,k`*`v__0__k,m`, k = 1 .. N), m = 1 .. M))-A__1*F*(sum(`C__2__i,k`*`psi__theta__k,j`, k = 1 .. N))+A__1*K__4*(sum(`C__1__j,k`*`psi__theta__k,i`, k = 1 .. N))/r__i+(1/2)*A__0*K__4*(sum(`C__1__j,k`*`v__0__k,i`, k = 1 .. N))/r__i

latex(eval(new,1));

\sum _{m=1}^{M} \left( \sum _{k=1}^{N}C_{\mbox {{\tt 1}}_{\mbox {{\tt

i,m}}}}\,C_{\mbox {{\tt 1}}_{\mbox {{\tt j,k}}}}\,v_{\mbox {{\tt 0}}_{
\mbox {{\tt k,m}}}} \right) -A_{0}\,H_{1}\,\sum _{m=1}^{M} \left(
\sum _{k=1}^{N}C_{\mbox {{\tt 1}}_{\mbox {{\tt i,m}}}}\,C_{\mbox {{
\tt 1}}_{\mbox {{\tt j,k}}}}\,v_{\mbox {{\tt 0}}_{\mbox {{\tt k,m}}}}
 \right) -A_{1}\,F\sum _{k=1}^{N}C_{\mbox {{\tt 2}}_{\mbox {{\tt i,k}}
}}\,\psi_{\theta_{\mbox {{\tt k,j}}}}+{\frac {A_{1}\,K_{4}\,\sum _{k=1
}^{N}C_{\mbox {{\tt 1}}_{\mbox {{\tt j,k}}}}\,\psi_{\theta_{\mbox {{
\tt k,i}}}}}{r_{i}}}+1/2\,{\frac {A_{0}\,K_{4}\,\sum _{k=1}^{N}C_{
\mbox {{\tt 1}}_{\mbox {{\tt j,k}}}}\,v_{\mbox {{\tt 0}}_{\mbox {{\tt
k,i}}}}}{r_{i}}}

 

Download end_ac.mw

seq(x[i], i = 1 .. 4);
 
                x[1], x[2], x[3], x[4]

seq(cat(x,i), i = 1 .. 4);

                     x1, x2, x3, x4

Here are the plots of Re(c1) and Im(c1).

I used two different formulas for A and S1 (since above you have used ones which do not correspond to the correct translation from Mathematica, that you asked about yesterday).

It looks as if the revised formulas produce results with zero imaginary component (and so might be what you were expecting?).

restart:

kernelopts(version);

`Maple 13.00, X86 64 LINUX, Apr 13 2009, Build ID 397624`

(1)

Digits:=15:

d1:=0.2:L1:=0.2:L2:=0.2:B1:=0.7:B:=1:beta:=0.01:

d2:=0.6:m:=0.1:k:=0.1:

h:=z->piecewise( z<=d1,    1,

                z<=d1+L1,   1-(gamma1/(2))*(1 + cos(2*(Pi/L1)*(z-d1-L1/2))),

                z<=B1-L2/2,  1 ,          

                z<=B1,  1-(gamma2/(2))*(1 + cos(2*(Pi/L2)*(z - B1))),

                z<=B1+L2/2,  1-(gamma2/(2))*(1 + cos(2*(Pi/L2)*(z - B1))),

               z<=B,    1):

A:=(-m^2/4)-(1/(4*k)):

S1:=(h(z)^2)/(4*A)-ln(A*h(z)^2+1)*(1+h(z)^2)/(4*A):

b1:=1/S1:

c1:=Int(b1,z=0..1,method=_Dexp,digits=15,epsilon=1e-4):

plot([seq(eval(Re(c1),gamma2=j),j in[0,0.02,0.06])],gamma1=0.02..0.1,
     legend = [gamma2 = 0.0,gamma2 = 0.02,gamma2 = 0.04],
     linestyle = [solid,dash,dot],
     color = [black,black,black],
     labels=[gamma1,'Re(c1)'],
     gridlines=false, axes=boxed);

 

plot([seq(eval(Im(c1),gamma2=j),j in[0,0.02,0.06])],gamma1=0.02..0.1,
     legend = [gamma2 = 0.0,gamma2 = 0.02,gamma2 = 0.04],
     linestyle = [solid,dash,dot],
     color = [black,black,black],
     labels=[gamma1,'Im(c1)'],
     gridlines=false, axes=boxed);

 

A:=(-m^2/4)-(1/4*k):

S1:=(h(z)^2)/4*A-ln(A*h(z)^2+1)*(1+h(z)^2)/4*A:

b1:=1/S1:

c1:=Int(b1,z=0..1,method=_Dexp,digits=15,epsilon=1e-4):

plot([seq(eval(Re(c1),gamma2=j),j in[0,0.02,0.06])],gamma1=0.02..0.1,
     legend = [gamma2 = 0.0,gamma2 = 0.02,gamma2 = 0.04],
     linestyle = [solid,dash,dot],
     color = [black,black,black],
     labels=[gamma1,'Re(c1)'],
     gridlines=false, axes=boxed);

 

plot([seq(eval(Im(c1),gamma2=j),j in[0,0.02,0.06])],gamma1=0.02..0.1,
     legend = [gamma2 = 0.0,gamma2 = 0.02,gamma2 = 0.04],
     linestyle = [solid,dash,dot],
     color = [black,black,black],
     labels=[gamma1,'Im(c1)'],
     gridlines=false, axes=boxed);

 

 

Download waseem_plot2.mw

These two approaches work in my Maple 2018.2,

restart;

kernelopts(version);

`Maple 2018.2, X86 64 LINUX, Nov 16 2018, Build ID 1362973`

vals:=[k=1.3806E-23, T=300, h__bar=1.05456E-34,
       L__z=6E-9, m__hhw=3.862216E-31,
       E__hh0=3.012136E-21, E__hh1=1.185628E-20]:

p := k*T*m__hhw*(ln(1+exp(E__fv-E__hh0)/(k*T))
                 +ln(1+exp(E__fv-E__hh1)/(k*T)))/(Pi*h__bar^2*L__z)
     = 0.3e25; #convert(0.3e25, rational,exact);

k*T*m__hhw*(ln(1+exp(E__fv-E__hh0)/(k*T))+ln(1+exp(E__fv-E__hh1)/(k*T)))/(Pi*h__bar^2*L__z) = 0.3e25

[solve({combine(p), E__fv<0},{E__fv})] assuming E__fv::real, op(vals);
eval(%,vals);

[{E__fv = E__hh0+ln(.5000000000*(-1.*exp(E__hh0-1.*E__hh1)*T*k-1.*k*T+((exp(E__hh0-1.*E__hh1))^2*T^2*k^2+4.*exp(E__hh0-1.*E__hh1)*exp(0.9424777962e25*h__bar^2*L__z/(k*T*m__hhw))*T^2*k^2-2.*exp(E__hh0-1.*E__hh1)*T^2*k^2+k^2*T^2)^(1/2))/exp(E__hh0-1.*E__hh1))}]

[{E__fv = -48.46001884}]

subsindets(combine(p),specfunc(ln),u->(ln(expand(op(u))))) assuming E__fv::real, op(vals);
eval(%,vals);
solve({%, E__fv<0}, {E__fv});

k*T*m__hhw*ln(1+exp(E__fv)/(exp(E__hh1)*k*T)+exp(E__fv)/(exp(E__hh0)*k*T)+(exp(E__fv))^2/(exp(E__hh0)*k^2*T^2*exp(E__hh1)))/(Pi*h__bar^2*L__z) = 0.3e25

0.7631009085e25*ln(1+0.4828818388e21*exp(E__fv)+0.5829371757e41*(exp(E__fv))^2) = 0.3e25

{E__fv = -48.46001884}

 

Download solve_M2018.2.mw

First 149 150 151 152 153 154 155 Last Page 151 of 336