## 9088 Reputation

17 years, 20 days

## Use Statistics[Fit]...

Here's the basic idea.
```restart;
y := a*x*sin(b*x)+c*x+d:
# create some data (X,Y) with specific parameter values
y1 := unapply(subs(a=1,b=2,c=3,d=4,y),x):
X := Vector([seq(0..10, 0.1)],datatype=float):
Y := Vector(map(y1,X),datatype=float):
# now, use Fit to extract the parameters from the data
with(Statistics):
yfit := Fit(subs(b=2,y), X,Y, x);
3.00000000016562396 x + 3.99999999924536232 + 0.99999999999148626 x sin(2 x)
y2 := unapply(yfit,x):
plot([y1,y2], 0..10);```
You'll notice I cheated and inserted the frequency term (b=2) into the expression passed to Fit. Ideally, Fit could solve for that parameter, however, I was not able to get it to return a reasonable value, despite trying various weighting schemes. However, I'm not particularly adept at using the Statistics package.

## try Matrix Browser...

You can use the Matrix Browser for this, after converting the Matrix to reals. If the magnitude of the values respresents the intensity, then map the abs function over the Matrix. For example:
```# Create generator for complex values in unit circle.
use gen = (rng -> RandomTools:-Generate(float(range=rng, method=uniform))) in
genZ := () -> evalc(polar(gen(0..1),gen(0..evalhf(2*Pi))));
end use;
# Create matrix of complex floats
img := Matrix(20,20,genZ):
# Compute magnitudes
mag := map(abs,img);```
In the standard gui, right click on the output, selecte "Browse", and then click the "Image" tab.

## A few problems....

There are a number of problems with the procedure assignment. First, the final assignment (dr = f6/2) does not do what you expect. Because dr is declared global and is the name of the procedure, it reassigns he procedure to the value f6/2. For example, consider
```tmp := proc(x) global tmp; tmp := x; end proc:
tmp(23);
23
tmp;
23
tmp(34);
23```
Instead you should just have the expression f6/2, which becomes the return value. Or use the statement return f6/2. There are other problems. You can use either sum, for symbolic summation, or add for summing over explicit values, however, summation is not an assigned Maple command (though integrate is a synonym for int). However, in your example, either choice has problems. Your inner do loop uses r as an index. As such, you won't be able to use r as in index for the final sum/add (it isn't clear why you are doing so, anyway). Because you are summing over an explicit range, you probably want to use add. However, it works differently than sum; you cannot do, for example
```s := i:
6```
and expect to get a reasonable answer. Yes, it worked in that case, however, that is because this was done at the global level. If you do the same thing in a procedure you get
```proc()
local s,i;
s := i;
end proc();
3*i;```
Giving a good explanation as to why it works at the global level is tricky; suffice it say that it is best not to use add (or seq or mul) in that fashion. The proper way is to explictly include the expression with the index in the call to add. In your case, that would be done with
`f4 := add(factorial(n+1)*((r+1-i)^n)*((-1)^i)/(factorial(i)*factorial(n+1-i)), i=0..r);`
.

## joy to ode...

