Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 163 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

A matter of taste?...

tomleslie 13876
September 28 2022
2 1

I personally prefer to have all procedures as "compartmentalised" as possible.

  1. If a procedure use commands from packages Foo and Bar then there should be a "uses Foo, Bar" statement within the procedure. This means that any commands from packages Foo and Bar are available within the procedure -but only within the procedure. Note that you cannot have a with() command inside a procedure, You won't get an error message, but it won't work!. Have a look at the first example in the attached,
  2. Don't use global variables - I will do almost anything to avoid these, just because it is asking for trouble. See examples in the attached

When using package commands at the "top-level", quite a lot depends on how many such commands will be used, how often, and from how many packages. One can either use the "long form" command, as in

PackageName:-PackageCommand()

or the "short form", as in

with(PackageName):
PackageCommand();

or set up one or more "aliases", whihc could be either

alias(PC=PackageName:-PackageCommand):
PC()

or

alias(PN=PackageName):
PN:-PackageCommand()

See the attached for more details.

One thing you should always be aware of when using mutiple packages is that the possibility of naming conflicts increases. As an example both the VectorCalculus() and LinearAlgebra() packages contain a Norm() command, but for the same input, will provide different answers, because (by delault). The former proides the Euclidean norm and the latter will provide the Infinity-norm, If you have loaded both packages using with() commands, then you will get the Norm() command from whichever one happened to be loaded last. (Again see the attached.) In such cases, use of long form package commands, becomes almost mandatory - probably with the second alias construct shown above to minimize subsequent typing.

> 

#
# 'with()' versus 'uses' in a procedure
#
  restart;
  foo1:=  proc(M)
               with(LinearAlgebra);
               return Determinant(M);
          end proc:
  foo1(<<1,2>|<3,4>>);
  foo2:=proc(M)
             uses LinearAlgebra;
             return Determinant(M);
        end proc:
  foo2(<<1,2>|<3,4>>);

Determinant(Matrix(2, 2, {(1, 1) = 1, (1, 2) = 3, (2, 1) = 2, (2, 2) = 4}))

  

-2

 (1)
> 

#
# Fun with globals - case 1.
#
# foo1(5) returns different answers depending
# on whether or not 'B' has been defined
#
  restart:
  foo1:=proc(a)
             return a+B;
        end proc:
  foo1(5);
  B:=3:
  foo1(5);
  B;
#
# Fun with globals - case 2.
#
# Changing a vartiable outside a procedure
# from within a procedure. This is EVIL
#
  restart:
  B:=3:
  foo1:= proc(a)
              global B;
              B:=B+a;
              return a
         end proc:
 foo1(5);
 B;
 

5+B

  

8

  

3

  

5

  

8

 (2)
> 

 

> 

#
# long form package calls, shortform package calls
# and aliases. For ways to achieve the same thing!
#
  restart:
  M:= <<1,2>|<3,4>>:
  LinearAlgebra:-Determinant( M );

  restart:
  M:=<<1,2>|<3,4>>:
  with(LinearAlgebra):
  Determinant(M);

  restart:
  M:=<<1,2>|<3,4>>:
  alias( LAD=LinearAlgebra:-Determinant ):
  LAD(M);

  restart:
  M:=<<1,2>|<3,4>>:
  alias( LA=LinearAlgebra ):
  LA:-Determinant(M);

-2

  

-2

  

-2

  

-2

 (3)
> 

  restart:
#
# This will produce the default Norm() from LinearAlgebra
#
  with(VectorCalculus):
  with(LinearAlgebra):
  Norm(<3,4>);

  restart:
#
# This will produce the default Norm() from VectorCalculus
#
  with(LinearAlgebra):
  with(VectorCalculus):
  Norm(<3,4>);
  restart:
#
# This will produce the default Norm() from VectorCalculus
#
  alias( LA=LinearAlgebra, VC=VectorCalculus ):
  LA:-Norm(<3,4>);
  VC:-Norm(<3,4>);

4

  

5

  

4

  

5

 (4)
> 

 

5

 (5)
> 

 

Download packFun.mw


 

My $0.02...

tomleslie 13876
September 28 2022
0 0

