## 6217 Reputation

10 years, 100 days

## See annotated code attached...

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:

## And...

what further data do you need????

## You say...

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);
 (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]               )         );
 (2)
 >

## Not from a plot command...

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

## This time...

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

## If I understand correctly...

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

 >
 >
 >
 >
 (1)
 >
 (2)
 >
 (3)
 >
 (4)
 >
 (5)
 >
 (6)
 >

## Residuals...

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

## Can't believe...

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

## I should have noted...

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

## The concepts...

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);
 (1)
 >

## It's a version issue...

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                         ):
 (1)
 > # # A few test cases for different grid spacings #   f(4,4);   f(4,6);   f(10,3);
 >

## Can't reproduce...

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().

## It depends on your usage...

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);
 >

## Using Rouben's suggestion...

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)= 2 , H > 0 , B > 0 , L__1 > 0; eq44b:= v__m(beta)= rhs(eq44a) * B * H; eq47:=subs([eq44b,eq38a],eq43);
 (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)));
 (2)
 > bc:=v__r(rhs(eq38a),theta,beta)=0;
 (3)
 > sol3:=simplify(dsolve([eq37a,bc], v__r(r, theta,beta)));
 (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));
 (5)
 > # # Evaluate this solution for specific B, H, L__1 # # This result is still ugly #   eval(sol4, [B=10, H=10, L__1=4]);
 (6)
 >