I don't have Maple8 locally installed, but tested the following with Maple 9, 9.5, 10, and 11 and got the same results (note change to BC, I just set D(Phi)(1) = 0, since that's what it evaluates to). N:=2: Xi:=1: mu:=0: m:=4: l:=2: sigma:=1: Tau:=0.7: M:=1/(2*(1-(Eta)^2)): a:=Eta: Eta:=Tau*sqrt(2/3): omega:=0.9*m*Eta: g:=r->sqrt((1-2*M/r^(2*N)+2*M*a^2/r^(2*N+2))^(-1)): h:=r->sqrt(r^2*(1+2*M*a^2/r^(2*N+2))): Omega:=r->2*M*a/(r^(2*N+2)+2*M*a^2): f:=r->r/(g(r)*h(r)): V0:=r->f(r)^2*sqrt(h(r))/(r^(N+1))*diff(((f(r)^2)*h(r)/r)*diff(sqrt(h(r))*r^N,r),r): V:=r->V0(r)+f(r)^2*mu^2-(omega-m*Eta)^2+(f(r)^2/(r^2))*(l*(l+2*N)-m^2*(1-(r^2/(h(r)^2)))+4*(1-sigma)*(h(r)^2/(r^2)-1)): YY:=diff(V(r),r): XX:=eval(YY,r=1): alpha1:=diff(g(r)^(-2),r): alpha:=r->alpha1/h(r): x:=r->int(g(r)/f(r),r): ODE := -f(r)*(g(r)^(-1))*diff(f(r)*(g(r)^(-1))*diff(Phi(r),r),r)+2*I*(omega-m*Eta)*f(r)*(g(r)^(-1))*diff(Phi(r),r)+Phi(r)*(V(r)+(omega-m*Eta)^2)=0: # Subject to the BC: BC := ( Phi(1)=1, D(Phi)(1) = 0 ): # I ask Maple to solve the ODE with respect to Phi(r) as follows: dsol := dsolve({ODE,BC},numeric); dsol := proc(x_rkf45) ... end proc dsol(2); [r = 2., Phi(r) = 20.1567348648848 - 7.44101965598550 I, d -- Phi(r) = 43.5895509783580 - 21.8247802302491 I] dr

## evaluate boundary condition at r=1...

Try changing second boundary condition assignment to:
`D(Phi)(r) = subs(r=1,(YY*(r-1))/(alpha^2-2*I*alpha*(omega-m*Eta)))`
It then runs on my box:
`dsol := dsolve({ODE,BC},numeric,Phi(r)):`

## partial fix...

