acer

32303 Reputation

29 Badges

19 years, 307 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

It looks like you have mistakenly split the do-loop code across Document Blocks (or Execution Groups, though there's a slight hint in your posting that you might be pasting into Document Blocks).

Don't press Enter/Return, while typing/pasting in the Block/Group with that do-loop code. That would execute what's entered so far, and split it.

You can use Shift-Enter to get to a new line, when editing in a Block/Group.

If you look at the example in the subsection "Getting the Number of Elements in a List" of Section 4.3 or the Programming Guide (which seems to be your example), then you can see that all the red code of the while-do loop is within the same Execution Group.

That is, the beginning while, and the end do, and all the loop's contents between those, must be within the same Execution Group (or Document Block). Otherwise you'd get those kinds of error messages.

Here's what it might look like, if accidentally split. Notice the left border of the Maple GUI can show Markers that denote the Block/Group boundaries.

Download doc_block_oops.mw

I dug around and found this in the vault (edited slightly to get rid of the infinity, which solve only sometimes removes).

This kind of thing has been asked before, both the basic query about conversion of a real range to inequality relations, as well as simplification of compound examples.

The basic examples are a common flavor, for such there is a dedicated command, as has been mentioned: convert(...,relation).

A related idea is that one can do some simplification for some more involved examples. Also, on output, multiple conditions and multiple solutions can be displayed using Or/And, or %or if one wants, or (last link in end note) union/intersect.

I put the RealRange inputs into 2D form, for fun.

restart;

 

G := proc(R)
   # extra stuff to beautify multiple output
   (u->ifelse(nops(u)>1,%or(op(u)),op(u)))(
      map(u->ifelse(u::{set,list} and nops(u)>1,
                    `%and`(op(u)),op(u)),
          [solve(subsindets(convert(R,relation),
                            And({`<`,`<=`},satisfies(u->has(u,infinity))),
                            ()->NULL),{x})]));
end proc:


Now, shown as 2D Input,

G(x::(RealRange(1, Open(4))))

%and(1 <= x, x < 4)

G(x::(RealRange(Open(2), Open(infinity))))

2 < x

G(x::(RealRange(2, 3)))

%and(2 <= x, x <= 3)

G(x::(RealRange(Open(-infinity), -3)))

x <= -3


Alternatively,

G(`in`(x, RealRange(1, Open(4))))

%and(1 <= x, x < 4)

G(`in`(x, RealRange(Open(2), Open(infinity))))

2 < x

G(`in`(x, RealRange(2, 3)))

%and(2 <= x, x <= 3)

G(`in`(x, RealRange(Open(-infinity), -3)))

x <= -3


Some fun,

G(`or`(x::(RealRange(2, Open(infinity))), x::(RealRange(-3, 3))))

-3 <= x

G(`and`(x::(RealRange(2, Open(infinity))), x::(RealRange(-3, 3))))

%and(2 <= x, x <= 3)

G(`or`(x::(RealRange(Open(-infinity), -2)), x::(RealRange(-1, 3))))

%or(x <= -2, %and(-1 <= x, x <= 3))

G(`or`(x::(RealRange(Open(-infinity), -2)), x::(RealRange(-3, 3))))

x <= 3

G(`and`(x::(RealRange(-infinity, -2)), x::(RealRange(-3, 3))))

%and(-3 <= x, x <= -2)

G(`or`(`in`(x, RealRange(-infinity, 6)), `in`(x, RealRange(-4, 7))))

x <= 7

G(`and`(`in`(x, RealRange(-infinity, -2)), `in`(x, RealRange(-3, 3))))

%and(-3 <= x, x <= -2)

G(Or(`in`(x, RealRange(-infinity, 6)), `in`(x, RealRange(-4, 7))))

x <= 7

G(And(`in`(x, RealRange(-infinity, -2)), `in`(x, RealRange(-3, 3))))

%and(-3 <= x, x <= -2)

Download some_2D_range_ineq.mw

[edit] It's it's of interest, here are a few old postings on related (or more involved) queries: 1, 2, 3

Given the input ode, why not go straight to,

   has(convert(ode,D),y(x));

without any replacement by Z^N terms.

The simplest (and most effective and efficient way, IMO) is to use a so-called operator form calling sequence of the evalf(Int(...)) or int(...,numeric) or int(...,0..6.0) commands, where the integrand is supplied as an operator/procedure f rather than a function call like f(x).

There are several ways to get it to work instead with something that looks like an unevaluated function call (of x, like f(x), say) as the integrand. But they are often unnecessarily more complicated than the above. I'll show a few such ways.

restart

with(Interpolation)

X := [2, 4, 8, 12, 14, 20]

Y := [1, 3, 11, 8, 19, 23]

Z__lin := LinearInterpolation(X, Y)

Z__spline := SplineInterpolation(X, Y)


These first ones are all pretty straightforward, easy, and terse, IMO.

These first four examples all use calling sequences of the command
with so-called operator-form for the integrand.

I'll do it for both your interpolations.

evalf(Int(Z__lin, 0 .. 6))

HFloat(13.999999999999517)

evalf(Int(Z__spline, 0 .. 6))

HFloat(14.878464818761161)

int(Z__lin, 0 .. 6, numeric)

HFloat(13.999999999999517)

int(Z__spline, 0 .. 6, numeric)

HFloat(14.878464818761161)

int(Z__lin, 0 .. 6.0)

HFloat(13.999999999999517)

int(Z__spline, 0 .. 6.0)

HFloat(14.878464818761161)


Now, for brevity, I'll just deal with just one of your interpolations,
while showing some so-called expression-form (function call) forms
for the integrand..

The next needs exactly two pairs of single quotes-quotes.

int((''Z__lin'')(x), x = 0 .. 6.0)

HFloat(13.999999999999517)


The next needs exactly three pairs of single right-quotes.

evalf(int(('''Z__lin''')(x), x = 0 .. 6))

HFloat(13.999999999999517)


The next MakeFunction (or unapply) require a special incantation.

fZ__lin := MakeFunction('Z__lin(x)', x, numeric)

int(fZ__lin(x), x = 0 .. 6.0)

HFloat(13.999999999999517)

evalf(int(fZ__lin(x), x = 0 .. 6))

HFloat(13.999999999999517)

NULL

Download Using_interpolation_ac.mw

You could add the dirfield option, eg.
   dirfield=[11,11]

And you can blend this with other options, eg.
   'arrows'='comet',dirfield=[11,11],arrowsize=1.2

You can save all the assigned names (including just user-defined names), but there's no easy way to save&restore all internal state including memoized results in remember tables and caches.

That's not to say that it's generally impossible (for, say, an expert), especially for Library routines. However there's no stock command for it.

Saving/restoring all the remembered results in all kernel built-ins (eg. evalf, series, etc) might be much harder and perhaps not possible.

Since there are several flavors to this, it would be better if you supplied a concrete, representative example.

Your eval calls attempt to replace names like f__1, f__2, etc, but your expression M actually contains f[1], f[2], etc.

plot_a_ac.mw

Here are a few ways, which you might adjust with options, etc.

My "trick" is to create a system for each param1 value, instead of trying to use it as a DS TransferFunction `parameter` in that routine's technical sense.

[edit: I added plots:-animate]

simple_bode_ac.mw

Are any of these the kind of thing you might be after?

If not, could you provide an example, with input, output, and desired "order" answer? (A worksheet might be best.)

restart;

A := asympt(exp(x^2)*(1-erf(x)),x);

1/(Pi^(1/2)*x)-(1/2)/(Pi^(1/2)*x^3)+(3/4)/(Pi^(1/2)*x^5)+O(1/x^7)

order(series(A,x));

-7

degree(op(indets(A,specfunc(O))[1]),x);

-7

B := MultiSeries:-series(eval(exp(x^2)*(1-erf(x)),x=1/x),
                         x,5):
eval(B,x=1/x);

1/(Pi^(1/2)*x)-(1/2)/(Pi^(1/2)*x^3)+(3/4)/(Pi^(1/2)*x^5)+O(1/x^7)

order(B);

7

Download asympt_order_q.mw

Using the global (top-level) variable Metric is quite poor.

For one thing, since a Matrix is a mutable data structure then some entry in the toplevel :-Metric might change. And then any already computed invMetric would be wrong. So every call to your module's computational routines would sensibly need to compare entries between :-Metric and some copy (so as to know whether to recompute invMetric. Yikes.

A more usual approach (without "resorting" to Objects) is to have an exported utility which sets both the stored local Metric and its inverse. Then you would simply invoke that where you would previously have redefined :-Metric in your earlier scheme of things.

restart

``

TM := module () local Metric, invMetric; export foo, bar, SetMetric;  Metric := Matrix(3, shape = symmetric, [[1, 0, 0], [0, 1, 0], [0, 0, 1]]); invMetric := LinearAlgebra:-MatrixInverse(Metric); SetMetric := proc (M::(Matrix(3))) if not EqualEntries(Metric, M) then print("updating Metric"); Metric := copy(M); invMetric := LinearAlgebra:-MatrixInverse(Metric) end if; return NULL end proc; foo := proc () print('Metric' = Metric) end proc; bar := proc () print('invMetric' = invMetric) end proc end module

TM:-foo()

Metric = Matrix(%id = 36893627907715573268)

TM:-bar()

invMetric = Matrix(%id = 36893627907715573868)

TM:-SetMetric(Matrix(3, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -2}))