Thhis is largely a matter of "taste", and perhaps(?) a competition between the code brevity (always good because it saves typing!) and code legibility (usually good because it is easier to figure out exactly what is being attempted.

The (perfectly correct) solution provided by Kitonum has gone for code brevity

If one considers code legibility, then the important realisation is that in neither of your examples do you want what would strictly be called an inequality plot, You want an inequality plot and one (or more) other plots. So an alternative is to generate all the plots you want with the appropriate plotting command, and combine them using the display() command

The attached does this for both examples - more typing than Kitonum's solution, but (maybe?) easier to understand  how it works

BTW As is often the case, the plots "render" much better within the Maple worksheet than on this site

> 

  restart:
  with(plots):

  expr1:= x^2+1:
  expr2:= (x-1)^2+(y-1)^2:
  xRange:= -5..8:
  yRange:= -6..6:
#
# First example
#
  display
  ( [ inequal
      ( { y > -1,
          y >= expr1,
          expr2 <= 16
        },
        x = -5 .. 8,
        y = -6 .. 6,
        optionsfeasible = (color = grey)
      ),
      implicitplot
      ( expr2 = 16,
        x = xRange,
        y = yRange,
        color=black
      )
    ],
    scaling=constrained,
    view=[xRange, yRange]
  );
#
# Second example
#
  display
  ( [ plot
      ( expr1,
        x = xRange,
        color=black
      ),
      implicitplot
      ( expr2 = 16,
        x = xRange,
        y = yRange,
        color=black
      ),
      inequal
      ( { y < 2,
          y >= expr1,
          expr2 <= 16
        },
        x = xRange,
        y = yRange,
        optionsfeasible = (color = grey)
      )
    ],
    scaling=constrained,
    view=[xRange, yRange]
  );

 

 
> 

 

Download combPlot.mw

Another way...

tomleslie 13876
September 27 2022
3 3

is shown in the attached - which basically produces a set of equations, such that the entries in the matrix An are equal to the corresponding entries in the identity matrix. Then the solve() command with the option 'allsolutions=true', returns n = 6*_Z1, with _Z1 assumed to be an integer.

See the attached

> 

  restart;
  with(LinearAlgebra):
  A:= < <0 | 1 | 0 | 0>,
        <0 | 0 | 1 | 0>,
        <1 | 0 | 0 | 0>,
        <0 | 0 | 0 | -1>
      >:
  IMat:= IdentityMatrix(4):
  solve( { entries
           ( MatrixPower(A, n)=~ IMat[],
             'nolist'
           )
         },
         allsolutions = true
       )[];
  about(_Z1);

n = 6*_Z1

  

Originally _Z1, renamed _Z1~:
  is assumed to be: integer
 

  
> 

 

> 

 

Download solMat.mw

 

You might find the attached helpful...

tomleslie 13876
September 23 2022
0 1

I have

  1. "cleaned-up" this worksheet extensively
  2. used exact arithmetic (rather than floats) wherever possible
  3. applied Acer;'s suggestion using simplify( fnormal( expr), zero) to remoe tiny imaginary part where necessary
  4. suppressed the result of intermediate calculations - just to make the worksheet more readable
  5. deleted execution groups which appeared to serve no purpose

Thee are still a couple of places where commands do not execute correctly -highlighted in red in the attached. In both instances, this is because the expression beeing evaluated contains an unspecified variable.

So for what it is worth, check out slits.mw

It is as simple...

tomleslie 13876
September 23 2022
1 0

as that shown in the attached.

> 

A := Matrix([[2, -5, 1], [0, 3, -2]]); B := Matrix([[0, 2, -1], [7, 1, 3]]); A+B

Matrix(2, 3, {(1, 1) = 2, (1, 2) = -5, (1, 3) = 1, (2, 1) = 0, (2, 2) = 3, (2, 3) = -2})

  

Matrix(2, 3, {(1, 1) = 0, (1, 2) = 2, (1, 3) = -1, (2, 1) = 7, (2, 2) = 1, (2, 3) = 3})

  

Matrix(%id = 36893488148485795220)

 (1)
> 

``

Download addMat.mw

Well I'd...

tomleslie 13876
September 23 2022
1 3

I'd have to dig out a textbook to find out what the "differential transform method" actually is - but since Maple's built-in numeric solvers have absolutely no problem producing a solution, I can't see why I would bother.

In the attached I have assumed that "infinity"=~10 and M=1. The former is probably reasonable since nothing much seems to change when eta>~5. The value of 'M' is a complete guess - so feel free to adjust entries in params to anything you want

> 

  restart:
  with(plots):
  odes:=  diff( f(eta), eta$3)+f(eta)*diff(f(eta), eta$2)-diff(f(eta),eta)^2-M*diff(f(eta),eta)=0,
          diff(g(eta),eta)-2*diff(f(eta),eta$2)-2*f(eta)*diff(f(eta),eta)=0;
  params:=[M=1, inf=10]:
  conds:= f(0)=0, g(0)=0, D(f)(0)=1, D(f)(inf)=0;
  sol:=dsolve(eval([odes, conds], params), numeric);

diff(diff(diff(f(eta), eta), eta), eta)+f(eta)*(diff(diff(f(eta), eta), eta))-(diff(f(eta), eta))^2-M*(diff(f(eta), eta)) = 0, diff(g(eta), eta)-2*(diff(diff(f(eta), eta), eta))-2*f(eta)*(diff(f(eta), eta)) = 0

  

f(0) = 0, g(0) = 0, (D(f))(0) = 1, (D(f))(inf) = 0

  

proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(25, {(1) = .0, (2) = .3559615719620481, (3) = .7177327124223771, (4) = 1.0856075113942896, (5) = 1.459919452367661, (6) = 1.8434216675992934, (7) = 2.238706437897399, (8) = 2.646931490780542, (9) = 3.0683362454300322, (10) = 3.4972043196777696, (11) = 3.9327507094693015, (12) = 4.375167527405896, (13) = 4.820471658546179, (14) = 5.266655003597093, (15) = 5.71372279100709, (16) = 6.161284885670424, (17) = 6.60894569614419, (18) = 7.056705287149855, (19) = 7.504531478794849, (20) = 7.952359501530675, (21) = 8.40018752606533, (22) = 8.841187217869862, (23) = 9.252338197962562, (24) = 9.637106306133425, (25) = 10.0}, datatype = float[8], order = C_order); Y := Matrix(25, 4, {(1, 1) = .0, (1, 2) = 1.0, (1, 3) = -1.4142135629642079, (1, 4) = .0, (2, 1) = .27968316962547096, (2, 2) = .6044682680758878, (2, 3) = -.8548472237007583, (2, 4) = -.7128409446596393, (3, 1) = .45085673745298804, (3, 2) = .36239228649394495, (3, 3) = -.5125000879572456, (3, 4) = -1.071943661977726, (4, 1) = .5547994909109478, (4, 2) = .21539503416613714, (4, 3) = -.30461458079948894, (4, 4) = -1.2614074290724118, (5, 1) = .6174002801908669, (5, 2) = .12686414795475193, (5, 3) = -.1794130017647044, (5, 4) = -1.3650885619743758, (6, 1) = .6549533736832511, (6, 2) = 0.7375605239653821e-1, (6, 3) = -.1043068139775268, (6, 4) = -1.4235239451709583, (7, 1) = .6772869712733923, (7, 2) = 0.4217157375297486e-1, (7, 3) = -0.59639617543090266e-1, (7, 4) = -1.4569391915469085, (8, 1) = .6903658572625241, (8, 2) = 0.2367523293262649e-1, (8, 3) = -0.3348184368945545e-1, (8, 4) = -1.4760445037561656, (9, 1) = .6978819304430172, (9, 2) = 0.13045896256020467e-1, (9, 3) = -0.18449694594204938e-1, (9, 4) = -1.4868690084099765, (10, 1) = .7020769235506761, (10, 2) = 0.7113274484681226e-2, (10, 3) = -0.10059704503693324e-1, (10, 4) = -1.4928614357065986, (11, 1) = .7043900229785282, (11, 2) = 0.38420501016844626e-2, (11, 3) = -0.5433500210849747e-2, (11, 4) = -1.4961505871925886, (12, 1) = .7056535711498219, (12, 2) = 0.2055112295195886e-2, (12, 3) = -0.29063962594502897e-2, (12, 4) = -1.4979428050269092, (13, 1) = .7063326093297881, (13, 2) = 0.10947923215180097e-2, (13, 3) = -0.15483093592918544e-2, (13, 4) = -1.4989046525226157, (14, 1) = .7066948548145358, (14, 2) = 0.5824792655130951e-3, (14, 3) = -0.8238038721453331e-3, (14, 4) = -1.499417415837864, (15, 1) = .706887857019484, (15, 2) = 0.3095046369407773e-3, (15, 3) = -0.4377794810384208e-3, (15, 4) = -1.4996905405205363, (16, 1) = .7069904874251435, (16, 2) = 0.16432445173196787e-3, (16, 3) = -0.23249120015165323e-3, (16, 4) = -1.4998357939841824, (17, 1) = .7070449797062479, (17, 2) = 0.8720736096277632e-4, (17, 3) = -0.1234689145759178e-3, (17, 4) = -1.4999129741477275, (18, 1) = .7070738956722984, (18, 2) = 0.4624074364800964e-4, (18, 3) = -0.6558518539263504e-4, (18, 4) = -1.4999540167689362, (19, 1) = .7070892189395271, (19, 2) = 0.2446979348149444e-4, (19, 3) = -0.34867538809472405e-4, (19, 4) = -1.4999758890699202, (20, 1) = .7070973129819845, (20, 2) = 0.1288503758192738e-4, (20, 3) = -0.18581842110460757e-4, (20, 4) = -1.4999876120958646, (21, 1) = .7071015546799427, (21, 2) = 0.6696886913585484e-5, (21, 3) = -0.9964464655586802e-5, (21, 4) = -1.4999939897927408, (22, 1) = .7071037080005216, (22, 2) = 0.33961704724957093e-5, (22, 3) = -0.54771878363046635e-5, (22, 4) = -1.499997545988323, (23, 1) = .7071047142936876, (23, 2) = 0.16513169257183738e-5, (23, 3) = -0.323709878419725e-5, (23, 4) = -1.4999996125871435, (24, 1) = .7071051424928092, (24, 2) = 0.6469795665018257e-6, (24, 3) = -0.20987272394271767e-5, (24, 4) = -1.5000010156984427, (25, 1) = .7071052536829966, (25, 2) = .0, (25, 3) = -0.15300543789426414e-5, (25, 4) = -1.5000021524112566}, datatype = float[8], order = C_order); YP := Matrix(25, 4, {(1, 1) = 1.0, (1, 2) = -1.4142135629642079, (1, 3) = 2.0, (1, 4) = -2.8284271259284157, (2, 1) = .6044682680758878, (2, 2) = -.8548472237007583, (2, 3) = 1.2089365362567133, (2, 4) = -1.3715752450945502, (3, 1) = .36239228649394495, (3, 2) = -.5125000879572456, (3, 3) = .7247845734050276, (3, 4) = -.698226167980914, (4, 1) = .21539503416613714, (4, 2) = -.30461458079948894, (4, 3) = .4307900692611768, (4, 4) = -.37022705099873965, (5, 1) = .12686414795475193, (5, 2) = -.1794130017647044, (5, 3) = .2537282975504501, (5, 4) = -.20217408254252994, (6, 1) = 0.7375605239653821e-1, (6, 2) = -.1043068139775268, (6, 3) = .14751210737439158, (6, 4) = -.1120000772617109, (7, 1) = 0.4217157375297486e-1, (7, 2) = -0.59639617543090266e-1, (7, 3) = 0.8434315131944053e-1, (7, 4) = -0.6215472016421087e-1, (8, 1) = 0.2367523293262649e-1, (8, 2) = -0.3348184368945545e-1, (8, 3) = 0.47350471308441354e-1, (8, 4) = -0.3427454242006564e-1, (9, 1) = 0.13045896256020467e-1, (9, 2) = -0.18449694594204938e-1, (9, 3) = 0.26091800144631153e-1, (9, 4) = -0.1869039866138809e-1, (10, 1) = 0.7113274484681226e-2, (10, 2) = -0.10059704503693324e-1, (10, 3) = 0.14226559548357533e-1, (10, 4) = -0.10131277274233616e-1, (11, 1) = 0.38420501016844626e-2, (11, 2) = -0.5433500210849747e-2, (11, 3) = 0.7684114789042607e-2, (11, 4) = -0.5454396902879145e-2, (12, 1) = 0.2055112295195886e-2, (12, 2) = -0.29063962594502897e-2, (12, 3) = 0.4110244681399333e-2, (12, 4) = -0.2912397858462812e-2, (13, 1) = 0.10947923215180097e-2, (13, 2) = -0.15483093592918544e-2, (13, 3) = 0.21896122815436127e-2, (13, 4) = -0.15500436843196444e-2, (14, 1) = 0.5824792655130951e-3, (14, 2) = -0.8238038721453331e-3, (14, 3) = 0.11649965054292464e-2, (14, 4) = -0.8243375443421578e-3, (15, 1) = 0.3095046369407773e-3, (15, 2) = -0.4377794810384208e-3, (15, 3) = 0.6190614292594163e-3, (15, 4) = -0.43798882298752264e-3, (16, 1) = 0.16432445173196787e-3, (16, 2) = -0.23249120015165323e-3, (16, 3) = 0.3287205211746788e-3, (16, 4) = -0.23263075185159962e-3, (17, 1) = 0.8720736096277632e-4, (17, 2) = -0.1234689145759178e-3, (17, 3) = 0.17451304228726467e-3, (17, 4) = -0.12361877562751237e-3, (18, 1) = 0.4624074364800964e-4, (18, 2) = -0.6558518539263504e-4, (18, 3) = 0.9261645438834314e-4, (18, 4) = -0.6577912528530556e-4, (19, 1) = 0.2446979348149444e-4, (19, 2) = -0.34867538809472405e-4, (19, 3) = 0.4912485303542096e-4, (19, 4) = -0.3513042329806194e-4, (20, 1) = 0.1288503758192738e-4, (20, 2) = -0.18581842110460757e-4, (20, 3) = 0.2602437423268316e-4, (20, 4) = -0.1894173331721604e-4, (21, 1) = 0.6696886913585484e-5, (21, 2) = -0.9964464655586802e-5, (21, 3) = 0.13742820211398585e-4, (21, 4) = -0.10458171014949486e-4, (22, 1) = 0.33961704724957093e-5, (22, 2) = -0.54771878363046635e-5, (22, 3) = 0.7269121834935969e-5, (22, 4) = -0.6151486204402128e-5, (23, 1) = 0.16513169257183738e-5, (23, 2) = -0.323709878419725e-5, (23, 3) = 0.3940287463506203e-5, (23, 4) = -0.41388896024576575e-5, (24, 1) = 0.6469795665018257e-6, (24, 2) = -0.20987272394271767e-5, (24, 3) = 0.2131000808773079e-5, (24, 4) = -0.32824893217319343e-5, (25, 1) = .0, (25, 2) = -0.15300543789426414e-5, (25, 3) = 0.10819094897710164e-5, (25, 4) = -0.3060108757885283e-5}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(25, {(1) = .0, (2) = .3559615719620481, (3) = .7177327124223771, (4) = 1.0856075113942896, (5) = 1.459919452367661, (6) = 1.8434216675992934, (7) = 2.238706437897399, (8) = 2.646931490780542, (9) = 3.0683362454300322, (10) = 3.4972043196777696, (11) = 3.9327507094693015, (12) = 4.375167527405896, (13) = 4.820471658546179, (14) = 5.266655003597093, (15) = 5.71372279100709, (16) = 6.161284885670424, (17) = 6.60894569614419, (18) = 7.056705287149855, (19) = 7.504531478794849, (20) = 7.952359501530675, (21) = 8.40018752606533, (22) = 8.841187217869862, (23) = 9.252338197962562, (24) = 9.637106306133425, (25) = 10.0}, datatype = float[8], order = C_order); Y := Matrix(25, 4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = -0.10082930375400029e-12, (1, 4) = .0, (2, 1) = -0.18004877667327648e-6, (2, 2) = 0.2546273850475269e-6, (2, 3) = -0.36009766321709544e-6, (2, 4) = 0.5650274213258057e-6, (3, 1) = -0.16755861718865572e-6, (3, 2) = 0.23696356800749333e-6, (3, 3) = -0.335117343874476e-6, (3, 4) = 0.3554712125671065e-6, (4, 1) = -0.10666884845538775e-6, (4, 2) = 0.1508523286330195e-6, (4, 3) = -0.21333778623288643e-6, (4, 4) = 0.15574220089535004e-6, (5, 1) = -0.50934146337702274e-7, (5, 2) = 0.720314005294041e-7, (5, 3) = -0.1018683264414776e-6, (5, 4) = 0.4494427229892633e-7, (6, 1) = -0.14196711505746563e-7, (6, 2) = 0.20076592932927846e-7, (6, 3) = -0.2839334998346839e-7, (6, 4) = -0.682693059827394e-8, (7, 1) = 0.52980878297142255e-8, (7, 2) = -0.7493544342990351e-8, (7, 3) = 0.10596426925955988e-7, (7, 4) = -0.27324089170223783e-7, (8, 1) = 0.12894605666859148e-7, (8, 2) = -0.18237104115747125e-7, (8, 3) = 0.25789738369894173e-7, (8, 4) = -0.3218921034388213e-7, (9, 1) = 0.13573326130629444e-7, (9, 2) = -0.19197603897477637e-7, (9, 3) = 0.271475895967633e-7, (9, 4) = -0.2970171155113244e-7, (10, 1) = 0.10985453216972764e-7, (10, 2) = -0.15538679029262038e-7, (10, 3) = 0.2197242582051388e-7, (10, 4) = -0.2440706556754974e-7, (11, 1) = 0.748612170343319e-8, (11, 2) = -0.10591073817527832e-7, (11, 3) = 0.14974574737952034e-7, (11, 4) = -0.1878253252982445e-7, (12, 1) = 0.43460618229992916e-8, (12, 2) = -0.6151981619329957e-8, (12, 3) = 0.8695573339613013e-8, (12, 4) = -0.14091553230151405e-7, (13, 1) = 0.20382412027984324e-8, (13, 2) = -0.28904123820353432e-8, (13, 3) = 0.40814508341588015e-8, (13, 4) = -0.10748563475531859e-7, (14, 1) = 0.5834613678934604e-9, (14, 2) = -0.835960182357947e-9, (14, 3) = 0.11739344672608685e-8, (14, 4) = -0.8672561581270383e-8, (15, 1) = -0.20169921482630112e-9, (15, 2) = 0.2705318613900875e-9, (15, 3) = -0.39365050611876587e-9, (15, 4) = -0.756366709532679e-8, (16, 1) = -0.5406482754895094e-9, (16, 2) = 0.7446912271411817e-9, (16, 3) = -0.10678939987483363e-8, (16, 4) = -0.7093408037075742e-8, (17, 1) = -0.6207517363351042e-9, (17, 2) = 0.8510747083508336e-9, (17, 3) = -0.12232324636803393e-8, (17, 4) = -0.6993808618976654e-8, (18, 1) = -0.5730122617879785e-9, (18, 2) = 0.7743891642501881e-9, (18, 3) = -0.11212757689915765e-8, (18, 4) = -0.70797360779954445e-8, (19, 1) = -0.477489776065208e-9, (19, 2) = 0.6271208877409757e-9, (19, 3) = -0.9216253647046788e-9, (19, 4) = -0.72392431617943276e-8, (20, 1) = -0.37744660114567547e-9, (20, 2) = 0.4694787617742575e-9, (20, 3) = -0.710117888319266e-9, (20, 4) = -0.74130781630948035e-8, (21, 1) = -0.2929420815501417e-9, (21, 2) = 0.3285438550398238e-9, (21, 3) = -0.5259613033791466e-9, (21, 4) = -0.7575455294801684e-8, (22, 1) = -0.23101825200587323e-9, (22, 2) = 0.21297015519175516e-9, (22, 3) = -0.38231827511526684e-9, (22, 4) = -0.7719036748239143e-8, (23, 1) = -0.1921976190896188e-9, (23, 2) = 0.12339770146716387e-9, (23, 3) = -0.28016240781922607e-9, (23, 4) = -0.7843284231519336e-8, (24, 1) = -0.17352042088890837e-9, (24, 2) = 0.5460541546875159e-10, (24, 3) = -0.21284296187236342e-9, (24, 4) = -0.7954455402860358e-8, (25, 1) = -0.17125925853219432e-9, (25, 2) = .0, (25, 3) = -0.1719713463799893e-9, (25, 4) = -0.8060467975552293e-8}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[25] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(5.650274213258057e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [4, 25, [f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta), g(eta)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[25] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[25] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(4, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(25, 4, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(4, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(25, 4, X, Y, outpoint, yout, L, V) end if; [eta = outpoint, seq('[f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta), g(eta)]'[i] = yout[i], i = 1 .. 4)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[25] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(5.650274213258057e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [4, 25, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[25] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[25] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(25, 4, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(4, {(1) = 0., (2) = 0., (3) = 0., (4) = 0.}); `dsolve/numeric/hermite`(25, 4, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 4)] end proc, (2) = Array(0..0, {}), (3) = [eta, f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta), g(eta)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [eta = res[1], seq('[f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta), g(eta)]'[i] = res[i+1], i = 1 .. 4)] catch: error end try end proc

 (1)
> 

  odeplot( sol,
           [ [eta, f(eta)],
             [eta, g(eta)],
             [eta, diff(f(eta), eta)]
           ],
           eta=0..10,
           color=[black, red, blue],
           legend=[f(eta), g(eta), diff(f(eta),eta)]
         );

 
> 

  
  
> 

 

Download odeProb.mw

An alternative...

tomleslie 13876
September 21 2022
1 0

Since the problem of reflection through a point has now occurred a couple of times, I have come up with a more "general" answer, whihc is to overload the geometry:-reflection() command with one which actually works!

See the additions

  myRefl:= module()
                  option package;
                  export reflection;
                  reflection:= proc( obj1, obj2, ptOrl)
                                     if   geometry:-form(ptOrl)=point2d
                                     then geometry:-rotation( obj1, obj2, Pi, clockwise, ptOrl):
                                     else geometry:-reflection( obj1, obj2, ptOrl):
                                     end if:
                               end proc:
           end module:
  with(myRefl):

to your original code, which now "works" - see the attached. Note that, in order for the overloaded command to take effect, it is important that the definition of myRefl module, and with(myRef1) command occurs after the with(geometry) command.

> 

  restart;
  with(plots):
  with(geometry):
  _EnvHorizontalName := x:
  _EnvVerticalName := y:
  myRefl:= module()
                  option package;
                  export reflection;
                  reflection:= proc( obj1, obj2, ptOrl)
                                     if   geometry:-form(ptOrl)=point2d
                                     then geometry:-rotation( obj1, obj2, Pi, clockwise, ptOrl):
                                     else geometry:-reflection( obj1, obj2, ptOrl):
                                     end if:
                               end proc:
           end module:
  with(myRefl):
  R := 11:
  r := 7:
  a := sqrt(R*r):
  b := 2:
  circle(C1, [point(P1, [0, 0]), R]):
  circle(C2, [point(P2, [R + 2*b + r, 0]), r]):
  ellipse(p, (x - R - b)^2/b^2 + y^2/a^2 = 1):
  draw([C1(color = yellow, filled = true),
        C2(color = red, filled = true), p(color = blue, filled = true),
        C1(color = black), C2(color = black), p(color = black)],
        axes = none, view = [-15 .. 35, -15 .. 15], scaling = constrained):
  alpha := arctan((R - r)/(R + 2*b + r));
  long := cos(alpha)*(R + 2*b + r);
  evalf(%);
  circle(C2, [point(P2, [long, r - R]), r]);
  rotation(p1, p, alpha, 'clockwise');
  detail(p1);
  point(A, 0, -R);
  point(B, long, -R);
  line(L1, [A, B]);
  point(cen, [(143*sqrt(5))/25, -(26*sqrt(5))/25]);
  reflection(L2, L1, cen);
  detail(L2);
  draw([C1(color = yellow, filled = true),
        C2(color = red, filled = true),
        p1(color = blue, filled = true),
        C1(color = black),
        C2(color = black),
        p1(color = black),
        L1(color = black),
        L2(color = black)
       ],
        axes = none,
        view = [-15 .. 35, -15 .. 15],
        scaling = constrained);

alpha := arctan(2/11)

  

long := 242*sqrt(5)*(1/25)

  

21.64513802

  

C2

  

p1

  

GeometryDetail(["name of the object", p1], ["form of the object", ellipse2d], ["center", [(143/25)*sqrt(5), -(26/25)*sqrt(5)]], ["foci", [[(143/25)*sqrt(5)-(2/25)*sqrt(73)*sqrt(5), -(26/25)*sqrt(5)-(11/25)*sqrt(73)*sqrt(5)], [(143/25)*sqrt(5)+(2/25)*sqrt(73)*sqrt(5), -(26/25)*sqrt(5)+(11/25)*sqrt(73)*sqrt(5)]]], ["length of the major axis", 2*sqrt(77)], ["length of the minor axis", 4], ["equation of the ellipse", (18/875)*y^2-(73/875)*x*y+(9333/38500)*x^2+(13/25)*y*sqrt(5)-(143/50)*x*sqrt(5)+165/4 = 0])

  

A

  

B

  

L1

  

cen

  

L2

  

GeometryDetail(["name of the object", L2], ["form of the object", line2d], ["equation of the line", -12584/125-(242/25)*y*5^(1/2)+(2662/25)*5^(1/2) = 0])

  

 
> 

 

Download geoRfl2.mw

 

I'm pretty sure...

tomleslie 13876
September 21 2022
0 0

that this is another manifestation of the "bug" you found previously - namely that reflections of geometric objects in a point seems to be very error-prone.

I can't remember if I informed you, but I did reportt this bug to Maple, they have admitted that there is an issue and hopefully it will be addressed.

Meanwhile I think the workaround in the attached will achieve what you want

> 

  restart:
  with(plots):
  with(geometry):
  _EnvHorizontalName := x:
  _EnvVerticalName := y:
  R := 11:
  r := 7:
  a := sqrt(R*r):
  b := 2:
  circle(C1, [point(P1, [0, 0]), R]):
  circle(C2, [point(P2, [R + 2*b + r, 0]), r]):
  ellipse(p, (x - R - b)^2/b^2 + y^2/a^2 = 1):
  alpha := arctan((R - r)/(R + 2*b + r)):
  long := cos(alpha)*(R + 2*b + r):
  circle(C2, [point(P2, [long, r - R]), r]):
  rotation(p1, p, alpha, 'clockwise'):
  point(A, 0, -R):
  point(B, long, -R):
  line(L1, [A, B]):  
  reflection( L2, L1, line(LZ, y=coordinates(center(p1))[2])):
  draw([ C1(color = yellow, filled = true),
         C2(color = red, filled = true),
         p1(color = blue, filled = true),
         C1(color = black),
         C2(color = black),
         p1(color =   black),
         L1(color = black),
         L2(color = black)
       ],
       axes = none,
       view = [-15 .. 35, -15 .. 15],
       scaling = constrained
     );

 
> 

 

Download geoRefl.mw

An even simpler way...

tomleslie 13876
September 21 2022
3 0

is shown in the attached

> 

  restart;
  A := {la, ra, rb, t, t10, t11, t12, t13, t15, t17, t18, t19, t21, t23, t24, t25, t26, t27, t28, t31, t32, t33, t34, t35, t36, t37, t39, t41, t42, t43, t44,   t45, t46, t48, t49, t5, t50, t52, t53, t54, t55, t56, t57, t58, t59, t61, t62, t63, t64, t65, t67, t68, t69, t7, t70, t9, xa, za}:
indets(A, suffixed(t, nonnegint)) ;

{t10, t11, t12, t13, t15, t17, t18, t19, t21, t23, t24, t25, t26, t27, t28, t31, t32, t33, t34, t35, t36, t37, t39, t41, t42, t43, t44, t45, t46, t48, t49, t5, t50, t52, t53, t54, t55, t56, t57, t58, t59, t61, t62, t63, t64, t65, t67, t68, t69, t7, t70, t9}

 (1)

NULL

Download tSuff.mw

You cannot use ShowSolution()...

tomleslie 13876
September 20 2022
2 1

to show the steps for any Maple calculation. It only works for a limited number of "model" problems within the StudentCalculus() package. For yu specific problem, things may be a little clearer if you add a simplify() command - this cerytainly provides a more understandable output expression. See the atached

> 

with(VectorCalculus)
> 

SetCoordinates(spherical[r, theta, phi])

spherical[r, theta, phi]

 (1)
> 

field1 := VectorField(`<,>`(0, f(r, theta, phi), 0))

Vector(3, {(1) = 0, (2) = f(r, theta, phi), (3) = 0})

 (2)
> 

pde := -(diff(field1, theta)) = 0

(Vector(3, {(1) = -(diff(f(r, theta, phi), theta))*cos(phi)^2*cos(theta)*sin(theta)+f(r, theta, phi)*cos(phi)^2*sin(theta)^2-(diff(f(r, theta, phi), theta))*cos(theta)*sin(phi)^2*sin(theta)+f(r, theta, phi)*sin(theta)^2*sin(phi)^2+(diff(f(r, theta, phi), theta))*sin(theta)*cos(theta)+f(r, theta, phi)*cos(theta)^2, (2) = -(diff(f(r, theta, phi), theta))*cos(phi)^2*cos(theta)^2+f(r, theta, phi)*cos(phi)^2*sin(theta)*cos(theta)-(diff(f(r, theta, phi), theta))*cos(theta)^2*sin(phi)^2+f(r, theta, phi)*sin(theta)*sin(phi)^2*cos(theta)-(diff(f(r, theta, phi), theta))*sin(theta)^2-f(r, theta, phi)*cos(theta)*sin(theta), (3) = 0})) = 0

 (3)
> 

simplify(%)

(Vector(3, {(1) = f(r, theta, phi), (2) = -(diff(f(r, theta, phi), theta)), (3) = 0})) = 0

 (4)
> 

``

Download Vdiff.mw

Another way...

tomleslie 13876
September 16 2022
1 0

is shown in the attached

> 

  restart;
  p:=x^6+2*x^3+x^2;
  trunP:= (poly, n)->select( j->degree(j,indets(p, name)[])<=n, poly):
  trunP(p,4);
  trunP(p,2);

x^6+2*x^3+x^2

  

2*x^3+x^2

  

x^2

 (1)
> 

 

Download trPol.mw

Do you need Maple for this?...

tomleslie 13876
September 16 2022
0 0

If you have

ode:=diff(y(x),x)=2-y(x);

and

initial condition: (dy/dx)(0)=0

then from the ode definition, at x=0 it must be true that y(0)=2. It is also fairly obvious that y(x)=2, satisfies the ode for all values of the independent variable 'x'.

This is the solution which Maple produces - see the attached

> 

 restart;
 ode:=diff(y(x),x)=2-y(x);
 init:=D(y)(0)=0;
 dsolve([ode,init]);

diff(y(x), x) = 2-y(x)

  

(D(y))(0) = 0

  

y(x) = 2

 (1)
> 

 

Download ODEtriv.mw

In the attached...

tomleslie 13876
September 14 2022
0 4

I have restricted myself to only two values in the command

ans1:= Array([seq( [j, doCalc(j)], j=0..0.1, 0.1)]):

Data from the three plots for each of the above two values is then written to a single matrix, which is then exported to Excel, or in a format sutable for reading directly by Mallab, using its readmatrix() command. Obviously you will have to change the file paths to something appropriate for your machine.

> 

restart:
> 

#kernelopts(version):

Procedure

 

``

> 

doCalc:= proc( xi )

                 # Import Packages
                 uses ArrayTools, Student:-Calculus1, LinearAlgebra,
                      ListTools, RootFinding, plots, ListTools:
                 local gamma__1:= .1093,
                       alpha__3:= -0.1104e-2,
                       k__1:= 6*10^(-12),
                       d:= 0.2e-3,
                       theta0:= 0.0001,
                       eta__1:= 0.240e-1,
                       alpha:= 1-alpha__3^2/(gamma__1*eta__1),
                       c:= alpha__3*xi*alpha/(eta__1*(4*k__1*q^2/d^2-alpha__3*xi/eta__1)),
                       theta_init:= theta0*sin(Pi*z/d),
                       n:= 30,
                       g, f, b1, b2, qstar, OddAsymptotes, ModifiedOddAsym,
                       qstarTemporary, indexOfqstar2, qstar2, AreThereComplexRoots,
                       soln1, soln2, qcomplex1, qcomplex2, gg, qq, m, pp, j, i,
                       AllAsymptotes, p, Efun, b, aa, F, A, B, Ainv, r, theta_sol, v, Vfun, v_sol,                           minp,nstar, soln3, Imagroots, Imagroot2, AreTherePurelyImaginaryRoots, k;

# Assign g for q and plot g, Set q as a complex and Compute the Special Asymptotes
AreThereComplexRoots := type(FAIL, 'truefalse');
g:= q-(1-alpha)*tan(q)+ c*tan(q):
f:= subs(q = x+I*y, g):
b1:= evalc(Re(f)) = 0:
b2:=y-(1-alpha)*tanh(y)-(alpha__3*xi*alpha/(eta__1*(4*k__1*y^2/d^2+alpha__3*xi/eta__1)))*tanh(y) = 0:
qstar:= (fsolve(1/c = 0, q = 0 .. infinity)):
OddAsymptotes:= Vector[row]([seq(evalf(1/2*(2*j + 1)*Pi), j = 0 .. n)]);

# Compute Odd asymptote

 ModifiedOddAsym:= abs(`-`~(OddAsymptotes, qstar));
 qstarTemporary:= min(ModifiedOddAsym);
 indexOfqstar2:= SearchAll(qstarTemporary, ModifiedOddAsym);
 qstar2:= OddAsymptotes(indexOfqstar2);

#compute x and y
AreThereComplexRoots:= type(true, 'truefalse');
  try
soln1:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = 0 .. infinity});
soln2:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = -infinity .. 0});
qcomplex1:= subs(soln1, x+I*y);
qcomplex2:= subs(soln2, x+I*y);
  catch:
AreThereComplexRoots:= type(FAIL, 'truefalse');
  end try;

AreTherePurelyImaginaryRoots:= type(true, 'truefalse'):
  try
Imagroots:= remove( `=`, I*~(Roots(b2, y)),0);
  catch:
AreTherePurelyImaginaryRoots:= type(FAIL, 'truefalse');
  end try;

## Compute the rest of the Roots
if (xi > 0) then
OddAsymptotes := Vector[row]([seq(evalf((1/2)*(2*j+1)*Pi), j = 0 .. n)]);
AllAsymptotes := sort(Vector[row]([OddAsymptotes, qstar]));
if AreThereComplexRoots then
gg := [qcomplex1, qcomplex2, op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])), seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
elif not AreThereComplexRoots then
gg := [op(Imagroots[1]), op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])), seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
end if;``NULL
else
OddAsymptotes := Vector[row]([seq(evalf((1/2)*(2*j+1)*Pi), j = 0 .. n)]);
AllAsymptotes := sort(Vector[row]([OddAsymptotes, qstar]));
if AreThereComplexRoots then
gg := [qcomplex1, qcomplex2, op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])), seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
elif not AreThereComplexRoots then
gg := [op(Imagroots), op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])), seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
end if;NULL``
end if:


# Remove the repeated roots if any & Redefine n

  qq:= MakeUnique(gg):
  m:= numelems(qq):

## Compute the `&tau;_n`time constants


 p:= [seq](gamma__1*alpha/(4*k__1*qq[i]^2/d^2-alpha__3*xi/eta__1), i= 1..m):

   
## Compute θ_n from initial conditions

   for i to m do
   Efun[i] := cos(qq[i]-2*qq[i]*z/d)-cos(qq[i]):
   end do:

## Compute integral coefficients for off-diagonal elements θ_n matrix

                printlevel := 2:
                for i to m do
                    for j from i+1 to m do
                        b[i, j] := int(Efun[i]*Efun[j], z = 0 .. d):
                        b[j, i] := b[i, j]:
                        aa[i, j] := b[i, j]:
                        aa[j, i] := b[j, i]:
                    end do :
                end do:

## Calculate integral coefficients for diagonal elements in theta_n matrix

  for i to m do
  aa[i, i] := int(Efun[i]^2, z = 0 .. d):
  end do:

## Calculate integrals for RHS vector

   for i to m do
   F[i] := int(theta_init*Efun[i], z = 0 .. d):
   end do:

## Create matrix A and RHS vector B

  A := Matrix([seq([seq(aa[i, j], i = 1 .. m)], j = 1 .. m)]):
  B := Vector([seq(F[i], i = 1 .. m)]):