One problem is your method for specifying the boundary conditions. The first should be specified as
` Phi(1)=1`
The second is problematic. Are you trying to specify a boundary condition at a symbolic value, p1? I don't believe that numeric dsolve will handle that (but am not sure). Another problem is the existence of the complex values in the boundary condition. There are also problems with some of your assignments. For example
`V:=r->V0(r)+f(r)^2*mu^2-(omega-m*Eta)^2+(f(r)^2/(r^2))*(l(l+2*N)-m^2*(1-(r^2/(h(r)^2)))+4*(1-sigma)*(h(r)^2/(r^2)-1));`
is suspect, the term l(l+2*N) is probably not what you intend (it always evaluates to 2, since l is assigned 2:
``` l := 2: l(l+23+a);
2;```
Probably you meant l*(l+2*N). The error could be just the cut/paste from the standard gui, but it is usually best to use explicit multiplication.

## Symbolic indexes to rtables...

The problems are the references to A[k,j] and M[k,j]. You cannot reference an rtable (a data structure that encompasses Matrices and Vectors) with a symbolic index. For example:
```V := Vector([1,2,3]):
sum(V[k], k=1..3);
If the sum is changed to an add, then the problem goes away because an integer value is substituted for the k before it is evaluated (add works at a lower level than sum):
```add(V[k], k=1..3);
6```
In your example, changing sub to add is useful, regardless, because add is usually better for summing over a given range. However, that alone won't fix the problem because the j term is symbolic. To fix that just wrap the whole expression in your do loop:
```for k from 1 to 5 do
for i from 1 to 10 do
end do;
end do:```
Now Maple is happy. Alas, that doesn't fix all the problems. The handling of u is mostly likely wrong (I'm not sure what you intend).

## one way...

```LegendreP(1,1,cos(theta));
LegendreP(1, 1, cos(theta))
simplify(%);
(1/2)                 (1/2)
(cos(theta) - 1)      (cos(theta) + 1)
expand(%^2);
2
-1 + cos(theta)
simplify(sqrt(%),symbolic);
I sin(theta)
```

## fix relations to define endpoint...

Note that you should change the relations to
`0 ≤ x ≤ a, a < x ≤ b`
Otherwise there is an indeterminancy when x = a (which Maple resolves, but is poor form).

## function composition...

Much as I enjoy doing homework, I prefer doing my own. However, I'll give you a few hints. First, I suspect that the asterisk used in [g*f] is intended to be the composition operator, not multiplication. Use @ (see ?@), the Maple composition operator, to do this, after assigning f and g as functional operators (see ?operators,functional). Here is an example of composing a function with itself
```f := x -> x^2 + 2:
h := f@f:
h(x);
(x^2+2)^2 + 2```
Note that one can also do
```h := f*f:
h(x);
(x^2+2)^2```
so the proper answer depends on the interpretation of the asterisk.

## use plotsetup...

Use plotsetup:
```plotsetup('gif', 'plotoutput' = "/tmp/my_plot.gif"):
plot(sin, 0..Pi);
```
To do this in a loop, you might try
```for i to 5 do
plotsetup('gif', 'plotoutput' = sprintf("/tmp/my_plot%d.gif", i));
plot( f[i], ... );
end do:```
where f[i] is the expression to be plotted.

## remeasured...

Here's the comparison without the algorithmic optimizations.
```               +-------- timing -------+--------- memory ---------+
Procedure      | maple (s) | total (s) | used (MiB) | alloc (MiB) | gctimes
---------------+-----------+-----------+------------+-------------+--------
tp_original    |    22.69  |    23.08  |   246.61   |     36.43   |    6
tp_evalhf      |     2.06  |     2.16  |     0.00   |      0.00   |    0
tp_autocompile |     0.15  |     0.38  |     0.55   |      0.00   |    0
tp_compile     |     0.09  |     0.09  |     0.00   |      0.00   |    0
---------------------------------------------------------------------------
```
tp_compile looks like
`tp_compile := Compile:-Compile( proc() ... end proc):`
It doesn't run any faster than tp_autocompile, the timing difference is that the compilation time is not included in the run. To measure those accurately I'd need to execute them several times. Declaring the local counters as integers doesn't have any measurable effect (measured with multiple executions).

## try autocompile...

You can also use the autocompile option, which is frequently faste than calling evalhf. Also, there are two small optimizations you can make, not really pertinent to the broader topic. First, don't compute the square root, just compare against 100.0^2. Second, break out of the inner loop if Temp3 ≥ 100.0^2. Here are the revisions:
```testproc3 := proc()
local tp;
tp := proc()
local i, j, k, total, Temp3, Temp1, Temp2;
total := 0;
for i to 100 do
Temp1 := i*i;
for j to 100 do
Temp2 := Temp1+j*j;
for k to 100 do
Temp3 := Temp2+k*k;
if Temp3 >= 10000.0 then
break;
end if;
total := total+1;
end do
end do
end do;
total/0.1e6
end proc:
evalhf(tp());
end proc:

testproc4 := proc()
local i, j, k, total, Temp3, Temp1, Temp2;
option autocompile;

total := 0;
for i to 100 do
Temp1 := i*i;
for j to 100 do
Temp2 := Temp1+j*j;
for k to 100 do
Temp3 := Temp2+k*k;
if Temp3 >= 10000.0 then
break;
end if;
total := total+1;
end do
end do
end do;
total/0.1e6
end proc:```
With testproc1 and testproc2 being yours and Jacque's procedures, a comparison of the four gives
```            +-------- timing -------+--------- memory ---------+
Procedure   | maple (s) | total (s) | used (MiB) | alloc (MiB) | gctimes
------------+-----------+-----------+------------+-------------+--------
testproc1   |    22.97  |    24.74  |   247.12   |     37.68   |    6
testproc2   |     2.06  |     2.05  |     0.00   |      0.00   |    0
testproc3   |     1.04  |     1.03  |     0.00   |      0.00   |    0
testproc4   |     0.27  |     0.55  |     3.27   |      2.56   |    0
------------------------------------------------------------------------
```

## simpler...

Slightly simpler is just
`map(sin, [3,4,5]);`
or
`evalf(map(sin,[3,4,5]));`
Another way to do this is
`[seq(sin(x), x=[3,4,5])];`
It builds a sequence and encloses it in a list.

## post truncated (less-than html problem?)...

The rest of your post is missing. If it includes a less-than sign, substitute the html equivalent `&lt;' (not including the quotes). Note that the range in the seq assignment to datalist should be 1..n, not 0..n. Also where is plot2?
 First 103 104 105 106 107 108 109 Page 105 of 111
﻿