Carl Love

Carl Love

28070 Reputation

25 Badges

13 years, 30 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

The first argument to animate needs to be the name of a plot command. Otherwise, how would it know that you want to animate a tube plot rather than, say, a space curve? The second argument is a list containing the arguments to that plot command, including the animation parameter. The third argument is a declaration of the animation parameter and its range. Then there are optional options, and options that apply to the plot as a whole.

plots:-animate(
     plots:-tubeplot,
     [[cos(t+phi), sin(t+phi), t], t= 0..6*Pi, numpoints=200],
     phi= 0..2*Pi,
     frames= 20,
     orientation= [40,70], shading= z, axes= box, style= patchnogrid
);

 

There are a vast array of gridline options described on the help pages ?plot,axis and ?plot,tickmarks. You can control the number of gridlines and their spacing (also color and linestyle). The spacing does not need to be even; you can place them at precise axis intersections. You can control both the major and minor gridlines.

Dumb example:

plot(0, -2..2, axis[1]= [gridlines= 20, color= red], axis[2]= [gridlines= 15, color= green]);

If after reading those help pages you still need help, please ask further.

There is no doubt a vast literature on this very important problem of finding the connected components of a graph. However, without researching any of that literature, I just now came up with this queue-based algorithm. So, there is good chance that there's a faster algorithm. But the below works. (There are also the existing commands GraphTheory:-ConnectedComponents and networks[components] that do the job, but I assumed that the goal was to write your own.)

You may well ask Since the order doesn't matter, why use a queue when you could just add/subtract items one at a time from a set? Good question. The answer is that the latter are very expensive operations in Maple; each add/subtract requires the reconstruction of the entire set. Queues (which are implemented as tables) work much faster.

I urge you not to use networks; it has been deprecated in favor of GraphTheory.

My procedure processes your table-based graph representation and returns a partition of its vertices into the connected components. I leave it to you to convert the blocks of this partition back into table-based graphs, which is trivial. If you have any trouble with that, let me know. My procedure can be easily modified to work with the GraphTheory representation of graphs. Let me know if you want that. Despite the title of your Question, my procedure is not limited to two components.

#Given a graph, return a partition of its vertices into the
#connected components.

ConnectedComponents:= proc(G::table)
local
     Q, #SimpleQueue of vertices to be visited in current component
     Comps:= table(), #set of components to return
     C,  #component currently being built
     v, #current vertex
     V:= {indices(G, 'nolist')}, #vertices remaining to be visited
     k #component index
;
     for k while V <> {} do
          #new component
          C:= table([V[1]=V[1]]);
          Q:= SimpleQueue(G[V[1]][]); #Queue the neighbors of new vertex.
          while not Q:-empty() do
               v:= Q:-dequeue();
               if assigned(C[v]) then next end if; #We already visited it.
               C[v]:= v; #Add this one to the current component.
               Q:-enqueue~(G[v]) #Queue all its neighbors.
          end do;
          #Convert to a set.
          C:= {indices(C, 'nolist')};
          V:= V minus C;
          Comps[k]:= C
     end do;
     {entries(Comps, 'nolist')}
end proc:      

ConnectedComponents(T);


 

Change the line

from i = 1 to N do

to

for i from 1 to N do

A dsolve(..., numeric) solution can be obtained very quickly. Just put the option numeric at the end of your dsolve command. Remove the evalf. To plot the solutions, do

plots:-odeplot(coz1, [x,X11(x)]);

and likewise for X22 and X33.

Let's define the error, e, as the absolute value of the difference between the true value and the computed approximation. Then for an iterative method we have a sequence of errors e[1], e[2], .... We say that the order of convergence is q if limit(e[n+1]/e[n]^q, n= infinity) approaches a positive constant u. (See the Wikipedia article "Rate of convergence".) It follows that for sufficiently large n, ln(e[n+1]) ~ ln(u) + q*ln(e[n]), where ~ denotes approximate equality. This is a linear equation in q, so we can use the e as the data for linear regression to approximate the "slope" q.