"updating Metric"

TM:-foo()

Metric = Matrix(%id = 36893627907715551708)

TM:-bar()

invMetric = Matrix(%id = 36893627907715552068)

Same Matrix, no need to update

TM:-SetMetric(Matrix(3, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -2}))

TM:-foo()

Metric = Matrix(%id = 36893627907715551708)

TM:-bar()

invMetric = Matrix(%id = 36893627907715552068)

TM:-SetMetric(Matrix(3, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -1/2}))

"updating Metric"

TM:-foo()

Metric = Matrix(%id = 36893627907715532428)

TM:-bar()

invMetric = Matrix(%id = 36893627907715532788)

NULL

Download 2024-12-30_Q_Module_Global_and_Local_ac.mw

Here's a somewhat crude use of dsolve's events option to get the points where the diff(x(t),t) hits zero. It could be made fancier.

If you reduce T then it finds less data point and runs faster.

restart

Settings(typesetdot = true)

eqn1 := (M+m)*(diff(x(t), `$`(t, 2)))+c__b*(diff(x(t), t))+m*l*(diff(varphi(t), `$`(t, 2))-(diff(varphi(t), t))^2*varphi(t))+K*x(t) = f*cos(phi)

eqn2 := m*l^2*(diff(varphi(t), `$`(t, 2)))+c__h*(diff(varphi(t), t))+m*l*(diff(x(t), `$`(t, 2)))+m*g*l*varphi(t) = 0

