tomleslie

4621 Reputation

10 Badges

9 years, 124 days

MaplePrimes Activity


These are answers submitted by tomleslie

Somewhat tongue-in-cheek!

Utilise the debugger

restart;
Test:=proc()
plot(x^j,x=-2..2):
end proc:
stopat(Test); #set watchpoint
for j from 1 to 3 do
    Test()
end do;

The first time "Test" is called the debugger window will appear - push it out of the way somewhere then just keep hitting "continue" to generate plots one at a time

Having slept on this problem (cos it bothered me) I came back to it

With a carefully sequenced set of substituions I can transform your system of two equations (one second-order) and three boundary conditions, into a system of two first-order equations and two boundary conditions, where the first-order equation can be written in the form suitable for application of the basic rkf45 method. The process is shown in the attached

rkf45.mw

I have produced this worksheet in a very pedantic manner, with lots of comments which explain exactly what is being done at each stage - the whole transformation process could be written much more succinctly, but I thought this was a case where clarity trumps brevity.

The only assumption I had to make was that x0>0.

The final few execution blocks in the worksheet show the production of the rkf45-based dsolve() procedure, and a couple of sample plots which can be produced from it. For comparison I have included the production of the rkf45_dae-based dsolve() procedure for your original equations/bcs, along with the plots which can be obtained from it

Somwhat interestingly the"untransformed", rkf45_dae-based procedure hits numerical problems at x~=2, whereas that based on rkf45 using the transformed equations does not

 

Type or (copy paste from maple) the above into ypur faourite text editor. If using copy/paste, amke sure ther are no extraneous characters (such as the maple prompt >). You are aiming for plain text.

Save it to somewhere meaningful (depending on where you park you maple work, you might want to set up a mapleUtilities subdirectory) with a meaningful fileName and the extension .mpl. So aver.mpl would be a reasonable choice

Any time you want to use it in any Maple worksheet, just issue the command read("aver.mpl")

The commonest way for this to fail is that Maple cannot find the file on the basis of its current path variable, so you may have to modify the read command to include either an absolute or relative path to the file, so the read command might end up something like

read( "C:/myMapleProJects/someUtilities/aver.mpl") with an absolute path.

Or if you always run your Maple projects in the directory "C:/myMapleProJects", then the relative path

read("./someUtilities/aver.mpl") would work.

 

Well, after fixing a few typos, such as deleting the line

for ode3 and ode4:

and inserting bcs2 into

