Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 165 days
Contact tomleslie

MaplePrimes Activity


These are replies submitted by tomleslie

You will notice...

tomleslie 13876
April 03 2020
0 0

@Adam Ledger 

that zip archives are allowed. All you have to do is 'zip' your data (ie .csv) file) with your worksheet and then upload here.

If for some reason that is too difficult for you, then chenge the file extension on your csv file to .txt' and upload that instead

This isn't exactly difficult

Hmmm...

tomleslie 13876
April 02 2020
0 0

@Kros 

When I said

Take my worksheet  pdeIss3.mw and run it.

which is the word you don't understand? The worksheet pdeIss3.mw contains the PDE/BCs

    PDE := diff(u(x, t), t) - diff(u(x, t), x$2) - u(x, t) + x - 2*sin(2*x)*cos(x) = 0;
  IBC := D[1](u)(Pi/2, t) = 1, u(0, t) = 0, u(x, 0) = x;

Before you do anything else, you should verify that my response solves your original problem!

In the worksheet you have recently supplied, you have changed the PDE/BCs to something completely different, ie

    PDE := diff(u(x, t), t,t) - diff(u(x, t), x$2) - 10*u(x, t) - 2*sin(2*x)*cos(x) = 0;
  IBC := D[1](u)(Pi/2, t) = 0, u(0, t) = 0, u(x, 0) = 0, D[2](u)(x, 0) = 0;

Why do you believe that just because Maple can analytically solve one PDE/BC combination,  it can analytically solve a completely different PDE/BC combination?

I have actually tried so solve this second  PDE/BC combination, and interrupted Maple after an hour, during which no solution was obtained. I would expect Maple to have more problems (maybe fail??) with this second PDE/BC combination simply because it is second-order in both time and space. The original PDE/BC combination is first order in time and second order in space. Increasing the order of PDEs generally(?) makes them much more difficult to solve.

I think that for this second PDE/BC combination you will be restricted to a numerical-only solution

Show the output...

tomleslie 13876
April 01 2020
0 0

@Kros 

Making statements like "nothing happens" tells me precisely nothing.

Take my worksheet  pdeIss3.mw and run it.

Don't edit it, don't touch it, don't do anything other than run it

And when I say run it I mean leave it running for at least 30mins.If it hasn't worked after 30mins, then kill the computation and upload the worksheet here

Instructions...

tomleslie 13876
April 01 2020
0 0

@Kros 

from here https://www.mapleprimes.com/posts/210295-The-Physics-Updates-How-To-Install

I missed the highlighted instructions in my previous post

  • Open Maple,
  • Click the MapleCloud icon in the upper-right corner to open the MapleCloud toolbar 
  • In the MapleCloud toolbar, open Packages
  • Find the Physics Updates package and click the install button, it is the last one under Actions
  • To check for new versions of Physics Updates, click the MapleCloud icon. If the Updates icon has a red dot, click it to install the new version

Flipping the sign...

tomleslie 13876
April 01 2020
0 0

@c10h20o 

makes a big difference!

Suggest you examine the possibilities in the attached very carefully

> 

  restart;
#
# OP's original expression, with parentheses fixed
#
  expr1:=(7.72-7.72*B)*(-7.717267500*a) = 662204.4444*B^2;
#
# Using 'implicitplot' with NEGATIVE values of 'a'
#
  p1:=plots:-implicitplot(expr1, a = -50000..0, B = 0 .. 1);
#
# Reflect this plot about the y-axis, essentially
# changing the sign of 'a'
#
  plottools:-reflect(p1, [[0,0],[0,1]]);

-7.717267500*(7.72-7.72*B)*a = 662204.4444*B^2

  

 

 
> 

#
# Or one can change the sign of 'a' in the original
# expression by deleting the highlighted minus sign
# in the original expression
#
#  (7.72-7.72*B)*(-7.717267500*a) = 662204.4444*B^2
#
# to get the following, and then use implicitplot()
# with a positive ranges for 'a'
#
  expr2:=(7.72-7.72*B)*(7.717267500*a) = 662204.4444*B^2;
  p2:=plots:-implicitplot(expr2, a = 0..50000, B = 0 .. 1);

7.717267500*(7.72-7.72*B)*a = 662204.4444*B^2

  

 
> 

#
# Or the OP has claimed that he solved for 'a'
# in terms of 'B' and then rotated the result
# anticlockwise by Pi/2, which would be
#
  p3:=plot( rhs(isolate(expr1, a)), B=0..0.85);
  plottools:-rotate(p3, Pi/2);

 

 
> 

 

Download Banda.mw

Someth8ing to try...

tomleslie 13876
April 01 2020
0 0

@Kros 
Because you have you have removed the code at the start of my (perfectly-functional) worksheets, ie this execution group

restart;
version();
Physics:-Version();

and the associated output, I cannot tell exactly which Maple version/subversion you are using. As you can see from the worksheets I have supplied, I am running Maple 2020, build version 1455132, with the Physics updates version 635

However looking at the worksheet you have posted, it is obvious that the command

pdsA:= pdsolve([PDE, IBC]);

is returning nothing, which is why the associated plotting command, ie