couple := [eqn1, eqn2]

M := 300; m := 200; l := 1; c__b := 300; K := 2000; f := 200; c__h := 0; g := 9.8

couple

[500*(diff(diff(x(t), t), t))+300*(diff(x(t), t))+200*(diff(diff(varphi(t), t), t))-200*(diff(varphi(t), t))^2*varphi(t)+2000*x(t) = 200*cos(phi), 200*(diff(diff(varphi(t), t), t))+200*(diff(diff(x(t), t), t))+1960.0*varphi(t) = 0]

T := 50000; `&omega;__n` := 10

50000

10

`&omega;__v` := `&omega;__n`*t/T; phi := int(%, t)

ics5 := x(0) = 0.1e-3, (D(x))(0) = 0, varphi(0) = 0.1e-3, (D(varphi))(0) = 0

solu12:=dsolve({couple[],ics5},[x(t),varphi(t)],output = listprocedure, numeric,method=rkf45,
               events=[[diff(x(t),t), halt]],
               maxfun =-1):

interface(warnlevel=0):
str := time[real]():
TT := 'TT': ii_up:=100000:
for ii from 1 to ii_up do
  solu12(T):
  TT[ii] := eval([1*omega__n*t(last)/T,x(t)(last)],solu12(last));
  try solu12(eventclear); catch: end try;
  if abs(TT[ii][1] - omega__n) < 1e-5 then
    high_ii:=ii; printf("%ld points", high_ii); ii:=ii_up; next;
  end if;
end do:
(time[real]()-str) * 'seconds';

79632 points

103.777*seconds

data := convert(TT,list): numelems(data);
plots:-pointplot(data, symbol=solidcircle, symbolsize=5);

79632

plots:-display(
  plot([seq(data[2*i], i=1..numelems(data)/2)], color=red),
  plot([seq(data[2*i-1], i=1..numelems(data)/2)], color=red)
);

 

 

Download saopinjifen1230_ac.mw

We can also overlay that with the odeplot. Here it is for N=10000 (similar effect for larger N, but it all takes longer). Continuing in that code,