## Calculate inverse of A and solve A*x=B
#Ainv := 1/A #r := MatrixVectorMultiply(Ainv, B)

              r := LinearSolve(A,B):

## Define Theta(z,t)
             theta_sol := add(r[i]*Efun[i]*exp(-t/p[i]), i = 1 .. m):

## Compute v_n for n times constant

             for i to m do
             v[i] := (-2*k__1*alpha__3*qq[i])*(1/(d^2*eta__1*alpha*gamma__1))+ alpha__3^2*xi/(2*(eta__1)^2*qq[i]*alpha*gamma__1)+xi/(2*eta__1*qq[i]):
             end do:

## Compute v(z,t) from initial conditions
             for i to m do
             Vfun[i] := d*sin(qq[i]-2*qq[i]*z/d)+(2*z-d)*sin(qq[i]):
             end do:

## Define v(z,t)
             v_sol := add(v[i]*r[i]*Vfun[i]*exp(-t/p[i]), i = 1 .. m):

##
             minp:=min( seq( Re(p[i]), i=1..m) ):
             member(min(seq( Re(p[i]), i=1..m)),[seq( Re(p[i]), i=1..m)],'nstar'):

## Return all the plots
            return theta_init,plot( [seq( Re(eval(theta_sol, t=j)), j in [1,3,5] )], z = 0 .. d), plot( [ seq( Re(eval( v_sol, t=j)), j in [1,3,5] ) ],z=0..d),theta_sol, v_sol, g, qq, Imagroot1, Imagroot2, minp, p[nstar], evalf(p);
  end proc:
> 

#

> 

#

Run the calculation for supplied value of 'xi

 

NULL

``

> 

ans:=[doCalc(-0.06)]:

r(q) graph

> 

with(plots):
plot(ans[6], q=0..evalf((5*Pi)),view=[DEFAULT,-20..20], labels =[q, r(q)], labelfont = ["TimesNewRoman", 14], labeldirections = ["horizontal", "vertical"]);

 
> 

with(plots): d:=0.0002:

Director angle plot for at d/2 as a funtion of time

> 

theta_f := subs(z=d/2, ans[4]):
plot(theta_f, t=0..10, view=[DEFAULT, DEFAULT], labels =[typeset(t," (seconds)"), typeset(theta('d/2', t),"(rad)")], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14], labeldirections = ["horizontal", "vertical"], axes=boxed);

 
> 

#extract theta_f as excel file
#ExcelTools:-Export(op([1, 1], plot(theta_f, t=0..100)));

Plot for many xi values

 
> 

#Convert and store xi and doCalc into an array
# ans1:= Array([seq( [j, doCalc(j)], j=-1..1, 0.1)]):
ans1:= Array([seq( [j, doCalc(j)], j=0..0.1, 0.1)]):

> 

#plots:-display( convert( ans1[8..10,3], list), axes=boxed, axis[1]=[tickmarks=[0="0", 0.00005="50", 0.00010="100", 0.00015=#"150", 0.00020="200"]], #axis[2]=  [tickmarks=[0="0", 0.00002="2e-05", 0.00004="4e-05", 0.00006="6e-05", 0.00008= "8e-05", #0.0001="1e-04"]], axes=boxed, labels =[typeset(z,#" (`&mu;m`)"), typeset(theta(z, t)," (rad)")], labelfont =["TimesNewRoman", 14], #axesfont = [14, 14]);
#plots:-display( convert( ans1[12..14,3], list), axes=boxed);
#plots:-display( convert( ans1[8..10,4], list), axes=boxed);
#plots:-display( convert( ans1[12..14,4], list), axes=boxed);