p1:= plot3d( rhs(pdsA),
                    x=0..Pi/2,
                    t=0..10,
                    color=red,
                    style=surface
                 );

fails.

Assuming that you are running Maple 2020, the most likely reason that your installation is not returning anything for the analytic solution of your pde, is that you are not running the latest and greatest version of the 'Physics' package. (Amongst other things, this package updates Maple's PDE solvers on a regular basis.)

The "simplest" way to install the latest 'Physics' updates is from the 'cloud' icon in the upper right of your Maple window, immediately below the standard "red cross" to close the application.

Click on this icon, select updates, and install the latest updates ofr the 'Physics' package. You may (or may not) be asked for your MapleSoft account password. Assuming that this installation is successful, restart Maple, and try either of the worksheets I have previously posted, ie

pdeIss2.mw
pdeIss3.mw

again

 

 

As a "table"...

tomleslie 13876
March 31 2020
0 0

@Kros 

Since both the numerical and analytic solution have been determined it is relatively easy to compute these for any given values of the independent variables 'x' and 't' and store the results in any convenient 2-D container.

In the attached, I have extended my previous worksheet to produce matrices of the solutions for both methods, together with a matrix of the 'error' between the two methods

When you ask

In this way can other types of equations be solved ?

I have no idea what this question means!

> 

  restart;
  version();
  Physics:-Version();
  

 User Interface: 1455132
         Kernel: 1455132
        Library: 1455132

  

1455132

  

`The "Physics Updates" version in the MapleCloud is 635 and is the same as the version installed in this computer, created 2020, March 21, 8:19 hours, found in the directory C:\Users\TomLeslie\maple\toolbox\Physics Updates\lib\`

 (1)
> 

#
# Define pde and bcs
#
  PDE := diff(u(x, t), t) - diff(u(x, t), x$2) - u(x, t) + x - 2*sin(2*x)*cos(x) = 0;
  IBC := D[1](u)(Pi/2, t) = 1, u(0, t) = 0, u(x, 0) = x;

diff(u(x, t), t)-(diff(diff(u(x, t), x), x))-u(x, t)+x-2*sin(2*x)*cos(x) = 0

  

(D[1](u))((1/2)*Pi, t) = 1, u(0, t) = 0, u(x, 0) = x

 (2)
> 

#
# Solve analytically (NB this takes a few minutes!)
# and plot the solution
#
  pdsA:= pdsolve([PDE, IBC]);
  p1:= plot3d( rhs(pdsA),
               x=0..Pi/2,
               t=0..10,
               color=red,
               style=surface
             );

`[Length of output exceeds limit of 1000000]`

  

`[Length of output exceeds limit of 1000000]`

  

u(x, t) = (1/8)*(-exp(-8*t)+1)*sin(3*x)+sin(x)*t+x

  

 
> 

#
# Solve numerically and plot the solution
#
  pdsN:= pdsolve(PDE, {IBC}, type = numeric):
  p2:= pdsN:-plot3d( u(x,t),
                     x=0..Pi/2,
                     t=0..10,
                     color=blue,
                     style=surface
                   );

_m8091742464

  

 
> 

#
# Display the analytic and numeric plots
# on the same graph - Pretty similar!
#
  plots:-display([p1,p2]);

 
> 

#
# Plot the absolute error between the analytic
# and numeric solutions
#
  p3:= pdsN:-plot3d( abs( u(x,t)-rhs(pdsA) ),
                     x=0..Pi/2,
                     t=0..10,
                     color=blue,
                     style=surface
                   );

 
> 

#
# Adjust the output display a little, just to make things
# easier to read. Does not affect calculation accuracy, this
# is just for convenience of display
#
  interface(rtablesize=12):
  interface(displayprecision=6):
> 

#
# Define a function which will return values of u(x,t)
# from the numerical pdsolve() procedure for any supplied
# values of the independent variables 'x' and 't'
#
  uNval:= (p,q)->eval( u(x,t),
                        pdsN:-value
                              ( t=q,
                                output=listprocedure
                              )
                     )(p):
#
# Compute the values for u(x,t) from the numerical
# pdsolve() procedure and store these in a matrix.
# Each row of the matrix corresponds to 't' values
# from 0 to 10 in steps of 1. Each column of the
# matrix corresponds to 'x' values from 0 to Pi/2
# in steps of Pi/20
#
  numM:= Matrix
         ( [ seq
             ( [ seq( uNval(i,j),
                      i = 0..evalf(Pi/2)-1.0e-09, evalf(Pi/20)
                    )
               ],
               j=0..10
             )
           ]
         );

Matrix(11, 10, {(1, 1) = 0., (1, 2) = .1570796, (1, 3) = .3141593, (1, 4) = .4712389, (1, 5) = .6283185, (1, 6) = .7853982, (1, 7) = .9424778, (1, 8) = 1.099557, (1, 9) = 1.256637, (1, 10) = 1.413717, (2, 1) = 0., (2, 2) = .3707308, (2, 3) = .7251787, (2, 4) = 1.049865, (2, 5) = 1.336329, (2, 6) = 1.582278, (2, 7) = 1.791454, (2, 8) = 1.972247, (2, 9) = 2.135382, (2, 10) = 2.291190, (3, 1) = 0., (3, 2) = .5275130, (3, 3) = 1.034879, (3, 4) = 1.504852, (3, 5) = 1.925391, (3, 6) = 2.290902, (3, 7) = 2.602185, (3, 8) = 2.865116, (3, 9) = 3.088403, (3, 10) = 3.280895, (4, 1) = 0., (4, 2) = .6843371, (4, 3) = 1.344665, (4, 4) = 1.959971, (4, 5) = 2.514634, (4, 6) = 2.999759, (4, 7) = 3.413198, (4, 8) = 3.758313, (4, 9) = 4.041785, (4, 10) = 4.270985, (5, 1) = 0., (5, 2) = .8412178, (5, 3) = 1.654563, (5, 4) = 2.415254, (5, 5) = 3.104090, (5, 6) = 3.708871, (5, 7) = 4.224503, (5, 8) = 4.651831, (5, 9) = 4.995511, (5, 10) = 5.261431, (6, 1) = 0., (6, 2) = .9981549, (6, 3) = 1.964572, (6, 4) = 2.870701, (6, 5) = 3.693759, (6, 6) = 4.418239, (6, 7) = 5.036100, (6, 8) = 5.545671, (6, 9) = 5.949581, (6, 10) = 6.252235, (7, 1) = 0., (7, 2) = 1.155149, (7, 3) = 2.274694, (7, 4) = 3.326312, (7, 5) = 4.283639, (7, 6) = 5.127862, (7, 7) = 5.847990, (7, 8) = 6.439833, (7, 9) = 6.903995, (7, 10) = 7.243396, (8, 1) = 0., (8, 2) = 1.312199, (8, 3) = 2.584926, (8, 4) = 3.782088, (8, 5) = 4.873733, (8, 6) = 5.837742, (8, 7) = 6.660173, (8, 8) = 7.334317, (8, 9) = 7.858753, (8, 10) = 8.234914, (9, 1) = 0., (9, 2) = 1.469306, (9, 3) = 2.895271, (9, 4) = 4.238027, (9, 5) = 5.464039, (9, 6) = 6.547877, (9, 7) = 7.472648, (9, 8) = 8.229124, (9, 9) = 8.813854, (9, 10) = 9.226789, (10, 1) = 0., (10, 2) = 1.626469, (10, 3) = 3.205728, (10, 4) = 4.694131, (10, 5) = 6.054558, (10, 6) = 7.258267, (10, 7) = 8.285417, (10, 8) = 9.124254, (10, 9) = 9.769300, (10, 10) = 10.219022, (11, 1) = 0., (11, 2) = 1.783689, (11, 3) = 3.516296, (11, 4) = 5.150399, (11, 5) = 6.645289, (11, 6) = 7.968914, (11, 7) = 9.098478, (11, 8) = 10.019705, (11, 9) = 10.725091, (11, 10) = 11.211612})

 (3)
> 

#
# Define a function which will return values of u(x,t)
# from the analytic pdsolve() procedure for any supplied
# values of the independent variables 'x' and 't'
#
  uAval:=unapply(rhs(pdsA), [x,t]):
#
# Compute the values for u(x,t) from the analytical
# pdsolve() procedure and store these in a matrix.
# Each row of the matrix corresponds to 't' values
# from 0 to 10 in steps of 1. Each column of the
# matrix corresponds to 'x' values from 0 to Pi/2
# in steps of Pi/20
#
  analM:= Matrix
          ( evalf~( [ seq
                      ( [ seq( uAval(i,j),
                               i = 0..evalf(Pi/2)-1.0e-09, evalf(Pi/20)
                             )
                        ],
                        j=0..10
                      )
                    ]
                  )
               );

Matrix(11, 10, {(1, 1) = 0., (1, 2) = .157080, (1, 3) = .314159, (1, 4) = .471239, (1, 5) = .628319, (1, 6) = .785398, (1, 7) = .942478, (1, 8) = 1.099557, (1, 9) = 1.256637, (1, 10) = 1.413717, (2, 1) = 0., (2, 2) = .370244, (2, 3) = .724269, (2, 4) = 1.048649, (2, 5) = 1.334946, (2, 6) = 1.580864, (2, 7) = 1.790109, (2, 8) = 1.971016, (2, 9) = 2.134245, (2, 10) = 2.290067, (3, 1) = 0., (3, 2) = .526697, (3, 3) = 1.033320, (3, 4) = 1.502681, (3, 5) = 1.922771, (3, 6) = 2.288000, (3, 7) = 2.599139, (3, 8) = 2.862016, (3, 9) = 3.085277, (3, 10) = 3.277718, (4, 1) = 0., (4, 2) = .683132, (4, 3) = 1.342337, (4, 4) = 1.956671, (4, 5) = 2.510556, (4, 6) = 2.995107, (4, 7) = 3.408156, (4, 8) = 3.753023, (4, 9) = 4.036333, (4, 10) = 4.265406, (5, 1) = 0., (5, 2) = .839566, (5, 3) = 1.651354, (5, 4) = 2.410662, (5, 5) = 3.098342, (5, 6) = 3.702214, (5, 7) = 4.217173, (5, 8) = 4.644029, (5, 9) = 4.987390, (5, 10) = 5.253094, (6, 1) = 0., (6, 2) = .996001, (6, 3) = 1.960371, (6, 4) = 2.864652, (6, 5) = 3.686127, (6, 6) = 4.409320, (6, 7) = 5.026190, (6, 8) = 5.535036, (6, 9) = 5.938446, (6, 10) = 6.240783, (7, 1) = 0., (7, 2) = 1.152435, (7, 3) = 2.269388, (7, 4) = 3.318643, (7, 5) = 4.273912, (7, 6) = 5.116427, (7, 7) = 5.835207, (7, 8) = 6.426042, (7, 9) = 6.889503, (7, 10) = 7.228471, (8, 1) = 0., (8, 2) = 1.308870, (8, 3) = 2.578405, (8, 4) = 3.772633, (8, 5) = 4.861697, (8, 6) = 5.823534, (8, 7) = 6.644224, (8, 8) = 7.317049, (8, 9) = 7.840560, (8, 10) = 8.216159, (9, 1) = 0., (9, 2) = 1.465304, (9, 3) = 2.887422, (9, 4) = 4.226624, (9, 5) = 5.449483, (9, 6) = 6.530641, (9, 7) = 7.453241, (9, 8) = 8.208055, (9, 9) = 8.791616, (9, 10) = 9.203848, (10, 1) = 0., (10, 2) = 1.621739, (10, 3) = 3.196439, (10, 4) = 4.680614, (10, 5) = 6.037268, (10, 6) = 7.237748, (10, 7) = 8.262258, (10, 8) = 9.099062, (10, 9) = 9.742673, (10, 10) = 10.191536, (11, 1) = 0., (11, 2) = 1.778173, (11, 3) = 3.505456, (11, 4) = 5.134605, (11, 5) = 6.625053, (11, 6) = 7.944854, (11, 7) = 9.071275, (11, 8) = 9.990068, (11, 9) = 10.693729, (11, 10) = 11.179224})

 (4)
> 

#
# Compute the difference between the matrix of numerical
# solutions and the matrix of analytical solutions. This
# is one measure of the "error" between the two solution
# methods
#
  analM-numM;

Matrix(11, 10, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (1, 6) = 0.1110223000e-15, (1, 7) = 0., (1, 8) = 0., (1, 9) = 0., (1, 10) = 0., (2, 1) = 0., (2, 2) = -0.487e-3, (2, 3) = -0.909e-3, (2, 4) = -0.1216433e-2, (2, 5) = -0.1383076e-2, (2, 6) = -0.1414418e-2, (2, 7) = -0.1345005e-2, (2, 8) = -0.1230946e-2, (2, 9) = -0.1137305e-2, (2, 10) = -0.1123447e-2, (3, 1) = 0., (3, 2) = -0.816e-3, (3, 3) = -0.1558834e-2, (3, 4) = -0.2171238e-2, (3, 5) = -0.2619747e-2, (3, 6) = -0.2902076e-2, (3, 7) = -0.3045684e-2, (3, 8) = -0.3100299e-2, (3, 9) = -0.3125648e-2, (3, 10) = -0.3177129e-2, (4, 1) = 0., (4, 2) = -0.1205289e-2, (4, 3) = -0.2327930e-2, (4, 4) = -0.3299680e-2, (4, 5) = -0.4078033e-2, (4, 6) = -0.4652070e-2, (4, 7) = -0.5041620e-2, (4, 8) = -0.5289963e-2, (4, 9) = -0.5451751e-2, (4, 10) = -0.5578783e-2, (5, 1) = 0., (5, 2) = -0.1651448e-2, (5, 3) = -0.3208671e-2, (5, 4) = -0.4592144e-2, (5, 5) = -0.5748675e-2, (5, 6) = -0.6657526e-2, (5, 7) = -0.7329828e-2, (5, 8) = -0.7801516e-2, (5, 9) = -0.8121429e-2, (5, 10) = -0.8337248e-2, (6, 1) = 0., (6, 2) = -0.2154144e-2, (6, 3) = -0.4201092e-2, (6, 4) = -0.6048681e-2, (6, 5) = -0.7631744e-2, (6, 6) = -0.8918530e-2, (6, 7) = -0.9910414e-2, (6, 8) = -0.1063507e-1, (6, 9) = -0.1113481e-1, (6, 10) = -0.1145264e-1, (7, 1) = 0., (7, 2) = -0.2713397e-2, (7, 3) = -0.5305235e-2, (7, 4) = -0.7669353e-2, (7, 5) = -0.9727318e-2, (7, 6) = -0.1143518e-1, (7, 7) = -0.1278348e-1, (7, 8) = -0.1379075e-1, (7, 9) = -0.1449202e-1, (7, 10) = -0.1492511e-1, (8, 1) = 0., (8, 2) = -0.3329227e-2, (8, 3) = -0.6521137e-2, (8, 4) = -0.9454217e-2, (8, 5) = -0.1203547e-1, (8, 6) = -0.1420756e-1, (8, 7) = -0.1594914e-1, (8, 8) = -0.1726867e-1, (8, 9) = -0.1819318e-1, (8, 10) = -0.1875477e-1, (9, 1) = 0., (9, 2) = -0.4001655e-2, (9, 3) = -0.7848843e-2, (9, 4) = -0.1140333e-1, (9, 5) = -0.1455628e-1, (9, 6) = -0.1723576e-1, (9, 7) = -0.1940749e-1, (9, 8) = -0.2106894e-1, (9, 9) = -0.2223841e-1, (9, 10) = -0.2294175e-1, (10, 1) = 0., (10, 2) = -0.4730702e-2, (10, 3) = -0.9288389e-2, (10, 4) = -0.1351676e-1, (10, 5) = -0.1728983e-1, (10, 6) = -0.2051988e-1, (10, 7) = -0.2315864e-1, (10, 8) = -0.2519168e-1, (10, 9) = -0.2662784e-1, (10, 10) = -0.2748618e-1, (11, 1) = 0., (11, 2) = -0.5516386e-2, (11, 3) = -0.1083982e-1, (11, 4) = -0.1579456e-1, (11, 5) = -0.2023619e-1, (11, 6) = -0.2406002e-1, (11, 7) = -0.2720270e-1, (11, 8) = -0.2963701e-1, (11, 9) = -0.3136160e-1, (11, 10) = -0.3238819e-1})

 (5)
> 

 


 

Download pdeIss3.mw

 

Yes...

tomleslie 13876
March 31 2020
0 0

@Kros 

You can put the analytic and numeric solutions on a single graph, and even plot the discrepancy between the two

See the attached

> 

  restart;
  version();
  Physics:-Version();
  

 User Interface: 1455132
         Kernel: 1455132
        Library: 1455132

  

1455132

  

`The "Physics Updates" version in the MapleCloud is 635 and is the same as the version installed in this computer, created 2020, March 21, 8:19 hours, found in the directory C:\Users\TomLeslie\maple\toolbox\Physics Updates\lib\`

 (1)
> 

#
# Define pde and bcs
#
  PDE := diff(u(x, t), t) - diff(u(x, t), x$2) - u(x, t) + x - 2*sin(2*x)*cos(x) = 0;
  IBC := D[1](u)(Pi/2, t) = 1, u(0, t) = 0, u(x, 0) = x;

diff(u(x, t), t)-(diff(diff(u(x, t), x), x))-u(x, t)+x-2*sin(2*x)*cos(x) = 0

  

(D[1](u))((1/2)*Pi, t) = 1, u(0, t) = 0, u(x, 0) = x

 (2)
> 

#
# Solve analytically (NB this takes a few minutes!)
# and plot the solution
#
  pdsA:= pdsolve([PDE, IBC]);
  p1:= plot3d( rhs(pdsA),
               x=0..Pi/2,
               t=0..10,
               color=red,
               style=surface
             );

u(x, t) = (1/8)*(-exp(-8*t)+1)*sin(3*x)+sin(x)*t+x

  

 
> 

#
# Solve numerically and plot the solution
#
  pdsN:= pdsolve(PDE, {IBC}, type = numeric);
  p2:= pdsN:-plot3d( u(x,t),
                     x=0..Pi/2,
                     t=0..10,
                     color=blue,
                     style=surface
                   );

_m1918883328

  

 
> 

#
# Display the analytic and numeric plots
# on the same graph - Pretty similar!
#
  plots:-display([p1,p2]);

 
> 

#
# Plot the absolute error between the analytic
# and numeric solutions
#
  p3:= pdsN:-plot3d( abs( u(x,t)-rhs(pdsA) ),
                     x=0..Pi/2,
                     t=0..10,
                     color=blue,
                     style=surface
                   );

 
> 

 


Download pdeIss2.mw

You still have...

tomleslie 13876
March 30 2020
0 0

@Zeineb 

12 equations in 16 unknowns?

What are u[4, 0], u[4, 1], u[4, 2], u[4, 3]? These quantities are not in your vector of unknowns, or in your lit of parameters.

See the attached

> 

restart; with(LinearAlgebra)
> 

for i to 3 do eq[i] := -alpha*(u[i+1, 1]+u[i-1, 1])-beta*(u[i, 2]+u[i, 0])+4*u[i, 1] end do

-alpha*(u[2, 1]+u[0, 1])-beta*(u[1, 2]+u[1, 0])+4*u[1, 1]

  

-alpha*(u[3, 1]+u[1, 1])-beta*(u[2, 2]+u[2, 0])+4*u[2, 1]

  

-alpha*(u[4, 1]+u[2, 1])-beta*(u[3, 2]+u[3, 0])+4*u[3, 1]

 (1)
> 

for i to 3 do eq[i+3] := -alpha*(u[i+1, 2]+u[i-1, 2])-beta*(u[i, 3]+u[i, 1])+4*u[i, 2] end do

-alpha*(u[2, 2]+u[0, 2])-beta*(u[1, 3]+u[1, 1])+4*u[1, 2]

  

-alpha*(u[3, 2]+u[1, 2])-beta*(u[2, 3]+u[2, 1])+4*u[2, 2]

  

-alpha*(u[4, 2]+u[2, 2])-beta*(u[3, 3]+u[3, 1])+4*u[3, 2]

 (2)
> 

for i to 3 do eq[i+6] := -alpha*(u[i-1, 0]+u[i+1, 0])-2*beta*u[i, 1]+4*u[i, 0] end do

-alpha*(u[0, 0]+u[2, 0])-2*beta*u[1, 1]+4*u[1, 0]

  

-alpha*(u[1, 0]+u[3, 0])-2*beta*u[2, 1]+4*u[2, 0]

  

-alpha*(u[2, 0]+u[4, 0])-2*beta*u[3, 1]+4*u[3, 0]

 (3)
> 

for i to 3 do eq[i+9] := -alpha*(u[i-1, 3]+u[i+1, 3])-2*beta*u[i, 2]+4*u[i, 3] end do

-alpha*(u[0, 3]+u[2, 3])-2*beta*u[1, 2]+4*u[1, 3]

  

-alpha*(u[1, 3]+u[3, 3])-2*beta*u[2, 2]+4*u[2, 3]

  

-alpha*(u[2, 3]+u[4, 3])-2*beta*u[3, 2]+4*u[3, 3]

 (4)
> 

eqsA := [seq(eq[j], j = 1 .. 12)]; eqsB := eval(eqsA, [alpha = 5, beta = 10, u[0, 0] = 100, u[0, 1] = 100, u[0, 2] = 100, u[0, 3] = 100]); numelems(eqsB); `union`(`~`[indets](eqsB, name)[]); numelems(%)

[-alpha*(u[2, 1]+u[0, 1])-beta*(u[1, 2]+u[1, 0])+4*u[1, 1], -alpha*(u[3, 1]+u[1, 1])-beta*(u[2, 2]+u[2, 0])+4*u[2, 1], -alpha*(u[4, 1]+u[2, 1])-beta*(u[3, 2]+u[3, 0])+4*u[3, 1], -alpha*(u[2, 2]+u[0, 2])-beta*(u[1, 3]+u[1, 1])+4*u[1, 2], -alpha*(u[3, 2]+u[1, 2])-beta*(u[2, 3]+u[2, 1])+4*u[2, 2], -alpha*(u[4, 2]+u[2, 2])-beta*(u[3, 3]+u[3, 1])+4*u[3, 2], -alpha*(u[0, 0]+u[2, 0])-2*beta*u[1, 1]+4*u[1, 0], -alpha*(u[1, 0]+u[3, 0])-2*beta*u[2, 1]+4*u[2, 0], -alpha*(u[2, 0]+u[4, 0])-2*beta*u[3, 1]+4*u[3, 0], -alpha*(u[0, 3]+u[2, 3])-2*beta*u[1, 2]+4*u[1, 3], -alpha*(u[1, 3]+u[3, 3])-2*beta*u[2, 2]+4*u[2, 3], -alpha*(u[2, 3]+u[4, 3])-2*beta*u[3, 2]+4*u[3, 3]]

  

[-5*u[2, 1]-500-10*u[1, 2]-10*u[1, 0]+4*u[1, 1], -5*u[3, 1]-5*u[1, 1]-10*u[2, 2]-10*u[2, 0]+4*u[2, 1], -5*u[4, 1]-5*u[2, 1]-10*u[3, 2]-10*u[3, 0]+4*u[3, 1], -5*u[2, 2]-500-10*u[1, 3]-10*u[1, 1]+4*u[1, 2], -5*u[3, 2]-5*u[1, 2]-10*u[2, 3]-10*u[2, 1]+4*u[2, 2], -5*u[4, 2]-5*u[2, 2]-10*u[3, 3]-10*u[3, 1]+4*u[3, 2], -500-5*u[2, 0]-20*u[1, 1]+4*u[1, 0], -5*u[1, 0]-5*u[3, 0]-20*u[2, 1]+4*u[2, 0], -5*u[2, 0]-5*u[4, 0]-20*u[3, 1]+4*u[3, 0], -500-5*u[2, 3]-20*u[1, 2]+4*u[1, 3], -5*u[1, 3]-5*u[3, 3]-20*u[2, 2]+4*u[2, 3], -5*u[2, 3]-5*u[4, 3]-20*u[3, 2]+4*u[3, 3]]

  

12

  

{u[1, 0], u[1, 1], u[1, 2], u[1, 3], u[2, 0], u[2, 1], u[2, 2], u[2, 3], u[3, 0], u[3, 1], u[3, 2], u[3, 3], u[4, 0], u[4, 1], u[4, 2], u[4, 3]}

  

16

 (5)
> 

"# #` OP seems uninterested (unaware?) of the unknowns `u[4,0],u[4,1],u[4,2],u[4,3] #` since (s)he generates the unknowns as` # #` vect_unkown:=`[u[1,0],u[2,0],u[3,0], u[1,1],u[2,1],u[3,1],u[1,2],u[2,2],u[3,2],u[1,3],u[2,3],u[3,3]]; # #` which can be written much more compactly as` # vect_unknown:=[seq( seq(u[i,j], i=1..3), j=0..3)]; # #` \Produce the coefficient matrix and the "right-hand-side"` #` vector. Note that the latter contains the unknowns` #` u`[4,0],u[4,1],u[4,2],u[4,3] of which the OP is unaware"?" # interface(rtablesize=12): A,b := GenerateMatrix(eqsB, vect_unknown)"

vect_unknown := [u[1, 0], u[2, 0], u[3, 0], u[1, 1], u[2, 1], u[3, 1], u[1, 2], u[2, 2], u[3, 2], u[1, 3], u[2, 3], u[3, 3]]

  

Matrix(%id = 18446744074426915286), Vector[column](%id = 18446744074426915166)

 (6)
> 

``

Download sysToMat.mw

Still not good enough...

tomleslie 13876
March 29 2020
0 0

@Zeineb 

If you use

indets(eqs, 'name')

this will return

{alpha, beta, u[0, 0], u[0, 1], u[0, 2], u[0, 3], u[1, 0], u[1, 1], u[1, 2], u[1, 3], u[2, 0], u[2, 1], u[2, 2], u[2, 3], u[3, 0], u[3, 1], u[3, 2], u[3, 3], u[4, 0], u[4, 1], u[4, 2], u[4, 3]}

However you define 'vect_unkown' as

[u[1, 0], u[2, 0], u[3, 0], u[1, 1], u[2, 1], u[3, 1], u[1, 2], u[2, 2], u[3, 2], u[1, 3], u[2, 3], u[3, 3]]

So what exactly are represented by the quantities alpha, beta, u[0, 0], u[0, 1], u[0, 2], u[0, 3]???

and interestingly...

tomleslie 13876
March 27 2020
0 0

@acer 

your results differ from my original post.

I screwed up when replacing '[]' parentheses.

Corrected(?) version is attached below - which at least now gives the same answers as yours

> 

  restart:
  with(RootFinding):
> 

  r[eD] := 4.5:
  rootfunction:= beta -> BesselY(1, beta)*BesselJ(1, beta*r[eD]) - BesselJ(1, beta)*BesselY(1, beta*r[eD]):
  a := Vector(40, fill = 0):
  StartPoint:=0.001:
  for i from 1 by 1 to 40 do
      a[i]:= NextZero(rootfunction, StartPoint):
      StartPoint:=a[i]:
  od:
  a;

_rtable[18446744074333913206]

 (1)

``

> 

  p__D:= (r__D, t__D) ->(2*((r__D^2)/4+t__D))/(r[eD]^2-1)-(r[eD]^2*ln(r__D))/(r[eD]^2-1)-(3*r[eD]^4-4*r[eD]^4*ln(r[eD])-2*r[eD]^2-1)/(4*(r[eD]^2-1)^2)+Pi*add(BesselJ(1,a[i]*r[eD])^2*(BesselJ(1,a[i])*BesselY(0,a[i]*r__D)-BesselY(1,a[i])*BesselJ(0,a[i]*r__D))/(a[i]*(BesselJ(1,a[i]*r[eD])^2-BesselJ(1,a[i])^2))*exp(-a[i]^2*t__D),  i = 1..40  );

proc (r__D, t__D) local i; options operator, arrow; 2*((1/4)*r__D^2+t__D)/(r[eD]^2-1)-r[eD]^2*ln(r__D)/(r[eD]^2-1)-(1/4)*(3*r[eD]^4-4*r[eD]^4*ln(r[eD])-2*r[eD]^2-1)/(r[eD]^2-1)^2+Pi*add(BesselJ(1, a[i]*r[eD])^2*(BesselJ(1, a[i])*BesselY(0, a[i]*r__D)-BesselY(1, a[i])*BesselJ(0, a[i]*r__D))*exp(-a[i]^2*t__D)/(a[i]*(BesselJ(1, a[i]*r[eD])^2-BesselJ(1, a[i])^2)), i = 1 .. 40) end proc

 (2)
> 

#
# Check a few values
#
  p__D(1,1);
  p__D(0.5,0.5);

.8021455142

  

1.384740257

 (3)

NULL

``

Download diffFunc2.mw

 

Yes...

tomleslie 13876
March 25 2020
0 0

@lcz 

It seems as if the GraphTheory package has had a significant upgrade in Maple 2020 - although a lot of it is (just?) layout+visualization

See the attached for a graph with a loop.


 

> 

  restart;
  with(GraphTheory):
  G:=Graph({[1,2],[1,5],[2,5],[2,3],[5,4],[3,4],[4,6]}):
  HasSelfLoop(G);
  DrawGraph(G,style=spring);
  G:=Graph({[1,1],[1,2],[1,5],[2,5],[2,3],[5,4],[3,4],[4,6]}):
  HasSelfLoop(G);
  DrawGraph(G,style=spring);

false

  

 

true

  

 
> 

 

Download gr2020.mw

 

Basically yes...

tomleslie 13876
March 25 2020
0 0

@Marliberty28 

It would be much simpler to provide the optimal solution if you used the big green up-arrow in the Mapleprimes toolbar to upload  the worksheet with your specific problem

The *best* method for obtaining the roots of any expression, may depend ((a little) on the precise nature of the expression

There is nothing particularly difficult about obtaining the roots of an ODE solution, although the precise method depends on whether you are solving the ODE analytically or numerically.

See the 'toy' examples in the attached

> 

  restart;
#
# Define a "toy" ODE with boundary condition
#
  ode:=diff(f(x),x)=cos(10*x);
  bc:=f(0)=0:
#
# Solve this analytically
#
  sol1:=unapply(rhs(dsolve([ode, bc])),x);
#
# Obtain the roots of the solution
#
  Student[Calculus1]:-Roots( sol1(x), x=-1 .. 1 );

diff(f(x), x) = cos(10*x)

  

proc (x) options operator, arrow; (1/10)*sin(10*x) end proc

  

[-(3/10)*Pi, -(1/5)*Pi, -(1/10)*Pi, 0, (1/10)*Pi, (1/5)*Pi, (3/10)*Pi]

 (1)
> 

#
# Solve the same ODE numerically
#
  sol2:=eval( f(x), dsolve([ode, bc], numeric, output=listprocedure));
#
# Obtain the roots of the solution
#
  Student[Calculus1]:-Roots( sol2(x), x=-1 .. 1 );

proc (x) local _res, _dat, _solnproc, _xout, _ndsol, _pars, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x) else _xout := evalf(x) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 26, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..63, {(1) = 1, (2) = 1, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 1, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = 0.5047658755841546e-2, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..1, {(1) = .0}, datatype = float[8], order = C_order)), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..1, {(1) = .1}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8]), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8])]), ( 8 ) = ([Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 1.0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..1, {(1, 1) = .0, (2, 0) = .0, (2, 1) = .0, (3, 0) = .0, (3, 1) = .0, (4, 0) = .0, (4, 1) = .0, (5, 0) = .0, (5, 1) = .0, (6, 0) = .0, (6, 1) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = f(x)]`; YP[1] := cos(10*X); 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = f(x)]`; YP[1] := cos(10*X); 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 26 ) = (Array(1..0, {})), ( 25 ) = (Array(1..0, {})), ( 24 ) = (0) ] )) ] ); _y0 := Array(0..1, {(1) = 0.}); _vmap := array( 1 .. 1, [( 1 ) = (1) ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); _i := false; if _par <> [] then _i := `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then _i := `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) or _i end if; if _i then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch: end try else try procname(procname("right")) catch: end try end if end if end if; return elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch: end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch: end try end if end if end if; return elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error end try; if _dat[17] <> _dtbl[1][17] then _dtbl[1][17] := _dat[17]; _dtbl[1][10] := _dat[10] end if; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; if type(_EnvDSNumericSaveDigits, 'posint') then _dat[4][26] := _EnvDSNumericSaveDigits else _dat[4][26] := Digits end if; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(1..2, {(1) = 18446744074403566270, (2) = 18446744074403566446}), (3) = [x, f(x)], (4) = []}); _solnproc := _dat[1]; _pars := map(rhs, _dat[4]); if not type(_xout, 'numeric') then if member(x, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x, ["last", 'last', "initial", 'initial', NULL]) then _res := _solnproc(convert(x, 'string')); if type(_res, 'list') then return _res[2] else return NULL end if elif member(x, ["parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [_res[2], seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x), 'string') = rhs(x); if lhs(_xout) = "initial" then if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else _res := _solnproc("initial" = ["single", 2, rhs(_xout)]) end if elif not type(rhs(_xout), 'list') then error "initial and/or parameter values must be specified in a list" elif lhs(_xout) = "initial_and_parameters" and nops(rhs(_xout)) = nops(_pars)+1 then _res := _solnproc(lhs(_xout) = ["single", 2, op(rhs(_xout))]) else _res := _solnproc(_xout) end if; if lhs(_xout) = "initial" then return _res[2] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [_res[2], seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x), 'string') = rhs(x)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _dat[3] end if; if procname <> unknown then return ('procname')(x) else _ndsol := 1; _ndsol := `tools/gensym`("f(x)"); eval(FromInert(_Inert_FUNCTION(_Inert_NAME("assign"), _Inert_EXPSEQ(ToInert(_ndsol), _Inert_VERBATIM(pointto(_dat[2][2])))))); return FromInert(_Inert_FUNCTION(ToInert(_ndsol), _Inert_EXPSEQ(ToInert(x)))) end if end if; try _res := _solnproc(_xout); _res[2] catch: error end try end proc

  

Warning, some roots are returned as numeric approximations

  

[-.9424777968, -.6283184968, -.3141593051, 0., .3141593051, .6283184968, .9424777968]

 (2)
> 

 


 

Download odeRoot.mw

the wrong way round!...

tomleslie 13876
March 24 2020
0 0

You should not

  1. define a function
  2. plot the function
  3. extract 'zeros' from the plot

What you should do is

  1. define a function
  2. extract 'zeros' from the function

If you uploaded the worksheet containing your function definition using the big green up-arrow in the Mapleprimes toolbar, then explicit help could be provided. SInce you can't/won't/don't uopload anything specific, the best I can suggest is that the RootFinding:-NextZero() command is probably a good place to start

So...

tomleslie 13876
March 23 2020
0 0

@nm 

Why is the default format for new worksheets set as Document Mode in you Options->Interface menu?

This makes no sense!

First 53 54 55 56 57 58 59 Last Page 55 of 207