## 7345 Reputation

10 years, 343 days

## Apply a little thought?...

solve() *seems* to work better with equalities than inequalities. So utilise this.

restart;
eqn:=-0.005*x^2+10.4*x-1230-624*sin(0.002*x);
rts:=select(type, [solve(eqn)], realcons);

will return two real values.

Now for large (absolute) values of 'x', your equation is dominated by the quadratic term - hence its value will be negative. This information, tells you that the value of 'eqn' is positive only between the two values of 'rts'.

## Alternatively...

One could use the 'CurveFitting" package - as in

restart:
#
# Specify data
#
xdata:= [21.2, 24.7, 20.8, 20.8, 20.3]:
ydata:= [16.6, 19.7, 16.4, 16.8, 16.9]:
#
# fit above data to a straight line using
# CurveFitting package
#
fit:= CurveFitting:-LeastSquares
( xdata,
ydata,
x
);
#
# Produce plot from the supplied points
#
p1:=  plots:-pointplot( xdata,
ydata,
symbol=solidcircle,
symbolsize=16,
color=blue
):
#
# Produce plot of the straight-line fit generated
# above
#
p2:=  plot( fit,
x=min(xdata)..max(xdata),
color=red
):
#
# Display both input points and fitted line. Fit
# is pretty rubbish
#
plots:-display([p1,p2]);

The above produces the same results as the solution posted by Carl

In both solutions, the fit might be considered as pretty cr*p

## Pay attention to datatypes...

Yorr original matrix Mn__11 (by default) has datatype=anything

In your final matrix, you specify datatype=float[8]. You appear to believe that (by magic?) Maple will convert datatype=anything to datatype=float[8]

So consider the possibility that in your first matrix, you populate it with 'strings' - perfectly fine because you have datatype=anything. Then in your second matrix, you are using an initialiser, which is a bunch of strings,  but you specify that entries are of datatype-float[8]. Now exactly what would you expect to happen????

I could force the behaviour you appear to want with

restart;
Mm__11 := Matrix
( 7,
7,
[  [  -a*b*I__0,               0,              0, -a*b*I__1,              0, 0, 0 ],
[               0,  -a*b*I__0,              0,               0, -a*b*I__1, 0, 0 ],
[               0,                0, -a*b*I__0,              0,               0, 0, 0 ],
[  -a*b*I__1,               0,               0, -a*b*I__2,              0, 0, 0 ],
[               0,  -a*b*I__1,               0,              0, -a*b*I__2, 0, 0 ],
[               0,               0,                0,              0,               0, 0, 0 ],
[               0,               0,                0,              0,               0, 0, 0 ]
]
):
#
# With the above definition, note that the stored
# datatype is "anything"
#
op([3,1], Mm__11);
#
# Now assign values
#
I__0:= 0.102033222998e-2:
I__1:= .80922342232909:
a:= .2223212330091:
b:= .293939202032292:
I__2:= 1:
#
# Define a new matrix, but since the original matrix
# ie Mm__11 has datatype "anything", it cannot be
# used as initialiser for a matrix with datatype
# float[8], unless an explicit type conversion is
# performed, as in the following
#
M__11 := Matrix
( convert~(Mm__11, float[8]),
datatype = float[8]
);

but this definition seems way too complicated - since I could achieve the same result with

restart;
I__0:= 0.102033222998e-2:
I__1:= .80922342232909:
a:= .2223212330091:
b:= .293939202032292:
I__2:= 1:
M := Matrix
( 7,
7,
[  [  -a*b*I__0,              0,              0, -a*b*I__1,              0, 0, 0 ],
[               0, -a*b*I__0,              0,              0, -a*b*I__1, 0, 0 ],
[              0,               0, -a*b*I__0,              0,              0, 0, 0 ],
[ -a*b*I__1,               0,              0, -a*b*I__2,              0, 0, 0 ],
[              0, -a*b*I__1,               0,              0, -a*b*I__2, 0, 0 ],
[              0,              0,               0,              0,              0, 0, 0 ],
[              0,              0,               0,              0,              0, 0, 0 ]
]
);

Am I missing something?

## hmmmm...

Well it is not generally recommended to do this, but

interface(warnlevel=0)

will suppress all warnings. See the help at ?interface

(and remember to turn them back  on again as soon as convenient!)

## One way...

but maybe not the best?

The following works for me: Change the file path/name to something appropriate for you

restart;
A:=Array():
k:=1:
while line <>0 do
L:= remove
( type,
StringTools:-Split( line, " "),
""
);
A(k):=[ StringTools:-Join( L[1..-3] ),
parse(L[-2]),
parse(L[-1])
];
k:=k+1;
od:
fclose(fid);
M:= Matrix(k-1,3, (i,j) -> A[i][j] );

## No such (guaranteed) thing...

Truly global optimisation for arbitrary functions is never guaranteed to be correct - unless there has been some major breakthrough in optimization theory which I haven't heard about.

Global optimization of convex functions is relatively trivial - but global optimization of non-convex functions.........?

Local optimization of non-convex functions is also trivial - but global......?

I suggest that you either stick with the commands in Maple's Optimization package (inbuilt and thus free), or download/install the DirectSearch package (again free) from the Maple Application Centre. The latter seems to have a few more in-built heuristics for checking global optimality

## Explanation (but no solution)...

I'm using UK english setup, so it is actually quite difficult for me the even enter accented  characters. However if one just considers the letter 'e' with an acute accent.