restart:

Digits:= 999:


#Given f, construct the Halley's method iterator:
Halley:= (f::algebraic, x::name)->
     unapply(simplify((d-> x - 2*f*d/(2*d^2 - f*diff(d,x)))(diff(f,x))), x)
:

#Given the iterator, a starting value, and a number of terms,
#construct the sequence of iterates.
Iterate:= proc(M::procedure, x0::complexcons, n::posint)
local k, S:= table([0=x0]);
     for k to n do S[k]:= evalf(M(S[k-1])) end do;
     convert(S, list)
end proc:

OrderOfConvergence:= (S::list(complexcons), x::complexcons)->
     (E-> diff(Statistics:-LinearFit([1,_], <E[..-2][]>, <E[2..][]>, _), _))((r-> evalf(ln(abs(r-x))))~(S))
:
 

f:= x*(x-1)^3:

M:= Halley(f,x);

proc (x) options operator, arrow; 6*x^3/(10*x^2-5*x+1) end proc

(1)

S1:= Iterate(M, 1.5, 42):  evalf[15]~(%);

[1.5, 1.26562500000000, 1.13786584182159, 1.07040105186645, 1.03559839406078, 1.01790287370617, 1.00897790388253, 1.00449563846931, 1.00224949966535, 1.00112517104490, 1.00056269096386, 1.00028137185960, 1.00014069252638, 1.00007034791261, 1.00003517436869, 1.00001758728745, 1.00000879366950, 1.00000439684119, 1.00000219842221, 1.00000109921151, 1.00000054960585, 1.00000027480295, 1.00000013740148, 1.00000006870074, 1.00000003435037, 1.00000001717519, 1.00000000858759, 1.00000000429380, 1.00000000214690, 1.00000000107345, 1.00000000053672, 1.00000000026836, 1.00000000013418, 1.00000000006709, 1.00000000003355, 1.00000000001677, 1.00000000000839, 1.00000000000419, 1.00000000000210, 1.00000000000105, 1.00000000000052, 1.00000000000026, 1.00000000000013]

(2)

That appears to be converging to the root at 1.

S2:= Iterate(M, -2, 42):  evalf[15]~(%);

[-2., -.941176470588235, -.343465682780596, -0.623834240219850e-1, -0.107834289178647e-2, -0.748310028149521e-8, -0.251417746805083e-23, -0.953540272808111e-70, -0.520197632139009e-209, -0.844610280780509e-627, -0.361510021297338e-1880, -0.283473373805794e-5640, -0.136674680739209e-16920, -0.153184722669639e-50761, -0.215673911246113e-152284, -0.601927372738376e-456853, -0.130852953777237e-1370558, -0.134431746828748e-4111675, -0.145766181215473e-12335026, -0.185832461374802e-37005079, -0.385048989544428e-111015238, -0.342530473186547e-333045714, -0.241128693989391e-999137142, -0.841197420201255e-2997411427, -0.357145388407826e-8992234280, -0.273329426779275e-26976702840, -0.122520968763409e-80930108520, -0.110352586446760e-242790325561, -0.806303973916063e-728370976685, -0.314519553989061e-2185112930054, -0.186678457533614e-6555338790162, -0.390331734766595e-19666016370487, -0.356822996199671e-58998049111461, -0.272589898063260e-176994147334383, -0.121529167269990e-530982442003149, -0.107694324522599e-1592947326009448, -0.749427629329443e-4778841978028346, -0.252545916854603e-14336525934085037, -0.966434230893684e-43009577802255112, -0.541586914806003e-129028733406765335, -0.953137900274108e-387086200220296005, -0.519539374908538e-1161258600660888014, -0.841408030930880e-3483775801982664042]

(3)

That appears to be converging to the root at 0.

Digits:= 15:

O1:= OrderOfConvergence(S1, 1);

HFloat(1.0006733239261278)

(4)

O2:= OrderOfConvergence(S2, 0);

