tomleslie

5348 Reputation

15 Badges

9 years, 293 days

MaplePrimes Activity


These are replies submitted by tomleslie

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

@brian bovril .

provided commands are issued in the following sequence

fopen(yourFileName)
writeTo(yourFileName)
fclose(yourFileName)

The file will be opened (and "locked") when the fopen() command is executeded. It will remain open (and "locked") until the fclose() command is executed.

Only thing I can think of is that you are trying to access the file  before the executing fclose().
 

You could do it either way!

If you are going to run this on multiple data sets, then creating a procedure within the package to perform the plots will probably save you some typing

Within the package, I'd probably separate the "fitting" and the "plotting" into two separate procedures - because sometimes one might want the former and not the latter.

See the attached, which is a slight extension of the worksheet I posted previously

  restart;
  Mat:= module()
               description "My Package";
               option package;
               export fitData, doPlots;
               fitData:=proc( xd::list, yd::list, fittype::name,
                              x::algebraic, c::algebraic
                            )
                            uses Statistics, CurveFitting:
                            if   fittype = modfit
                            then return Fit(c, xd, yd, x);
                            elif fittype = lsqfit
                            then return LeastSquares( xd, yd, x, curve=c)
                            fi:
               end proc;
               doPlots:=proc( xd::list, yd::list, cur::algebraic )
                            uses plots:
                            display
                            (  [ plot
                                 ( cur,
                                   x=min(xl)..max(xl),
                                   color=blue,
                                   title=typeset("Fit = ", cur)
                                 ),
                                 plot
                                 ( xl, yl,
                                   style=point,
                                   symbol=solidcircle,
                                   symbolsize=20,
                                   color=red
                                 )
                               ],
                               size=[600,600]
                            );
                        end proc;
        end module:
  with(Mat):
#
# test data
#
  xl:=[1,2,3,4]: yl:=[2,4,5,9]:
#
# Use Statistics:-Fit() with both linear
# and quadratic models
#
  f1:= fitData(xl, yl, 'modfit', x, 'a*x+b');
  doPlots(xl, yl, f1);
  f2:= fitData(xl, yl, 'modfit', x, 'a*x^2+b*x+c');
  doPlots(xl, yl, f2);
#
# Use CurveFitting:-LeastSquares() with
# both linear and quadratic models, then
# plot the fit
#  
  f3:= fitData(xl, yl, 'lsqfit', x, 'a*x+b');
  doPlots(xl, yl, f3);
  f4:= fitData(xl, yl, 'lsqfit', x, 'a*x^2+b*x+c');
  doPlots(xl, yl, f4);

HFloat(2.200000000000001)*x-HFloat(0.5000000000000031)

 

 

HFloat(0.49999999999999917)*x^2-HFloat(0.29999999999999444)*x+HFloat(1.9999999999999911)

 

 

-1/2+(11/5)*x

 

 

2-(3/10)*x+(1/2)*x^2

 

 

 

Download fitStuff2.mw

@pooyan1990 

It is possible to produce a 'symbolic' solution for the ODE - see the attached - but maybe this is too complicated to be useful?

This particular problem seems to fall into the category of those which should be solved numerically!


 

restart;

#B:=10;
#H:=10;
#L__1:=4;

eq0:=R__i=(1/2)*B+L__1;
eq01:=R__f=C+(1/2)*B+L__1;
eq02:=theta__1=arctan((H/2)/(B/2 + L__1));
eq03:=theta__2=Pi - arctan((H/2)/(B/2 - L__1));
eq43:= v__beta(r,theta,beta)= v__m(beta)*(1-(r^2/((R__max(theta,beta))^2)));
eq39:=R(beta)=R__i+2*beta*(R__f-R__i)/Pi;
eq39a:=simplify(subs([eq0,eq01], eq39));
eq38:=R__max(theta,beta)=R__max(theta,0)*(R(beta)/R__i);
#
# Replace eq41 with the simple requirement that
# R__max(theta,0)=phi(theta - Rouben's suggestion
#
eq41_orig:=R__max(theta,0)=piecewise(theta<=rhs(eq02),(B/2+L__1)/cos(theta),theta<=rhs(eq03),(H/2)/sin(theta),rhs(eq03)<theta,(B/2-L__1)/cos(Pi-theta));

