tomleslie

5358 Reputation

17 Badges

9 years, 301 days

MaplePrimes Activity


These are replies submitted by tomleslie

Use the big green up-arrow in the Mapleprimes toolbar to post your worksheet. This makes our life a lot easier

From the information you have supplied, there appear to be multiple undefined quantities - eg sigma, rho, g, h, nu, whihc will have to be given numeric values before fsolve() can be expected to work.

It *may* be better to use the 'complex' option with fsolve() for your problem - see the help at ?fsolve.details

@Kitonum 

Shouldn't 'v' be assigned [1, 2, 2] rather than [1, 2, 3] - results in cB=[6, 12, 13]

I have no knowledge of the "Clenshaw derivative" or how it is computed - and I don't plan to read the accompanying documentation to find out. I have "simplified" your code - hopefully without changing the underlying logic.

In your original code there are a few things which don't make much sense. For example, you define Vectors using the initialization

Vector(z..Nm+1+z);

Bad news - the start index for vectors is always 1. To avoid this issue  have redefined these as 1-D arrays, whihc can use any contiguous sange of integers as indexes.

There are a couple of other places (highlighted in the attached) where you rely on the fact that uninitialized quantities will evaluate to zero. This *may* be correct, but the "recursive" nature of the associated calculations means that all values will be determined by these uninitialized quantities. I'd call this seriously high-risk programming

Anyhow, for what it is worth, you might want to check out the attached

  restart:
  Clenshaw_Dx_1D:= proc(z,C,Nm,s)
                        local k,
                              A:= Array(z..Nm+1+z),
                              B:= Array(z..Nm+1+z);
                      #
                      # First iteration of the following loop
                      # both A[k+1] and A[k+2] will evaluate
                      # to zero
                      #
                      # Second iteration of this loop A[k+2] will
                      # evaluate to zero
                      #
                      # This *may* be "desired" behaviour, but just
                      # seems *odd*
                      #
                        for k from Nm-1+z by -1 to 1+z do
                            A[k]:= C[k] + 2*s*A[k+1] - A[k+2]:
                        od:
                      #
                      # First iteration of the following loop
                      # both A[k+1], B[k+1] and B[k+2] will
                      # evaluate to zero
                      #
                      # Second iteration of this loop B[k+2] will
                      # evaluate to zero
                      #
                      # This *may* be "desired" behaviour, but just
                      # seems *odd*
                      #
                        for k from Nm-1+z by -1 to 1+z do
                            B[k]:= 2*A[k+1] + 2*s*B[k+1] - B[k+2]:
                        od:
                        return A[z+1]/2 + s*B[1+z] - B[2+z]:
                   end:

  Chebycoeff1D:= proc(express, Nn, A, B)
                      local C2:= NULL,
                            Cfac:= Array(1..Nn),
                            k, K;
                      for k from 1 to Nn do
                          Cfac[k]:= express( cos(Pi/Nn*(k-0.5))*A+B );
                      od;              
                      for K from 1 to Nn do
                          C2:= C2, (2/Nn)*add
                                          ( Cfac[k]*cos(Pi*(K-1)/Nn*(k-0.5)),
                                            k=1..Nn
                                          );
                      od:
                      return [C2]:
                 end:

  nn:= 1.0:
  Lc:= 0.0:
  Rc:= 1.0:
  xM:= 10: # A guess!
  func:= x -> nn*sin(2*Pi*x)+0.5*nn*sin(Pi*x):
  C:= Array( Chebycoeff1D(func, xM+1, 0.5*(Rc-Lc),0.5*(Rc+Lc))):
  Clen1D_Dx:=Clenshaw_Dx_1D(0, C, xM+1, 0.0);

.9252046157

(1)

 


 

Download Clen.mw

Use the big green up-arrow in the Mapleprimes toolbar to attach your worksheet