> 

#ExcelTools:-Export(op([1, 1], plots:-display( convert( ans1[12..14,3], list))));

> 

myData:= `<|>`( seq( seq( plottools:-getdata(ans1[i,3])[j,3],j=1..3), i=1..2)
              );

_rtable[36893488148134438060]

(3.1)
> 

 #
 # Export the above myData to Excel
 #
   ExcelTools:-Export( myData, "C:/Users/TomLeslie/Desktop/dat1.xlsx");

> 

 #
 # Export the above myData to Matlab - this can be read from within
 # Matlab by using the command
 #
 #             M=readmatrix('C:/Users/TomLeslie/Desktop/dat2')
 #
   ExportMatrix("C:/Users/TomLeslie/Desktop/dat2", myData, target=Matlab, mode=ascii);

54139

(3.2)

NULL

Download dataOrg.mw

Keep it simple...

tomleslie 13876
September 14 2022
1 0

As in the attached

> 

  p1 := proc(a, b, debug::truefalse := false)
             print("p1", debug);
             a + b;
        end proc;
  p2 := proc(a, b, debug::truefalse := false)
             print("p2", debug);
             p1(a, b, debug);
        end proc;
  p2(1, 2, true);

proc (a, b, debug::truefalse := false) print("p1", debug); a+b end proc

  