eq41:=R__max(theta,0)=phi(theta);

eq38a:= simplify(subs([eq41,eq39a,eq0], eq38));
eq44 := v__m(beta) = 1/(2*int((1 - r^2/R__max(theta, beta)^2)*r, [r = 0 .. R__max(theta, beta), theta = 0 .. Pi]));
eq44a:=simplify(subs(eq38a,eq44)) assuming B/L__1 >= 2 , H > 0 , B > 0 , L__1 > 0;
eq44b:= v__m(beta)= rhs(eq44a) * B * H;
eq47:=subs([eq44b,eq38a],eq43);

eq0 := R__i = (1/2)*B+L__1

 

eq01 := R__f = C+(1/2)*B+L__1

 

eq02 := `&theta;__1` = arctan(H/(2*((1/2)*B+L__1)))

 

eq03 := `&theta;__2` = Pi-arctan(H/(2*((1/2)*B-L__1)))

 

eq43 := `v__&beta;`(r, theta, beta) = v__m(beta)*(1-r^2/R__max(theta, beta)^2)

 

eq39 := R(beta) = R__i+2*beta*(R__f-R__i)/Pi

 

eq39a := R(beta) = (1/2)*B+L__1+2*beta*C/Pi

 

eq38 := R__max(theta, beta) = R__max(theta, 0)*R(beta)/R__i

 

R__max(theta, 0) = piecewise(theta <= arctan((1/2)*H/((1/2)*B+L__1)), ((1/2)*B+L__1)/cos(theta), theta <= Pi-arctan((1/2)*H/((1/2)*B-L__1)), (1/2)*H/sin(theta), Pi-arctan((1/2)*H/((1/2)*B-L__1)) < theta, -((1/2)*B-L__1)/cos(theta))

 

R__max(theta, 0) = phi(theta)

 

R__max(theta, beta) = ((B+2*L__1)*Pi+4*beta*C)*phi(theta)/((B+2*L__1)*Pi)

 

v__m(beta) = 2/(int(R__max(theta, beta)^2, theta = 0 .. Pi))

 

v__m(beta) = 2*(B+2*L__1)^2*Pi^2/((B*Pi+4*C*beta+2*L__1*Pi)^2*(int(phi(theta)^2, theta = 0 .. Pi)))

 

v__m(beta) = 2*(B+2*L__1)^2*Pi^2*B*H/((B*Pi+4*C*beta+2*L__1*Pi)^2*(int(phi(theta)^2, theta = 0 .. Pi)))

 

v__beta(r, theta, beta) = 2*(B+2*L__1)^2*Pi^2*B*H*(1-r^2*(B+2*L__1)^2*Pi^2/(((B+2*L__1)*Pi+4*beta*C)^2*phi(theta)^2))/((B*Pi+4*C*beta+2*L__1*Pi)^2*(int(phi(theta)^2, theta = 0 .. Pi)))

(1)

eq37:=diff(v__r(r, theta, beta),r)*r*(R__f-r*cos(theta))+v__r(r,theta,beta)*(R__f-2*r*cos(theta))+r*diff(v__beta(r,theta,beta),beta)=0;
eq37a:=simplify(expand(subs(eq47, eq37)));

(diff(v__r(r, theta, beta), r))*r*(R__f-r*cos(theta))+v__r(r, theta, beta)*(R__f-2*r*cos(theta))+r*(diff(v__beta(r, theta, beta), beta)) = 0

 

(-2*((1/2)*r*(r*cos(theta)-R__f)*(diff(v__r(r, theta, beta), r))+v__r(r, theta, beta)*(r*cos(theta)-(1/2)*R__f))*((B+2*L__1)*Pi+4*beta*C)^5*phi(theta)^2*(int(phi(theta)^2, theta = 0 .. Pi))-16*(B+2*L__1)^2*B*C*(((B+2*L__1)*Pi+4*beta*C)^2*phi(theta)^2-2*r^2*Pi^2*(B+2*L__1)^2)*r*H*Pi^2)/((int(phi(theta)^2, theta = 0 .. Pi))*((B+2*L__1)*Pi+4*beta*C)^5*phi(theta)^2) = 0

(2)

bc:=v__r(rhs(eq38a),theta,beta)=0;