BISSEC := proc (P, U, V)\
              #
              # added  packages which this procedure uses
              # so that the latter is "self-contained"
              #
                uses LinearAlgebra, plots;

                local a, b, eq1, eq2, M1, M2, t, PU, PV, bissec1, bissec2;
                a := (P-U)/Norm(P-U, 2)+(P-V)/Norm(P-V, 2);
                M1 := P+a*t;
                b := (P-U)/Norm(P-U, 2)-(P-V)/Norm(P-V, 2);
                M2 := P+b*t;
              #
              # Given M1, M2, defined as above, what is intended
              # by M1[1], M1[2], M2[1], M2[2] in the equations
              # below
              #

                eq1 := op(eliminate({x = M1[1], y = M1[2]}, t));
                eq2 := op(eliminate({x = M2[1], y = M2[2]}, t));
                P := convert(P, list);
                U := convert(U, list);
                V := convert(V, list);
                PU := plot([P, U]);
                PV := plot([P, V]);
                bissec1 := implicitplot(op(eq1[2]), x = 0 .. 5, y = 0 .. 10, color = red);
                bissec2 := implicitplot(op(eq2[2]), x = 0 .. 5, y = 0 .. 10, color = green);
                display([bissec1, bissec2, PU, PV], scaling = constrained)
           end proc:

 

@imparter 

what further data do you need????

@imparter 

I want to plot the 2 dimension plot for different values of M . I have three coupled difference scheme .> i am attaching the codes  

Except I have no idea what you want to plot against what!

In the attached, I demonstrate that

  1. In your worksheet you *seem* to generate 270 equations in 272 unknowns.
  2. The unknowns are  C[1..9, 1..10],  T[1..9,1..10], U[1..9,1..10], M, and Gr.
  3. If I provide values for 'M' and 'Gr', this reduces to 270 equations in 270 unknowns - which Maple solves

But what, if anything, should be done with this data????

restart;
# Parameter values:
 Pr:=0.71:E:=1:A:=0:Sc:=0.02: K:=1:

a := 0: b := 1: N := 9:
h := (b-a)/(N+1): k := (b-a)/(N+1):

 lambda:= 1/h^2:  lambda1:= 1/k^2:
# Initial conditions
for i from 0 to N do
  U[i, 0] := h*i+1:
end do:


for i from 0 to N do
  T[i, 0] := h*i+1:
end do:
for i from 0 to N do
  C[i, 0] := h*i+1:
end do:

# Boundary conditions
for j from 0 to N+1 do
  U[0, j] := exp(A*j*lambda);
  U[N+1, j] := 0;
  T[0, j] := j*lambda1;
  T[N+1, j] := 0;
  C[0, j] := j*lambda1;
  C[N+1, j] := 0
end do:


#Discretization Scheme
for i to N do
  for j from 0 to N do
    eq1[i, j]:= lambda1*(U[i, j+1]-U[i, j]) = (Gr/2)*(T[i, j+1]+T[i,j])+(Gr/2)*(C[i, j+1]+C[i,j])+(lambda^2/2)*(U[i-1,j+1]-2*U[i,j+1]+U[i+1,j+1]+U[i-1,j]-2*U[i,j]+U[i+1,j])-(M/2)*(U[i,j+1]+U[i,j]) ;
    eq2[i, j]:= lambda1*(T[i, j+1]-T[i, j]) = (1/Pr)*(lambda^2/2)*(T[i,j+1]-2*T[i,j+1]+T[i+1,j+1]+T[i-1,j]-2*T[i,j]+T[i+1,j])+(E*lambda^2)*((U[i+1,j]-U[i,j])^2);
    eq3[i, j]:= lambda1*(C[i, j+1]-C[i, j]) = (1/Sc)*(lambda^2/2)*(C[i,j+1]-2*C[i,j+1]+C[i+1,j+1]+C[i-1,j]-2*C[i,j]+C[i+1,j])+(K/2)*((C[i,j+1]+C[i,j]))  
  end do
end do:
 

#
# Deeermine the unknowns in the system
#
  `union`(  seq(seq( indets( eq1[i,j], name), i=1..N), j=0..N),
            seq(seq( indets( eq2[i,j], name), i=1..N), j=0..N),
            seq(seq( indets( eq3[i,j], name), i=1..N), j=0..N)
          );
#
# And how many unknowns
#
   numelems(%);
#
# And the number of equations
#
  numelems(eq1)+numelems(eq2)+numelems(eq3);

