pagan

5062 Reputation

23 Badges

14 years, 184 days

 

 

"A map that tried to pin down a sheep trail was just credible,

 but it was an optimistic map that tried to fix a the path made by the wind,

 or a path made across the grass by the shadow of flying birds."

                                                                 - _A Walk through H_, Peter Greenaway

 

MaplePrimes Activity


These are replies submitted by pagan

@PatrickT I suspect some misunderstanding.

Firstly, in your earlier versions you made assignments such as sol:=dsol along with the comment that this produced a copy of dsol. But that is not so. That assignment does not make a copy of dsol, or of the three procs within dsol. (More on copying, below...)

Even if you assign the value of (read: the procedure assigned to) dsol to various other names, you still have to reset the initial conditions (over again) of the earlier names each time that you wish to go back and recalculate with them.

I don't see any cheap way around having to reset the initial conditions each time you wish to recalculate. You don't have to reinvoke dsolve, I think, but you do have to reinstate the various initial conditions. For example,

restart:

dsol :=
  dsolve(
      [ diff(u(t),t)=u(t)-v(t), diff(v(t),t)=u(t)-w(t), diff(w(t),t)=v(t)-w(t)
      , u(0)=2, v(0)=1, w(0)=2 ]
      , [u(t),v(t),w(t)]
    , 'type' = numeric
    , 'output' = listprocedure
  ):

solA:=dsol:
solA(initial=[2,1,2]):
solA(-135.0);

[t(-135.0) = -135.0, u(t)(-135.0) = HFloat(3.9077543506678216), 

  v(t)(-135.0) = HFloat(4.992270555376236), 

  w(t)(-135.0) = HFloat(4.0845162047084)]

solB:=dsol:
solB(initial=[6,7,8]):
solB(-135.0);

 [t(-135.0) = -135.0, u(t)(-135.0) = HFloat(8.08460488495403), 

   v(t)(-135.0) = HFloat(7.17680411000912), 

   w(t)(-135.0) = HFloat(6.092199225055099)]

solA(initial=[2,1,2]):
solA(-135.0);

[t(-135.0) = -135.0, u(t)(-135.0) = HFloat(3.9077543506678216), 

  v(t)(-135.0) = HFloat(4.992270555376236), 

  w(t)(-135.0) = HFloat(4.0845162047084)]

Are you seeing behaviour where this kind of approach is not good enough, and still has "contamination"?

In light of the above, there's not much point to having different names solA, solB, etc, with this particular scheme above. But... see at end for an alternative which does merit the various names.

I don't think that you can get around this calling the `copy` command on each of the three procs within dsol either. I mean, by doing solAu:=copy(eval(u(t),dsol)), etc, or with various judiciously placed calls like forget(solB), etc. And that's not only because of the way Maple uniquifies procs (ie. "simpls them"). You can subsop(1=..,eval(u(t),dsol)) and get various instances with actual different addresses, and still have same behaviour.

If you want to have various named procedure, for convenience of reexecution, then you can write them as simple procedures which wrap around preliminary calls to reset the initial conditions. Eg.

solA:=proc(tt) dsol(initial=[2,1,2]); dsol(tt); end proc:

solB:=proc(tt) dsol(initial=[6,7,8]); dsol(tt); end proc:

solA(-135.0);

[t(-135.0) = -135.0, u(t)(-135.0) = HFloat(3.9077543506678216), 

  v(t)(-135.0) = HFloat(4.992270555376236), 

  w(t)(-135.0) = HFloat(4.0845162047084)]

solB(-135.0);

 [t(-135.0) = -135.0, u(t)(-135.0) = HFloat(8.08460488495403), 

   v(t)(-135.0) = HFloat(7.17680411000912), 

   w(t)(-135.0) = HFloat(6.092199225055099)]

solA(-135.0);

[t(-135.0) = -135.0, u(t)(-135.0) = HFloat(3.9077543506678216), 

  v(t)(-135.0) = HFloat(4.992270555376236), 

  w(t)(-135.0) = HFloat(4.0845162047084)]