Warning, model is not of full rank

 

HFloat(2.999999999999998)

(5)

 

``

The warning is due to several of iterates having order of convergence exactly equal to 3 (as far as floating-point computation is concerned). The warning can be ignored.

Download Halley.mw

Welcome. I moved your Question from the Posts category to the Questions category.

Try typing this. Type it explicitly; do not use palettes or menus.

int(exp(2*x), x);

Using the hint, the second initial value is obviously 21. Indeed, I don't see that the syntax allows you to call Roots(..., method= secant) or Secant without specifying a second initial value. So, I don't know what "secant method" command you are currently using. You should also lower your tolerance. Here's my worksheet:

restart:
Digits:= 18: #Close to double precision.
f:= mul(x-k, k= 1..20) - 1e-8*x^19:
r1:= Student:-NumericalAnalysis:-Roots(
     f, x= [20,21],
     method= secant,
     maxiterations= 2^10, tolerance= 10^(2-Digits)
);


#Check residual
eval(f, x=r1);
      17.9

#Compare with fsolve results.
R:= [fsolve(f, complex, fulldigits)]:
#Find fsolve root closest to r1:
Min:= infinity:
for r in R do
     m:= abs(r-r1);
     if m < Min then  r2:= r; Min:= m  end if
end do:
 
r2;
      20.2402751265279542


eval(f, x= r2);

#So, the secant method beat fsolve for accuracy!

abs(r1-r2);

Download Roots.mw

I thought that this was already answered in the past week. Your function has no discontinuities, which is obvious if you plot it. If you don't understand this, then you need to study continuity more. I'd be happy to explain it. (There is a command discont for finding discontinuities, but it says that it cannot work on this function. I think that it is confounded by the max.) Your function does however have points of nondifferentiability. Both and diff will find these without needing further prompting. These are very powerful symbolic commands.

f:= max(x^2, sqrt(abs(x))):
diff(f,x);

We'll need to take this step by step. The first thing is to correct the obvious syntax errors. You have

legend= sprintf("Total dose...

but you did not close the quotation marks.

People here would be more likely to debug your code if you used 1D input. Most of us OGs find working with 2D input very tiresome. 

Generally, the plots, if simple, are printed in the reverse order that they are listed. So, try

plots:-display([pts2, pts1, seq(p1[i], i= 3..1, -1)]);

However, this is confounded when the plots contain different types of material. The types are surfaces, curves, points, text, and axes. These types are printed in that order, and then the printing order within each type is the reverse of the order that they appear in the display. These rules apply to inline plots. Directly exported plots may be different; I'll need to research that.

It can be done using just columnwise indexing, like this:

(x,y):= 'lista[..,k]' $ k= 1..2;

or

(x,y):= seq(lista[..,k], k= 1..2);

No, your approach is not correct. Your limits of integration are wrong. [1/2 <= x*y <= 2, 1 <= x <= 3] means that the limits of integration should be [y= 1/(2*x)..2/x, x= 1..3]. The integral should be entered as

int(exp(1/(x*y))/(y^2*(x+1)^2), [y= 1/(2*x)..2/x, x= 1..3]);

Note that you should enter this using a single int to avoid an irrelevant warning message.

Note that 1/x*y is not the same as 1/(x*y). Indeed, the former is the same as y/x and the latter is the same as 1/x/y or 1/y/x. So I made it exp(1/(x*y)) because the other forms are not elementary integrals. This parentheses placement is just a guess on my part; you'll need to check it.

Maple's differential operator, D, can find the derivative of quite complicated procedurally defined functions.

f:=x->max(x^2, sqrt(abs(x))):
D(f);

eqn:= diff(y(x),x$4)-20*diff(y(x),x$3)+158*diff(y(x),x$2)-580*diff(y(x),x)+841*y(x);
eval(eqn, y(x)= x*exp(5*x)*cos(2*x));

The above returns 0, which verifies the solution.

First 276 277 278 279 280 281 282 Last Page 278 of 395