{Gr, M, C[1, 1], C[1, 2], C[1, 3], C[1, 4], C[1, 5], C[1, 6], C[1, 7], C[1, 8], C[1, 9], C[1, 10], C[2, 1], C[2, 2], C[2, 3], C[2, 4], C[2, 5], C[2, 6], C[2, 7], C[2, 8], C[2, 9], C[2, 10], C[3, 1], C[3, 2], C[3, 3], C[3, 4], C[3, 5], C[3, 6], C[3, 7], C[3, 8], C[3, 9], C[3, 10], C[4, 1], C[4, 2], C[4, 3], C[4, 4], C[4, 5], C[4, 6], C[4, 7], C[4, 8], C[4, 9], C[4, 10], C[5, 1], C[5, 2], C[5, 3], C[5, 4], C[5, 5], C[5, 6], C[5, 7], C[5, 8], C[5, 9], C[5, 10], C[6, 1], C[6, 2], C[6, 3], C[6, 4], C[6, 5], C[6, 6], C[6, 7], C[6, 8], C[6, 9], C[6, 10], C[7, 1], C[7, 2], C[7, 3], C[7, 4], C[7, 5], C[7, 6], C[7, 7], C[7, 8], C[7, 9], C[7, 10], C[8, 1], C[8, 2], C[8, 3], C[8, 4], C[8, 5], C[8, 6], C[8, 7], C[8, 8], C[8, 9], C[8, 10], C[9, 1], C[9, 2], C[9, 3], C[9, 4], C[9, 5], C[9, 6], C[9, 7], C[9, 8], C[9, 9], C[9, 10], T[1, 1], T[1, 2], T[1, 3], T[1, 4], T[1, 5], T[1, 6], T[1, 7], T[1, 8], T[1, 9], T[1, 10], T[2, 1], T[2, 2], T[2, 3], T[2, 4], T[2, 5], T[2, 6], T[2, 7], T[2, 8], T[2, 9], T[2, 10], T[3, 1], T[3, 2], T[3, 3], T[3, 4], T[3, 5], T[3, 6], T[3, 7], T[3, 8], T[3, 9], T[3, 10], T[4, 1], T[4, 2], T[4, 3], T[4, 4], T[4, 5], T[4, 6], T[4, 7], T[4, 8], T[4, 9], T[4, 10], T[5, 1], T[5, 2], T[5, 3], T[5, 4], T[5, 5], T[5, 6], T[5, 7], T[5, 8], T[5, 9], T[5, 10], T[6, 1], T[6, 2], T[6, 3], T[6, 4], T[6, 5], T[6, 6], T[6, 7], T[6, 8], T[6, 9], T[6, 10], T[7, 1], T[7, 2], T[7, 3], T[7, 4], T[7, 5], T[7, 6], T[7, 7], T[7, 8], T[7, 9], T[7, 10], T[8, 1], T[8, 2], T[8, 3], T[8, 4], T[8, 5], T[8, 6], T[8, 7], T[8, 8], T[8, 9], T[8, 10], T[9, 1], T[9, 2], T[9, 3], T[9, 4], T[9, 5], T[9, 6], T[9, 7], T[9, 8], T[9, 9], T[9, 10], U[1, 1], U[1, 2], U[1, 3], U[1, 4], U[1, 5], U[1, 6], U[1, 7], U[1, 8], U[1, 9], U[1, 10], U[2, 1], U[2, 2], U[2, 3], U[2, 4], U[2, 5], U[2, 6], U[2, 7], U[2, 8], U[2, 9], U[2, 10], U[3, 1], U[3, 2], U[3, 3], U[3, 4], U[3, 5], U[3, 6], U[3, 7], U[3, 8], U[3, 9], U[3, 10], U[4, 1], U[4, 2], U[4, 3], U[4, 4], U[4, 5], U[4, 6], U[4, 7], U[4, 8], U[4, 9], U[4, 10], U[5, 1], U[5, 2], U[5, 3], U[5, 4], U[5, 5], U[5, 6], U[5, 7], U[5, 8], U[5, 9], U[5, 10], U[6, 1], U[6, 2], U[6, 3], U[6, 4], U[6, 5], U[6, 6], U[6, 7], U[6, 8], U[6, 9], U[6, 10], U[7, 1], U[7, 2], U[7, 3], U[7, 4], U[7, 5], U[7, 6], U[7, 7], U[7, 8], U[7, 9], U[7, 10], U[8, 1], U[8, 2], U[8, 3], U[8, 4], U[8, 5], U[8, 6], U[8, 7], U[8, 8], U[8, 9], U[8, 10], U[9, 1], U[9, 2], U[9, 3], U[9, 4], U[9, 5], U[9, 6], U[9, 7], U[9, 8], U[9, 9], U[9, 10]}

 

272

 

270

(1)