proc (a, b, debug::truefalse := false) print("p2", debug); p1(a, b, debug) end proc

  

"p2", true

  

"p1", true

  

3

 (1)
> 

 

Download procpass2.mw

Hmmmm...

tomleslie 13876
September 14 2022
1 11
  1. Your "equations" are supplied as a Vector() - why?
  2. The first two terms in each equation appear to be scalars, and these are combined using non-commutative multiplication - the `.` character. Is this for real?
  3. Each equation then adds the (non-commutative) product of two (apparent) scalars to a single element vector. Again, is this for real?

If you just want to solve three "simple" linear equations, then consider the attached.

> 

  restart;
  with(LinearAlgebra):
  eqs:= [ 0.85*(phi__cs*beta__s*f__con * D__pile*w__1) -5/6*H - 5/8*P,
          0.85*(phi__cs*beta__s*f__con * D__pile*w__2) +5/6*H - 5/8*P,
          0.85*(phi__cs*beta__s*f__con * D__pile*w__3) +1/2*H - 3/8*P
        ];
  unks:=Vector([w__1, w__2, w__3]):
  A, b:=GenerateMatrix(convert~(eqs,rational), convert(unks, list) );
  ans:=zip(`=`, unks, simplify(LinearSolve(A, b)))

eqs := [.85*`&phi;__cs`*`&beta;__s`*f__con*D__pile*w__1-5*H*(1/6)-5*P*(1/8), .85*`&phi;__cs`*`&beta;__s`*f__con*D__pile*w__2+5*H*(1/6)-5*P*(1/8), .85*`&phi;__cs`*`&beta;__s`*f__con*D__pile*w__3+(1/2)*H-3*P*(1/8)]

  

A, b := Matrix(3, 3, {(1, 1) = (17/20)*`&phi;__cs`*`&beta;__s`*f__con*D__pile, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = (17/20)*`&phi;__cs`*`&beta;__s`*f__con*D__pile, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = (17/20)*`&phi;__cs`*`&beta;__s`*f__con*D__pile}), Vector(3, {(1) = (5/6)*H+(5/8)*P, (2) = -(5/6)*H+(5/8)*P, (3) = -(1/2)*H+(3/8)*P})

  

Vector[column](%id = 36893488148078529100)

 (1)
> 

 

Download linSol.mw

 

First 15 16 17 18 19 20 21 Last Page 17 of 207