Hopefully this may be of some use to you.

@Joe Riel Yes, unapply(..,numeric) is nice and terse. (And, it's not necessary to have both y00 and y0. A single, appropriate proc would do. I figured that he wanted it clearly laid out, and might eventually use a more complicated example.)

restart:
> y00 := unapply(t,t,'numeric'):                                              
> y0:=unapply(y00(t-1),t,'numeric'):                                       
> p:=dsolve({D(x)(t)=y0(t), x(t0)=x0}, x(t), numeric, parameters=[t0, x0]):
> p(parameters=[0,y0(0)]);                                                 
                              [t0 = 0., x0 = -1.]

> p(4);                                                                    
                       [t = 4., x(t) = 3.00000000000000]

> p(7.5);
                      [t = 7.5, x(t) = 19.6250000000000]

Is it clear why the 'known' option wasn't needed, in order to use y0(t) in the ode system?

And output=listprocedure could be used in order to avoid use of op([2,2],..) in his second use of dsolve.

> p:=dsolve({D(x)(t)=y0(t), x(t0)=x0}, x(t), numeric,
>           parameters=[t0, x0], output=listprocedure):
> X:=eval(x(t),p):        
> X(parameters=[0,y0(0)]);
                              [t0 = 0., x0 = -1.]

> X(7.5);                 
                               19.6250000000000

Etc.

@Joe Riel Yes, unapply(..,numeric) is nice and terse. (And, it's not necessary to have both y00 and y0. A single, appropriate proc would do. I figured that he wanted it clearly laid out, and might eventually use a more complicated example.)

restart:
> y00 := unapply(t,t,'numeric'):                                              
> y0:=unapply(y00(t-1),t,'numeric'):                                       
> p:=dsolve({D(x)(t)=y0(t), x(t0)=x0}, x(t), numeric, parameters=[t0, x0]):
> p(parameters=[0,y0(0)]);                                                 
                              [t0 = 0., x0 = -1.]

> p(4);                                                                    
                       [t = 4., x(t) = 3.00000000000000]

> p(7.5);
                      [t = 7.5, x(t) = 19.6250000000000]

Is it clear why the 'known' option wasn't needed, in order to use y0(t) in the ode system?

And output=listprocedure could be used in order to avoid use of op([2,2],..) in his second use of dsolve.

> p:=dsolve({D(x)(t)=y0(t), x(t0)=x0}, x(t), numeric,
>           parameters=[t0, x0], output=listprocedure):
> X:=eval(x(t),p):        
> X(parameters=[0,y0(0)]);
                              [t0 = 0., x0 = -1.]

> X(7.5);                 
                               19.6250000000000

Etc.

Thx for changing the title from "Prameter" to "Parameter".

This is one of the archetypal (avoidable) inefficient programming techniques in Maple -- to build a list by concatenation which is a worse complexity class than using `seq`.

The significant consequence is that the timing disadvantage is a function of the number of elements created. The performance of list concatenation gets even worse relative to that of using `seq` as the problem size increases.

restart:
st:=time():
  K:=[seq(i^2,i=1..10^4)]:
time()-st;
                               0.

st :=time():
L := []:
for j to 10^4 do L := [op(L), j^2] end do:
time()-st;
                             0.960

evalb(K=L);
                              true

restart:
st:=time():
  K:=[seq(i^2,i=1..10^5)]:
time()-st;
                             0.020

restart:
st :=time():
L := []:
for j to 10^5 do L := [op(L), j^2] end do:
time()-st;
                             93.780

Above, the timing of list concatenation got worse by a factor of approximately 100 when the problem size was increased by a factor of 10.

But `seq` slows down by not too much more than the factor by which the problem size increases (up until much larger sizes, at least).

restart:
st:=time():
  K:=[seq(i^2,i=1..10^6)]:
time()-st;
                             0.170

restart:
st:=time():
  K:=[seq(i^2,i=1..10^7)]:
time()-st;
                             1.850

The inefficiency is not only for timing. It is also for the total memory allocation.

This is one of the archetypal (avoidable) inefficient programming techniques in Maple -- to build a list by concatenation which is a worse complexity class than using `seq`.

The significant consequence is that the timing disadvantage is a function of the number of elements created. The performance of list concatenation gets even worse relative to that of using `seq` as the problem size increases.

restart:
st:=time():
  K:=[seq(i^2,i=1..10^4)]:
time()-st;
                               0.

st :=time():
L := []:
for j to 10^4 do L := [op(L), j^2] end do:
time()-st;
                             0.960

evalb(K=L);
                              true

restart:
st:=time():
  K:=[seq(i^2,i=1..10^5)]:
time()-st;
                             0.020

restart:
st :=time():
L := []:
for j to 10^5 do L := [op(L), j^2] end do:
time()-st;
                             93.780

Above, the timing of list concatenation got worse by a factor of approximately 100 when the problem size was increased by a factor of 10.

But `seq` slows down by not too much more than the factor by which the problem size increases (up until much larger sizes, at least).

restart:
st:=time():
  K:=[seq(i^2,i=1..10^6)]:
time()-st;
                             0.170

restart:
st:=time():
  K:=[seq(i^2,i=1..10^7)]:
time()-st;
                             1.850

The inefficiency is not only for timing. It is also for the total memory allocation.

Using `solve` instead of `fsolve` is inefficient for finding approximations of roots of univariate polynomials which Maple cannot factor explicitly.

The reason is explained here.

Basically, calling evalf on the indexed RootOfs will involve a call to fsolve which finds all the roots. This will be done for each indexed RootOf that gets evalf'd. In the worst case using `solve` instead of `fsolve` for such a task can involve fsolve being called n times instead of just once, where n is the degree of the univariate polynomial.

Some timings:

restart:
Digits:=12:
st:=time():
for i from 1 to 100 do
  forget(evalf):
  evalf([solve(z^5-5*z^4-5*z-5,z)]);
end do:
time()-st;
                             1.329

restart:
Digits:=12:
st:=time():
for i from 1 to 100 do
  forget(evalf):
  [fsolve(z^5-5*z^4-5*z-5,z,complex)];
end do:
time()-st;
                             0.169

restart:
Digits:=12:
st:=time():
for i from 1 to 100 do
  forget(evalf):
  [solve(z^5-5*z^4-5*z-5.0,z)];
end do:
time()-st;
                             1.611

The results from fsolve are also more accurate than the result from solve above.

Using `solve` instead of `fsolve` is inefficient for finding approximations of roots of univariate polynomials which Maple cannot factor explicitly.

The reason is explained here.

Basically, calling evalf on the indexed RootOfs will involve a call to fsolve which finds all the roots. This will be done for each indexed RootOf that gets evalf'd. In the worst case using `solve` instead of `fsolve` for such a task can involve fsolve being called n times instead of just once, where n is the degree of the univariate polynomial.

Some timings:

restart:
Digits:=12:
st:=time():
for i from 1 to 100 do
  forget(evalf):
  evalf([solve(z^5-5*z^4-5*z-5,z)]);
end do:
time()-st;
                             1.329

restart:
Digits:=12:
st:=time():
for i from 1 to 100 do
  forget(evalf):
  [fsolve(z^5-5*z^4-5*z-5,z,complex)];
end do:
time()-st;
                             0.169

restart:
Digits:=12:
st:=time():
for i from 1 to 100 do
  forget(evalf):
  [solve(z^5-5*z^4-5*z-5.0,z)];
end do:
time()-st;
                             1.611

The results from fsolve are also more accurate than the result from solve above.

There is not yet an attachment to this post.

@serena88 If your "data" is just an expression sequence, or a list or a set, then you can use the `nops` command to get the number of entries. If the "data" is just an expression sequence (ie. items in a sequence, printed with commas but otherwise not encapsulated in another structure) then you can use `nops` after wrapping it in [..] to make a list.

> s:=3,5,7,9;               
                                s := 3, 5, 7, 9

> L:=[s];                   
                               L := [3, 5, 7, 9]

> S:={s};                   
                               S := {3, 5, 7, 9}

> nops(L);
                                      4

> nops(S);         
                                      4

> V:=Vector[row]([s]);
                               V := [3, 5, 7, 9]

> op(1,V);
                                       4

In Maple 15 and 16 you can also use the `numelems` command to get the number of entries of a variety of data structures. This is nice, in that it's just one command to remember (as opposed to `op`, `nops`, `Dimension`, etc).

> numelems(L), numelems(S), numelems(V);
                                    4, 4, 4

> T:='T':
> T[1]:=3: T[2]:=5: T[3]:=7: T[4]:=9:
> eval(T);
                      table([1 = 3, 2 = 5, 3 = 7, 4 = 9])

> numelems(T);
                                       4

> M:=Matrix([[s]]);      
                            M := [3    5    7    9]

> numelems(M);
                                       4

@serena88 If your "data" is just an expression sequence, or a list or a set, then you can use the `nops` command to get the number of entries. If the "data" is just an expression sequence (ie. items in a sequence, printed with commas but otherwise not encapsulated in another structure) then you can use `nops` after wrapping it in [..] to make a list.

> s:=3,5,7,9;               
                                s := 3, 5, 7, 9

> L:=[s];                   
                               L := [3, 5, 7, 9]

> S:={s};                   
                               S := {3, 5, 7, 9}

> nops(L);
                                      4

> nops(S);         
                                      4

> V:=Vector[row]([s]);
                               V := [3, 5, 7, 9]

> op(1,V);
                                       4

In Maple 15 and 16 you can also use the `numelems` command to get the number of entries of a variety of data structures. This is nice, in that it's just one command to remember (as opposed to `op`, `nops`, `Dimension`, etc).

> numelems(L), numelems(S), numelems(V);
                                    4, 4, 4

> T:='T':
> T[1]:=3: T[2]:=5: T[3]:=7: T[4]:=9:
> eval(T);
                      table([1 = 3, 2 = 5, 3 = 7, 4 = 9])

> numelems(T);
                                       4

> M:=Matrix([[s]]);      
                            M := [3    5    7    9]

> numelems(M);
                                       4

@Markiyan Hirnyk The fact that ?Table brings up the help page for the `table` command is because Table is a help alias for table. It doesn't mean anything more that that. You just see (Table) in the Help browser because of this. A help alias is a convenience, to make help lookups less unfriendly. You won't see Table describe as being equivalent to the command itself, in the actual help page.

It's like Solve. You can query ?Solve. But there's no top-level command Solve related to solve. There's a help alias, so that ?Solve produces what might be a useful help query hit.

There is no Table command or structure, just as there is no Solve command. The help system just recongnizes both as mere aliases.

Or were you joking, here?

@Markiyan Hirnyk The fact that ?Table brings up the help page for the `table` command is because Table is a help alias for table. It doesn't mean anything more that that. You just see (Table) in the Help browser because of this. A help alias is a convenience, to make help lookups less unfriendly. You won't see Table describe as being equivalent to the command itself, in the actual help page.

It's like Solve. You can query ?Solve. But there's no top-level command Solve related to solve. There's a help alias, so that ?Solve produces what might be a useful help query hit.

There is no Table command or structure, just as there is no Solve command. The help system just recongnizes both as mere aliases.

Or were you joking, here?

I like the new feature as it improves some plots whose default view was bad, for many previous releases. A typical example is tan(x), where the presence of the asymptotes caused the displayed y-range to be far too great.

If you really want backwards compatibility so that all your old worksheets produce the same (in this view regard) results from the `plot` command then you could try something drastic so as to reverse the default behavior.

The attached worksheet has been run in Maple 16.

backwards.mw

[edit note: I changed the word "fixes" to "improves", and "wrong" to "bad", in the hope that fewer people will misinterpret my opinion.]

The is a consequence of the new Smarter View of Maple 16.

Is this a workaround?

plot(exp(x),x=-3..3,smartview=false);
3 4 5 6 7 8 9 Last Page 5 of 81