#
# So if one supplies values for 'Gr' and 'M', and
# (assuming the equations are consistent), one ought
# to be able to get values for C[1..9, 1..10],
# T[1..9,1..10], and U[1..9,1..10]
#
# As an example below, choos Gr=1.0 and M=2, then the
# following obtains a 'solution` afer a minute or so
#
  fsolve( eval( [ seq(seq(eq1[i,j], i=1..N),j=0..N),
                  seq(seq(eq2[i,j], i=1..N),j=0..N),
                  seq(seq(eq3[i,j], i=1..N),j=0..N)
                ],
                [Gr=1.0, M=2]
              )
        );

{C[1, 1] = -2.987037136, C[1, 2] = 104.9301573, C[1, 3] = 93.04802705, C[1, 4] = 208.8954429, C[1, 5] = 189.0834863, C[1, 6] = 312.8607833, C[1, 7] = 285.1189458, C[1, 8] = 416.8261238, C[1, 9] = 381.1544054, C[1, 10] = 521.7839323, C[2, 1] = -1.988668177, C[2, 2] = .9877141813, C[2, 3] = 101.9154786, C[2, 4] = -6.878246563, C[2, 5] = 213.7438979, C[2, 6] = -22.66763497, C[2, 7] = 333.4958547, C[2, 8] = -46.38056031, C[2, 9] = 461.1713486, C[2, 10] = -77.02415965, C[3, 1] = -1.989942067, C[3, 2] = 1.988549585, C[3, 3] = -2.987634581, C[3, 4] = 106.8486107, C[3, 5] = -115.6688137, C[3, 6] = 331.3105414, C[3, 7] = -355.8688962, C[3, 8] = 691.2085747, C[3, 9] = -739.4222793, C[3, 10] = 1203.370365, C[4, 1] = -1.991256663, C[4, 2] = 1.990181692, C[4, 3] = -1.990419759, C[4, 4] = .9920094061, C[4, 5] = 103.8256555, C[4, 6] = -217.4182788, C[4, 7] = 546.5632778, C[4, 8] = -900.1701899, C[4, 9] = 1588.973464, C[4, 10] = -2324.842817, C[5, 1] = -1.992611984, C[5, 2] = 1.991814992, C[5, 3] = -1.992013193, C[5, 4] = 1.992012631, C[5, 5] = -2.990025164, C[5, 6] = 108.7636006, C[5, 7] = -328.0424876, C[5, 8] = 876.2897633, C[5, 9] = -1777.870973, C[5, 10] = 3368.833309, C[6, 1] = -1.994008043, C[6, 2] = 1.993449504, C[6, 3] = -1.993607935, C[6, 4] = 1.993607549, C[6, 5] = -1.993607676, C[6, 6] = .9959923528, C[6, 7] = 105.7327162, C[6, 8] = -431.6105694, C[6, 9] = 1305.429451, C[6, 10] = -3079.256053, C[7, 1] = -1.995444858, C[7, 2] = 1.995085244, C[7, 3] = -1.995203986, C[7, 4] = 1.995203744, C[7, 5] = -1.995203839, C[7, 6] = 1.995203839, C[7, 7] = -2.992422269, C[7, 8] = 110.6754873, C[7, 9] = -544.0638699, C[7, 10] = 1851.787366, C[8, 1] = -1.996922445, C[8, 2] = 1.996722229, C[8, 3] = -1.996801346, C[8, 4] = 1.996801217, C[8, 5] = -1.996801280, C[8, 6] = 1.996801280, C[8, 7] = -1.996801280, C[8, 8] = .9999795850, C[8, 9] = 107.6366656, C[8, 10] = -648.4507657, C[9, 1] = -1.998440821, C[9, 2] = 1.998360475, C[9, 3] = -1.998400017, C[9, 4] = 1.998399968, C[9, 5] = -1.998400000, C[9, 6] = 1.998400000, C[9, 7] = -1.998400000, C[9, 8] = 1.998400000, C[9, 9] = -2.994825118, C[9, 10] = 113.5799071, T[1, 1] = 2.051943501, T[1, 2] = 97.95738627, T[1, 3] = 103.8443444, T[1, 4] = 194.9704255, T[1, 5] = 207.3004680, T[1, 6] = 290.4890615, T[1, 7] = 312.4331577, T[1, 8] = 384.3165832, T[1, 9] = 419.3874991, T[1, 10] = 477.1327350, T[2, 1] = 3.051261099, T[2, 2] = .3156375825, T[2, 3] = 99.52103631, T[2, 4] = 4.421410864, T[2, 5] = 192.9771636, T[2, 6] = 13.29159212, T[2, 7] = 280.4121624, T[2, 8] = 29.79061451, T[2, 9] = 358.7128970, T[2, 10] = 58.01439391, T[3, 1] = 3.063349007, T[3, 2] = 1.200903643, T[3, 3] = 2.399487524, T[3, 4] = 95.84724691, T[3, 5] = -86.33294602, T[3, 6] = 275.7136622, T[3, 7] = -255.4296074, T[3, 8] = 530.0402551, T[3, 9] = -491.4054127, T[3, 10] = 843.1739982, T[4, 1] = 3.074188563, T[4, 2] = 1.101289867, T[4, 3] = 3.400748677, T[4, 4] = -.9818812947, T[4, 5] = 99.32446116, T[4, 6] = -184.1274183, T[4, 7] = 459.2626794, T[4, 8] = -709.3437173, T[4, 9] = 1233.664839, T[4, 10] = -1713.063713, T[5, 1] = 3.083762040, T[5, 2] = .9831881216, T[5, 3] = 3.451517590, T[5, 4] = -.1477966969, T[5, 5] = 3.082964308, T[5, 6] = 93.70915571, T[5, 7] = -269.2848378, T[5, 8] = 717.2550908, T[5, 9] = -1404.727662, T[5, 10] = 2608.437721, T[6, 1] = 3.092051461, T[6, 2] = .8389149178, T[6, 3] = 3.506547473, T[6, 4] = -.2950015032, T[6, 5] = 4.109573569, T[6, 6] = -2.343669730, T[6, 7] = 99.44920275, T[6, 8] = -366.0630373, T[6, 9] = 1075.035860, T[6, 10] = -2452.732028, T[7, 1] = 3.099038592, T[7, 2] = .6581520192, T[7, 3] = 3.556674117, T[7, 4] = -.5029570870, T[7, 5] = 4.215897147, T[7, 6] = -1.659336274, T[7, 7] = 4.150112060, T[7, 8] = 91.23612469, T[7, 9] = -445.1216297, T[7, 10] = 1496.921080, T[8, 1] = 3.104704940, T[8, 2] = .4269665478, T[8, 3] = 3.583859979, T[8, 4] = -.8146662578, T[8, 5] = 4.308447820, T[8, 6] = -2.026152581, T[8, 7] = 5.214593930, T[8, 8] = -4.179125799, T[8, 9] = 99.92394335, T[8, 10] = -540.9496436, T[9, 1] = 3.109031750, T[9, 2] = .1265218117, T[9, 3] = 3.555703644, T[9, 4] = -1.304792227, T[9, 5] = 4.323309568, T[9, 6] = -2.652987312, T[9, 7] = 5.296816450, T[9, 8] = -4.016085222, T[9, 9] = 5.529854633, T[9, 10] = 88.66581839, U[1, 1] = .8066155842, U[1, 2] = .9462319949, U[1, 3] = .9416200435, U[1, 4] = .9198276005, U[1, 5] = 1.003947662, U[1, 6] = .9387422967, U[1, 7] = 1.038188225, U[1, 8] = .9745224131, U[1, 9] = 1.062854789, U[1, 10] = 1.015991239, U[2, 1] = .6076183125, U[2, 2] = .8810244983, U[2, 3] = .8549869184, U[2, 4] = .8077689863, U[2, 5] = .9618237121, U[2, 6] = .8226672556, U[2, 7] = 1.013487897, U[2, 8] = .8711931233, U[2, 9] = 1.045566943, U[2, 10] = .9316577535, U[3, 1] = .3967886714, U[3, 2] = .8331786289, U[3, 3] = .7305446300, U[3, 4] = .7130097599, U[3, 5] = .8654089652, U[3, 6] = .7187222790, U[3, 7] = .9203879519, U[3, 8] = .7740626852, U[3, 9] = .9456216441, U[3, 10] = .8480115348, U[4, 1] = .1678668208, U[4, 2] = .8119724757, U[4, 3] = .5774625596, U[4, 4] = .6266175580, U[4, 5] = .7639214519, U[4, 6] = .5767309170, U[4, 7] = .8701229573, U[4, 8] = .5705146484, U[4, 9] = .9568269459, U[4, 10] = .5700495724, U[5, 1] = -0.8577241316e-1, U[5, 2] = .8281443464, U[5, 3] = .3821399738, U[5, 4] = .5835477418, U[5, 5] = .6018197140, U[5, 6] = .5117176572, U[5, 7] = .6886090460, U[5, 8] = .5528802846, U[5, 9] = .6688378493, U[5, 10] = .7053786590, U[6, 1] = -.3712533649, U[6, 2] = .8941781292, U[6, 3] = .1278339947, U[6, 4] = .6033522398, U[6, 5] = .3772524594, U[6, 6] = .4873340111, U[6, 7] = .5097288711, U[6, 8] = .4305241541, U[6, 9] = .6470393924, U[6, 10] = .2960555259, U[7, 1] = -.6963434390, U[7, 2] = 1.024841208, U[7, 3] = -.2066582915, U[7, 4] = .7126784925, U[7, 5] = 0.5845611829e-1, U[7, 6] = .5594771709, U[7, 7] = .2145857759, U[7, 8] = .4822840690, U[7, 9] = .2973915372, U[7, 10] = .5229030098, U[8, 1] = -1.069610010, U[8, 2] = 1.237794494, U[8, 3] = -.6483286081, U[8, 4] = .9473653435, U[8, 5] = -.3990023358, U[8, 6] = .7801707953, U[8, 7] = -.2360001550, U[8, 8] = .6746732293, U[8, 9] = -.1177000901, U[8, 10] = .5839162594, U[9, 1] = -1.500593481, U[9, 2] = 1.554293259, U[9, 3] = -1.231550049, U[9, 4] = 1.355300085, U[9, 5] = -1.056875740, U[9, 6] = 1.224710837, U[9, 7] = -.9309659318, U[9, 8] = 1.129739573, U[9, 9] = -.8316431190, U[9, 10] = 1.056090461}