v__r(((B+2*L__1)*Pi+4*beta*C)*phi(theta)/((B+2*L__1)*Pi), theta, beta) = 0

(3)

sol3:=simplify(dsolve([eq37a,bc], v__r(r, theta,beta)));

v__r(r, theta, beta) = -8*(B+2*L__1)^2*B*(((B+2*L__1)*Pi+4*beta*C)*phi(theta)+r*Pi*(B+2*L__1))*C*r*(((B+2*L__1)*Pi+4*beta*C)*phi(theta)-r*Pi*(B+2*L__1))*H*Pi^2/((r*cos(theta)-R__f)*(Int(phi(theta)^2, theta = 0 .. Pi))*((B+2*L__1)*Pi+4*beta*C)^5*phi(theta)^2)

(4)

#
# Having obtained a solution in terms of phi(theta)
# substitute the original piecewise definition from
# eq41_orig
#
# The result is as ugly as sin!
#
  sol4:=eval( sol3, phi(theta)=rhs(eq41_orig));

sol4 := v__r(r, theta, beta) = -8*(B+2*L__1)^2*B*(((B+2*L__1)*Pi+4*beta*C)*piecewise(theta <= arctan(H/(2*((1/2)*B+L__1))), ((1/2)*B+L__1)/cos(theta), theta <= Pi-arctan(H/(2*((1/2)*B-L__1))), H/(2*sin(theta)), Pi-arctan(H/(2*((1/2)*B-L__1))) < theta, -((1/2)*B-L__1)/cos(theta))+r*Pi*(B+2*L__1))*C*r*(((B+2*L__1)*Pi+4*beta*C)*piecewise(theta <= arctan(H/(2*((1/2)*B+L__1))), ((1/2)*B+L__1)/cos(theta), theta <= Pi-arctan(H/(2*((1/2)*B-L__1))), H/(2*sin(theta)), Pi-arctan(H/(2*((1/2)*B-L__1))) < theta, -((1/2)*B-L__1)/cos(theta))-r*Pi*(B+2*L__1))*H*Pi^2/((r*cos(theta)-R__f)*(Int(piecewise(theta <= arctan(H/(2*((1/2)*B+L__1))), ((1/2)*B+L__1)/cos(theta), theta <= Pi-arctan(H/(2*((1/2)*B-L__1))), H/(2*sin(theta)), Pi-arctan(H/(2*((1/2)*B-L__1))) < theta, -((1/2)*B-L__1)/cos(theta))^2, theta = 0 .. Pi))*((B+2*L__1)*Pi+4*beta*C)^5*piecewise(theta <= arctan(H/(2*((1/2)*B+L__1))), ((1/2)*B+L__1)/cos(theta), theta <= Pi-arctan(H/(2*((1/2)*B-L__1))), H/(2*sin(theta)), Pi-arctan(H/(2*((1/2)*B-L__1))) < theta, -((1/2)*B-L__1)/cos(theta))^2)

(5)

#
# Evaluate this solution for specific B, H, L__1
#
# This result is still ugly
#
  eval(sol4, [B=10, H=10, L__1=4]);

v__r(r, theta, beta) = -(259200*((4*C*beta+18*Pi)*piecewise(theta <= arctan(5/9), 9/cos(theta), theta <= Pi-arctan(5), 5/sin(theta), Pi-arctan(5) < theta, -1/cos(theta))+18*r*Pi))*C*r*((4*C*beta+18*Pi)*piecewise(theta <= arctan(5/9), 9/cos(theta), theta <= Pi-arctan(5), 5/sin(theta), Pi-arctan(5) < theta, -1/cos(theta))-18*r*Pi)*Pi^2/((r*cos(theta)-R__f)*(Int(piecewise(theta <= arctan(5/9), 9/cos(theta), theta <= Pi-arctan(5), 5/sin(theta), Pi-arctan(5) < theta, -1/cos(theta))^2, theta = 0 .. Pi))*(4*C*beta+18*Pi)^5*piecewise(theta <= arctan(5/9), 9/cos(theta), theta <= Pi-arctan(5), 5/sin(theta), Pi-arctan(5) < theta, -1/cos(theta))^2)

(6)

 


 

Download toughODE.mw

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