eval({ode3,ode4,bcs1},w=G[k]

And greatly simplifying some of the rest of your code in the same way that I did in my answer to your previous question, the answer which you do not want seems to be the correct answer. See the attached worksheet

tp2.mw

If you really, really know that the answer is

[-0.2], [0.51553], [0.4000]

[-0.1], [0.57000], [0.4371]

[0.], [0.62756], [0.4764]

[0.1], [0.68811], [0.5176]

[0.2], [0.75153], [0.5609]

then there is an error in the problem set up somewhere - incorrect parameter values?? bcs???

 

whilst this approach might work for the trivial problem which you gave, it is doomed to failure whenever you have more than one abs() term in your expression(s).

Why do you think that

because maple may has difficulties when dealing with absolute values

What evidence do you have for this statement?

What you want to plot is a little unconveentional. The attached worksheet produces a couple of plots.

The first is the conventional f(eta) vs eta, and theta(eta) vs eta

The second is theta(eta) vs L, for eta = 0,1,2,3, cos I wasn't sure which value of eta you were interested in

trickPlot.mw

Not sure how the first part of your question relates to the second: converting your integrand to polar coordinates is trivial, but not so for the limits.

Given the quoted limits, your question only makes sense if the integration wrt x is performed first, so the solution to your integral would be given by

int( int(9*x*y^2, x=-sqrt(a^2-y^2)..0), y=0..1);

So why do you want polar coordinates???

Since your polynomials are univariate, I suggest you examine the PolynomialTools[CoefficientList] command. Use of this command ensures that polynomial coefficients appear in defined order (with an option to reverse).

My solution for your problem is shown in the attached

coeffs.mw

@Kitonum 

No real problem with the above but


for k to 9999 do
     q:=k;
     if add
        ( irem(q, 10, 'q')^3,
           k = 1..length(q)
        )
        >= k
     then n:=k
     fi;
od:
n;

runs about 15x faster

restart;
  with(Statistics):
  with(Student[Statistics]):
  with(plots):
#
# generate some data
#
  N1:=NormalRandomVariable[Normal](0,1):
  N2:=NormalRandomVariable[Normal](0.5, 0.5):
  A:=Sample(N1,100):
  B:=Sample(N2,100):
#
# Produce two histograms
#
  P1:=Histogram(A, style=polygon, color="Red"):
  P2:=Histogram(B, style=polygon, color="blue"):
#
# display them
#
  display(P1,P2);

VectorCalculus[CrossProduct](a,b), requires that a, b be 3-D vectors, yours are 2-D. Easiest thing to do is just add a third entry and set it equal to zero as in

restart;
r:=(u,v)->Vector[column]([u*cos(v), u*sin(v),0]);
ru:=diff~(r(u,v),u);
rv:=diff~(r(u,v),v);
VectorCalculus[CrossProduct](ru,rv);

Fairly easy to filter the array with

 

val:=NULL;
Z:=Array(1..5);
Z[1]:=val;
Z[2]:=5;
Z[3]:=val;
Z[4]:=val;
Z[5]:=8;
data:=[seq(`if`(Z[n]<>val,[n,Z[n]],NULL), n=1..5)];
t1:=plot(data, style=point, symbol = solidcircle, color= red):
plots[display]({t1});

Whatever is assigned to val (NULL, a numeric, it doesn't matter) will be omitted from the data plotted

In the code loop

for i from 1 to 2*r+1 do
CCC[i]:=simplify(Tr(X).H[i]-Tr(X0).H[i]-3*(Tr(X).Dtau1.CC[i].Tr(Dtau2).U));
od;

H[i] is a scalar but Tr(X) and Tr(X0) are vectors so the matrix multiplication operation '.' is invalid. Use '*' instead.

Furthermore

Tr(X)*H[i] is a 1*3 vector
Tr(X0)*H[i] is a 1*3 vector, but
-3*(Tr(X).Dtau1.CC[i].Tr(Dtau2).U)) is a scalar -and so cannot be combined with the preceding vectors. If you mean this term to be add to each element of the preceding vectors then use

CCC[i]:=simplify(Tr(X)*H[i]-Tr(X0)*H[i]-~3*(Tr(X).Dtau1.CC[i].Tr(Dtau2).U));

Couldn't check the subsequent code becuase the next execution group contains


J:=3*(Tr(X).X+Tr(U).U);
 MM:=Minimize(J,{seq(CCC[i]=0,i=1..N*M)});
 MINIMUM_IS:=op(1,MM);


Both N and M are undefined, so the seq won't execute - based on the generation of CCC in the preceding loop, 2*r+1 looks like a good bet for the last argument to seq, but I decided that this was piling up too much guesswork about what you wished to achieve

 

I find that if I use

dataplot(L0[()..(), 1], L0[()..(), 2], axes = box, size = [600, 320] ,symbol=point);

in Maple2015, then this gives the same result (at least "visually") as

plot(L0[()..(), 1], L0[()..(), 2], axes = box, size = [600, 320]);

in Maple 18 (apart from the fact that they use different colors by default - trivial to fix).

The dataplot command is new in Maple2015 but AFAIK, it is meant to be a "simpler" interface to the same underlying commands. All of which means that the behaviour of the underlying plot command (ie the ones you used) has been changed - I don't know why, so I think you have a legitimate gripe and if you don't get a better answer here than this one, then I suggest you send it to the Maple techies as a bug.

A few (probably irrelevant) comments

  1. Your data file contain roughly 10 identical x-values with different y-values: I say roughly 10 because it seems to vary from about 6 to about 11. when first examining your data I thought this might be an issue and upsetting the plot, but it doesn't *seem* to be. I tried plotting the second column of you matrix against row number: Maple 18 and Maple2015 still gave substantially different results, so I don't think this is an issue
  2. interesting way to select columns of a matrix the L0[()..(),1]. Normally I'd just use L0[..,1]

Couldn't really follow what you code was trying to do. An implemenatation of the logical process you went through in the statement of your question would be

   restart;
#
# Initialise 4 bags with nothing in any of them
# and an index for the "current bag
#
   bags:=table([b1=0, b2=0, b3=0, b4=0]):
   bi:=0:
#
# loop round all sums which have to be achieved
#
   for j from 1 to 15 do
     #
     # check whether j cannot be obtained as a sum
     # of existing bags
     #
     if member( j,
                      map
                      ( add,
                        combinat[choose] ( [entries( bags, 'nolist' ) ] )
                      )
                   )=false
     then
         #
         # j must go in the next available bag
         #
            bi:= bi+1:
            bags[ b||bi ]:= j:
     fi;
   od;
#
# display result
#
    op(bags);

First 99 100 101 102 103 104 105 Page 101 of 106