(2)

 


 

Download bigSys.mw

@Stretto 

You can of curse use sommething from the CurveFitting() package to "fit" the data to a model function, or if youu don't know the model function, just use a 'Spline' fit

@Maple_lover1 

with constraints on variables, and checks for "additional" solutions


 

  restart;

  A:=Matrix(2, 2, [[-0.0001633261895*z[1, 2]^2 + 0.0002805135275*z[1, 2]*z[2, 2] - 0.0001200583046*z[2, 2]^2 + 0.0006934805795*z[1, 1]^2 - 0.001190280265*z[1, 1]*z[2, 1] + 0.00007689977894*z[1, 1]*z[1, 2] - 0.00009937418547*z[1, 1]*z[2, 2] + 0.0005090615773*z[2, 1]^2 - 0.00003303758400*z[2, 1]*z[1, 2] + 0.00005683264925*z[2, 1]*z[2, 2] + 0.7021232886*z[1, 1] - 0.3171553245*z[1, 2] - 0.08291569324*z[2, 1] + 0.04647270631*z[2, 2] - 0.1436869545, 0.0002939068385*z[2, 1]^2 + 0.4237544799*z[1, 1] - 0.03129537402*z[1, 2] - 0.06276282411*z[2, 1] + 0.02529757039*z[2, 2] + 0.0004003811990*z[1, 1]^2 + 0.0002177682527*z[1, 1]*z[1, 2] - 0.0006872086309*z[1, 1]*z[2, 1] - 0.0001976167183*z[1, 1]*z[2, 2] - 0.0001764013184*z[2, 1]*z[1, 2] + 0.0001600777394*z[2, 1]*z[2, 2] - 0.1237363898], [0.00006763201108*z[2, 1]*z[1, 2] - 0.0001020812322*z[1, 2]*z[2, 2] - 0.00001554113990*z[2, 1]*z[2, 2] - 0.00003577693711*z[1, 1]*z[1, 2] + 0.0004330743651*z[1, 1]*z[2, 1] - 0.00001941220415*z[1, 1]*z[2, 2] - 0.01736180925 + 0.5623450996*z[2, 1] - 0.2353707048*z[2, 2] - 0.0003226356619*z[1, 1]^2 + 0.00007598605473*z[1, 2]^2 - 0.0001392051452*z[2, 1]^2 + 0.00003283047567*z[2, 2]^2 + 0.04653058230*z[1, 1] - 0.03026711709*z[1, 2], -0.00008037012799*z[2, 1]^2 + 0.03994641178*z[1, 1] - 0.02291248064*z[1, 2] + 0.3140461555*z[2, 1] + 0.01853659924*z[2, 2] - 0.0001862737861*z[1, 1]^2 - 0.0001013147396*z[1, 1]*z[1, 2] + 0.0002500356011*z[1, 1]*z[2, 1] + 0.00005403916772*z[1, 1]*z[2, 2] + 0.00008206914192*z[2, 1]*z[1, 2] - 0.00004377396755*z[2, 1]*z[2, 2] - 0.01370765196]]):