Maple seems to return the 'length' of this character as 2. I suspect that the reason for this is that if one checks the 'byte' representation of this character using

convert("é", bytes);

This returns the two-byte list

[195, 169]

If I convert this to hexadecimal, with

map(convert, [195, 169], hex) I get [C3, A9]

Now checking the (Unicode) UTF-8 encoding scheme, one finds that the two-byte sequence [C3, A9] corresponds to UTF Code Point U+00E9 which is identified as " LATIN SMALL LETTER E WITH ACUTE "

So I can only conclude that Maple is using Uniicode, and in UTF-8, any character not in the 'old' 7-bit ASCII, is encoded as a multibyte sequence - and Maple's 'length' command is just counting the number of bytes!

It is relatively easy to plot Phi(t) versus t without "discretizing it". See the attached

odeSol.mw

The function Phi seems to be periodii, with period ~8. You want

"I would like to plot the probability density function of a state variable obtained from solving differential equations"

The existence of a probability density function, implies the existence of a random variable. What is you random variable and what is its relation to the solution of your ODE

## Basic theory...

You have a second order ODE. The general solution therefore contains two arbitrary constants. These arbitrary constants can be eliminated if you provide two boundary/initial conditions

## Hmmmmm...

Well based on this question (and your previous one), I can't help think that you are in a deep hole and all that I am doing is assisting with the digging!

However if you really, really, really want to do this, then

varnames:=['S','T','U','V','W'];
diff(f(map(entries, varnames, nolist) ), U[0,0]);

ought to work

## Nonsense...

Rr[x,y] refers to the x-th, y-th entry of the indexable entity Rr.

But it seems that Rr is a (non-indexable) procedure which accepts two arguments!!!!!

So what can the OP mean - once again reduced to guesswork!!!

Perhaps what is required is to apply the procedure Rr() to all entries in an indexable entity - as in the following  'toy' example

A:=Matrix(2,2,[[1,2],[3,4]]);
B:=Matrix(2,2,[[5,6],[7,8]]);
Tr:=proc(x)
return x^2:
end proc;
#
# One way
#
Rr~(A);
Rr~(B);
#
# Another way
#
map(Rr, A);
map(Rr, B);

On the other hand perhaps what is required is a procedure which takes two values, and can be mapped across two indexable entities as in

restart;
A:=Matrix(2,2,[[1,2],[3,4]]);
B:=Matrix(2,2,[[5,6],[7,8]]);
Rr:= proc(x,y)
return x^2+y^2:
end proc;
zip(Rr, A, A);
zip(Rr, B, B);
zip(Rr, A, B);

Who knows?

## Several ways...

The optimum choice probably depends on what you plan on doing with the columns, once you have them.

with(LinearAlgebra):
#
# generate a random amtrix, just for illustration
#
M:=RandomMatrix(10,10);
#
# Generate a list of columns in a few different ways
#
Column(M, 1..10);
seq( M[..,j], j=1..10);
for j from 1 by 1 to 10 do
M[..,j];
od;

## Well...

Don't undertstand the 'for each' construct - but you don't need it when definiing a vector with an embedded function initialiser.

You do however need an explicit value for 'N' in order to define a Vector()

The initialiser function within the Vector definition should be '->' not '-->'

The following works for me

N:=10;
G1:=Vector[column](N,i->f(i,1));
G2:=Vector[column](N,i->f(i,2));

or if you are planning on defining a lot of vectors , then you could use

N:=10;
for j from 1 by 1 to 2 do
G||j:=Vector[column](N,i->f(i,j));
od;

## You are doing it the difficult way!!!!...

Since your orginal post stated that you were " simulating the movement of 4 astronomical bodies with Newtons equation", I guessed that this involved solving a system of ODEs numerically. Technically dsolve() returns a procedure (or a list of procedures, or an Array, depending on the options you set in the dsolve() command. Values for all problem variables can be be obtained by appropriate calls to this procedure.

I have added a few execution groups to your worksheet which performs these calls and plots the relevant spacecurves. I have deliberately done this in a step-by-step manner to illustrate the process. I could have achieved the same functionality more efficiently (and with less code to type!), but it would probably be more difficult for you to understand.

Check out the attached

planets.mw

## Quick (and dirty) tip...

When solve() produces some horrendously long expressions, making it difficult to figure out whihc one(s) might be interesting, it is sometimes useful to use complexplot(), as in

restart;
eqn:= (1/4)*x*(3-2*x)^2/(1-x)^3=18325.73901+exp(36.58420421-10902.09286/T);
plots:-complexplot( [solve(eqn,x)],
T=300..800,
axes=boxed
);

which will plot all the solutions - real and complex. From the output of the above it is farly obvious that over the specified range of the cariable 'T',  one solution lies on the real axis and two don't. If you know Maple's default colour choice for such a plot (darkRed, blue, green) then you can tell immediately that the second solution returned by solve() is the real one. So if the real solution is the one you want then

plot( solve(eqn,x)[2],T=300..800);

will plot it.

You can then delete the complexplot() command from your worksheet to leave

restart;
eqn:= (1/4)*x*(3-2*x)^2/(1-x)^3=18325.73901+exp(36.58420421-10902.09286/T);
plot( solve(eqn,x)[2],T=300..800);

and, hey presto you are a genius!

 First 111 112 113 114 115 116 117 Last Page 113 of 148
﻿