solu12(eventdisable={1}):
p1 := plots:-odeplot(solu12, [omega__v, x(t)], 0 .. T, numpoints=50000, color="Niagara Azure"):
plots:-display(
  p1,
  plot([seq(data[2*i], i=1..numelems(data)/2)], color=red),
  plot([seq(data[2*i-1], i=1..numelems(data)/2)], color=red),
  size=[700,400]
);


saopinjifen1230_ac2.mw

Perhaps you could use WorksheetToMapleText from the Worksheet package.

As an example,

currentdir("/home/harpo/skits");
read "slapstick/boguscalculus.mw";

In other words, if you've set currentdir(...) then you don't have to take the trouble to form the fully qualified location (by, say, concatenation).

That's the main point of being able to set the current working directory; it allows you to use relative paths, without having to form them into absolute paths (wrt the root directory).

ps. As a side note, if also you happen to want to obtain the home directory programmatically, you could issue the Maple command kernelopts("homedir") or if it's set as an environment variable in your OS shell, getenv("HOME") .

Does the command currentdir() return a location which is writable by you, when executed by your script approach?

If not, you could set it to a location writable by you, or try passing the full location in the first argument passed to Export. Eg,
   "C:/Users/JoyDivisionMan/LineDensity.jpg"
or what have you, depending on your OS and account details.

nb. We don't know what Maple version and OS you have.

I changed your Post into a Question.

Here is one way of correcting the several mistakes in use of the dsolve result and printf.

I retained your dsolve usage and odeplots.

 

 

restart; with(plots); _local(gamma)

sys := {diff(i(t), t) = lambda*s(t)-gamma*i(t), diff(r(t), t) = gamma*i(t), diff(s(t), t) = -lambda*s(t)}

ic := {i(0) = 1, r(0) = 0, s(0) = 999999}

gamma := .1

lambda := .2

NULL

sol := dsolve(`union`(sys, ic), numeric)

``

display([odeplot(sol, [t, s(t)], 0 .. 60, color = red), odeplot(sol, [t, i(t)], 0 .. 60, color = blue), odeplot(sol, [t, r(t)], 0 .. 60, color = green)], labels = ["Time", "Number of People"], legend = ["Susceptible", "Infected", "Recovered"])

results := map(proc (u) options operator, arrow; map(proc (uu) options operator, arrow; lhs(uu) = round(rhs(uu)) end proc, u) end proc, [seq(sol(tval), tval = 0 .. 50)])

printf("%-10s %-15s %-15s %-15s\n", "Day", "Susceptible", "Infected", "Recovered"); printf("---------------------------------------------\n"); for entry in results do printf("%-10d %-15d %-15d %-15d\n", eval(t, entry), eval(s(t), entry), eval(i(t), entry), eval(r(t), entry)) end do

Day        Susceptible     Infected        Recovered      
---------------------------------------------
0          999999          1               0              
1          818730          172214          9056           
2          670319          296822          32859          
3          548811          384014          67175          
4          449329          441982          108689         
5          367879          477303          154818         
6          301194          495235          203571         
7          246597          499977          253427         
8          201896          494865          303239         
9          165299          482542          352160         
10         135335          465088          399577         
11         110803          444136          445061         
12         90718           420952          488330         
13         74273           396516          529210         
14         60810           371574          567616         
15         49787           346686          603527         
16         40762           322269          636969         
17         33373           298620          668006         
18         27324           275950          696726         
19         22371           254396          723234         
20         18316           234039          747645         
21         14996           214922          770083         
22         12277           197052          790671         
23         10052           180414          809534         
24         8230            164976          826794         
25         6738            150694          842568         
26         5517            137514          856969         
27         4517            125378          870106         
28         3698            114224          882078         
29         3028            103991          892981         
30         2479            94617           902905         
31         2029            86040           911931         
32         1662            78201           920137         
33         1360            71046           927594         
34         1114            64519           934367         
35         912             58571           940517         
36         747             53154           946099         
37         611             48225           951164         
38         500             43741           955759         
39         410             39664           959926         
40         335             35960           963704         
41         275             32596           967129         
42         225             29541           970234         
43         184             26769           973047         
44         151             24253           975596         
45         123             21971           977905         
46         101             19902           979997         
47         83              18025           981892         
48         68              16324           983608         
49         55              14782           985162         
50         45              13385           986570         

 

Download sir_model_ac.mw

ps. It's potentially confusing to people if you use the Maple 2025.0 Beta release for your attachments here./

First 9 10 11 12 13 14 15 Last Page 11 of 336