#
# Extract equations from above matrix and define
# constraint ranges for variables
#
  eqs:=[ entries(A, 'nolist') ]:
  rngs:={ seq
          ( j=0..1,
            j in `union`(indets~([entries(A,nolist)])[])
          )
        };

#
# Obtain a solution  of eqs with constraints
#
  sol1:= fsolve( eqs,
                 rngs
               );
#
# Check to see if there are any other solutions
# by telling fsolve() to 'avoid' the one previously
# obtained
#
  sol2:= fsolve(  eqs,
                  rngs,
                  avoid = {sol1}
               );
#
# Double-check by running DirectSearch:-SolveEquations()
# with the 'AllSolutions' option. Just o see if there are
# any other solution, anywhere
#
# OP (probably?) won't be able to do this
#
  DirectSearch:-SolveEquations( eqs,
                                rngs,
                                AllSolutions=true
                              );

rngs := {z[1, 1] = 0 .. 1, z[1, 2] = 0 .. 1, z[2, 1] = 0 .. 1, z[2, 2] = 0 .. 1}

 

sol1 := {z[1, 1] = .3117132485, z[1, 2] = .2328518749, z[2, 1] = 0.2064174947e-1, z[2, 2] = 0.7118281938e-2}

 

sol2 := fsolve([-0.1633261895e-3*z[1, 2]^2+0.2805135275e-3*z[1, 2]*z[2, 2]-0.1200583046e-3*z[2, 2]^2+0.6934805795e-3*z[1, 1]^2-0.1190280265e-2*z[1, 1]*z[2, 1]+0.7689977894e-4*z[1, 1]*z[1, 2]-0.9937418547e-4*z[1, 1]*z[2, 2]+0.5090615773e-3*z[2, 1]^2-0.3303758400e-4*z[2, 1]*z[1, 2]+0.5683264925e-4*z[2, 1]*z[2, 2]+.7021232886*z[1, 1]-.3171553245*z[1, 2]-0.8291569324e-1*z[2, 1]+0.4647270631e-1*z[2, 2]-.1436869545, 0.6763201108e-4*z[2, 1]*z[1, 2]-0.1020812322e-3*z[1, 2]*z[2, 2]-0.1554113990e-4*z[2, 1]*z[2, 2]-0.3577693711e-4*z[1, 1]*z[1, 2]+0.4330743651e-3*z[1, 1]*z[2, 1]-0.1941220415e-4*z[1, 1]*z[2, 2]-0.1736180925e-1+.5623450996*z[2, 1]-.2353707048*z[2, 2]-0.3226356619e-3*z[1, 1]^2+0.7598605473e-4*z[1, 2]^2-0.1392051452e-3*z[2, 1]^2+0.3283047567e-4*z[2, 2]^2+0.4653058230e-1*z[1, 1]-0.3026711709e-1*z[1, 2], 0.2939068385e-3*z[2, 1]^2+.4237544799*z[1, 1]-0.3129537402e-1*z[1, 2]-0.6276282411e-1*z[2, 1]+0.2529757039e-1*z[2, 2]+0.4003811990e-3*z[1, 1]^2+0.2177682527e-3*z[1, 1]*z[1, 2]-0.6872086309e-3*z[1, 1]*z[2, 1]-0.1976167183e-3*z[1, 1]*z[2, 2]-0.1764013184e-3*z[2, 1]*z[1, 2]+0.1600777394e-3*z[2, 1]*z[2, 2]-.1237363898, -0.8037012799e-4*z[2, 1]^2+0.3994641178e-1*z[1, 1]-0.2291248064e-1*z[1, 2]+.3140461555*z[2, 1]+0.1853659924e-1*z[2, 2]-0.1862737861e-3*z[1, 1]^2-0.1013147396e-3*z[1, 1]*z[1, 2]+0.2500356011e-3*z[1, 1]*z[2, 1]+0.5403916772e-4*z[1, 1]*z[2, 2]+0.8206914192e-4*z[2, 1]*z[1, 2]-0.4377396755e-4*z[2, 1]*z[2, 2]-0.1370765196e-1], {z[1, 1], z[1, 2], z[2, 1], z[2, 2]}, {z[1, 1] = 0 .. 1, z[1, 2] = 0 .. 1, z[2, 1] = 0 .. 1, z[2, 2] = 0 .. 1}, avoid = {{z[1, 1] = .3117132485, z[1, 2] = .2328518749, z[2, 1] = 0.2064174947e-1, z[2, 2] = 0.7118281938e-2}})

 

Matrix(%id = 18446744074887932494)

(1)

#
# Insert solution values into matrix
# in the required order
#
  Z:=Matrix(2,2,(i,j)->eval(z[i,j], sol1 ));

Matrix(2, 2, {(1, 1) = .3117132485, (1, 2) = .2328518749, (2, 1) = 0.2064174947e-1, (2, 2) = 0.7118281938e-2})

(2)

 


 

Download matSol2.mw

@syhue 

  1. Use of the &^ operator is restricted to the case of  modular exponentiation, which you don't appear to be doing? Thus something like rand( (2&^50 mod 1234)..(2&^51 mod 1234))() makes perfect sense, but rand( 2&^50..2&^51)() makes no sense at all (and doesn't work)
  2. To choose a random entry from your factors, one way is to
    1. Construct your factors as a list
    2. Define a random number generator which will produce integers between 1 and the number of elements in the list - essentially a "random index" into the list
    3. Select the element with this "random index"

See the attached

NULL

NULL

with(RandomTools)

r := rand(2^50 .. 2^51); p := nextprime(r()); q := nextprime(r())

1822418686289087

 

1359743522261933

(1)

n := p*q

2478022003530687861300173425171

(2)

a := (p-1)/(q-1)

911209343144543/679871761130966

(3)

r := numer(a); s := denom(a)

911209343144543

 

679871761130966

(4)

w := s*(p-1)

1239011001765342339568982437076

(5)

f := [seq(j[1], `in`(j, ifactors(w)[2]))]; u := f[(rand(1 .. numelems(f)))()]

[2, 2777, 8429, 122819, 38928371, 2767779257]

 

38928371

(6)

 

Download resp.mw

 

 

 

@acer 

I usually just check the output residuals - if these are "small enough" I reckon I'm good to go. I have noticed in the past that it will output the best(?) solution it can find, and report residuals rather large residuals

@acer 

I missed the lexicographic order possibility -  my bad

I actually submitted an SCR about this problem some years ago (can't remember when).

The issue does not occur in Maple 18, but was "introduced" in Maple 2015 and has been in every version since, so people have been tripping over this one for ~4 years

Although I did the original worksheet in Maple 2016 (becuase that's what the OP has), the behaviour in Maple 2019 is identical :-(

of namespace', function overloading and variable scope are basic features of any programming language. If you really do not understand these concepts then I respectfully suggest that you keep reading the attached toy example until enlightenment occurs

  restart;
#
# Define som simple "top-level" functions
# and check how they evaluate
#
  lambda:= x->x^2:
  phi:= x->x^3:
  pi:= x->x^4:
  lambda(2), phi(2), pi(2);
#
# Now load the NumberTheory package which contains
# its own definitions for the functions lambda, phi,
# and pi. These will 'overload' the previous
# definitions, so check how these now evaluate
#
  with(NumberTheory):
  lambda(2), phi(2), pi(2);
#
# If you wish to retrieve/use the function definitions
# performed at the top-level, rather than the overloaded
# versions introduced by the NumberTheory package, then
# prefix the function names with ':-', as in
#
  :-lambda(2), :-phi(2), :-pi(2);

4, 8, 16

 

1, 1, 1

 

4, 8, 16

(1)

 

Download nameSpace.mw

@Teep 

Attached is a minor revision of my original post which runs with no problems in Mpale 2018

  restart:
  kernelopts(version);
  xmin:=1: xmax:=3: ymin:=2.5: ymax:=4.5:
  f:= (nx, ny)-> plots:-pointplot
                        ( [ seq
                            ( seq
                              ( [i, j],
                                i=xmin..xmax, (xmax-xmin)/nx
                              ),
                              j=ymin..ymax, (ymax-ymin)/ny
                            )
                          ],
                          symbol=solidcircle,
                          symbolsize=20,
                          color=red,
                          scaling=constrained
                        ):

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

(1)

#
# A few test cases for different grid spacings
#
  f(4,4);
  f(4,6);
  f(10,3);

 

 

 

 

Download doGrids18.mw

1 2 3 4 5 6 7 Last Page 1 of 149