Carl Love

Carl Love

28065 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

I assume that you intended for g to be considered a function of x, y, z. In that case, you should change Diff(g, x) to D(g)(x) or Diff(g(x,y,z), x), and likewise for y and z.

When you simplify the expressions in these matrices, they become incredibly long---essentially, they are expanded as polynomials in the ~ 10-20 variables. The only change I made was to remove the line kappa:= simplify(kappa). With this change, I was able to eval the K=8 case (the 254x20) in 0.25 seconds, and I did the K=10 (1022 x 24) case in 2.5 seconds and the K=11 (2046 x 26) in 8.5 seconds. The mykappa itself also runs much, much faster without the simplify. So, there, I gave you a factor of about 1000 improvement.

Not simplifying should also give you greater accuracy in the floating-point evaluation---fewer operations means less round-off error.

With some research in Wikipedia, I found an algorithm for counting the Eulerian circuits (not paths) in directed graphs (digraphs). It's called the BEST algorithm, BEST being the initials of the four authors' surnames. A digraph has an Eulerian circuit iff for every vertex v, indegree(v) = outdegree(v). When these are equal, we will just call it deg(v). Denote by N(G) the number of spanning trees of G. Then the number of Eulerian circuits is

N(G)*mul((deg(v)-1)!, v in G)

(For the graph in question, clearly all the factorials are 1.) Maple has a command for counting the spanning trees: GraphTheory:-NumberOfSpanningTrees. Applying the BEST algorithm in this case, we get 

restart:
macro(GT= GraphTheory):
G1:= GT:-Graph( {[C,G], [G,C], [T,G], [G,T], [A,T], [C,A], [T,C]}):
GT:-NumberOfSpanningTrees(G1)*mul((v/2-1)!, v in GT:-DegreeSequence(G1));
                               3

So, there are three Eulerian circuits.

(The division by 2 in the factorial is because DegreeSequence adds together the InDegree and OutDegree.)

There is no bug. Option inplace, the default, means that LinearSolve does not return a value directly; instead, it alters the input matrix.

(Dammit, the worksheet uploader is broken right now! I have to upload this worksheet in plaintext without the output. You'll have to download the worksheet to see the output.)

restart:
macro(LA= LinearAlgebra, LAM= LinearAlgebra:-Modular):
# Generate a random example linear system AX=B. The purpose of these while loops is to make sure that
# the example is not trivial.
det:= 0:
while det = 0 do
     A:= LAM:-Create(2, 3, 3, random, integer[]);
     det:= LAM:-Determinant(2, A)
end do:
N:= 0:
while N=0 do
     B:= LAM:-Create(2, 3, 1, random, integer[]);
     N:= LA:-Norm(B)
end do:

Aug:= < A | B >;
              
# By default, Modular:-LinearSolve works inplace, which means that it does not return a value directly.
# Instead, it alters the input matrix.
LAM:-LinearSolve(2, Aug, 1);
Aug;
             
X:= Aug[..,-1];
          
# Verify solution.
LAM:-Multiply(2, A, X);
          
B;
              
# Contrast with solution obtained with inplace= false.
X:= LAM:-LinearSolve(2, < A | B >, 1, inplace= false);


#Download the worksheet:  LA_Modular.mw

This works: I explicitly construct the relevant polygon from the data provided by arc:

M:= op(1, plottools:-arc([0,0], 1, Pi/6..3*Pi/4)):
plots:-polygonplot(< < M[1,1] | 0 >, M, < M[-1,1] | 0 > >);

Let me know how that goes for you.

Regarding the help: Some of the options listed on the page ?plot,options apply only to the plot command itself. I think that those options are adaptive, discont, filled, resolution, sample, and smartview. And option filledregions applies only to the three plots commands contourplot, implicitplot, and listcontplot.

Assuming that the B's are nonnegative (if the B's can be negative, you'll need to explain the significance of that):

A:= [3,5,7,2,11,2,1,10,6,15]:
B:= [0.23, 0.11, 0.05, 0.05, 0.04, 0.06, 0.10, 0.20, 0.06, 0.10]:
Statistics:-Histogram(zip(`$`, A, round~(1000./`+`(B[]) *~ B)), discrete);

How about this?

restart:
with(GraphTheory):
G1:= Graph(
     [A, AT, T, CA, TC, GT, TG, C, CG, G, GC],
     {
          [C,CG], [CG,G], [G,GC],
          [GC,C], [G,GT], [GT,T],
          [T,TG], [TG,G], [A,AT], [AT,T],
          [C,CA], [CA,A], [T,TC], [TC,C]
     }
):

SetVertexPositions(
     G1,
     [
          [0,2], [1.5,2], [3,2],
          [0,1], [1.5,1], [3,1], [3.5,1],
          [0,0], [1.5,0], [3,0],
                 [1.5,-.5]
     ]
);
HighlightVertex(G1, [A,T,C,G], red);
HighlightVertex(G1, [AT,CA,TC,GT,TG,CG,GC], grey);

DrawGraph(G1);

I am guessing that by "labelled" you mean that you want the contours themselves labelled. That would perhaps be possible with some carefully placed textplots. But perhaps the following is good enough: I've placed the labels on the z-axis and some eye-tracking lines from the labels to the contours.

Contours:= [10*k $ k =-3..2]:
P1:= plots:-contourplot3d(
     x^3-y^2, x= -3..3, y= -3..2,
     axes= box, shading= zhue, filled,
     contours= Contours,
     tickmarks= [default$2, Contours],
     transparency= 0.5, thickness= 3
):
P2:= seq(
     plots:-pointplot3d(
          [[-3,2,c], [3,2,c]],
          style= line, color= gray
     ),
     c= Contours
):
plots:-display([P2,P1]);

The following command will create the graph. Having Maple draw the graph exactly the same way as you did is another matter. But Maple's graph is isomorphic to yours.
restart:
with(GraphTheory):
G1:= Graph( {[C,G], [G,C], [T,G], [G,T], [A,T], [C,A], [T,C]});
   
DrawGraph(G1);

graph.mw

The trick is to introduce a new independent variable for the derivative and to link it to the other independent variables. But pdsolve(..., numeric) is very fussy about how you introduce it. I tried three ways that didn't work before I got the working way below.

(**)

restart:

Abbreviations:

(**)

U:= u(y,t):  T:= theta(y,t):  UY:= Uy(y,t):

I've substituted Uy(y,t) for diff(u(y,t), y) in the first PDE.

(**)

pde1:= diff(Uy(y,t),y) = S*Uy(y,t) + diff(U,t) + M*U + U/K - T;

diff(Uy(y, t), y) = S*Uy(y, t)+diff(u(y, t), t)+M*u(y, t)+u(y, t)/K-theta(y, t)

No change in PDE 2.

(**)

pde2:= diff(T,y$2) = Pr*(diff(T,t) + S*diff(T,y));

diff(diff(theta(y, t), y), y) = Pr*(diff(theta(y, t), t)+S*(diff(theta(y, t), y)))

I added a PDE to link the new dependent variable to the others.

(**)

pde3:= diff(U,y) = Uy(y,t);

diff(u(y, t), y) = Uy(y, t)

(**)

IBC:= {
     u(y,0)=0, theta(y,0)=0, u(0,t)=0,
     theta(0,t)=t0, u(1,t)=0, theta(1,t)=0
};

{theta(0, t) = t0, theta(1, t) = 0, theta(y, 0) = 0, u(0, t) = 0, u(1, t) = 0, u(y, 0) = 0}

Parameters (I changed your parameter t to t0 because t is already used as an independent variable):

(**)

Pr:= 7:  S:= 0.6:  M:= 1:  K:= 0.5:  t0:= 0.2:

(**)

Sol:= pdsolve({pde1,pde2,pde3}, IBC, numeric);

module () local INFO; export plot, plot3d, animate, value, settings; option `Copyright (c) 2001 by Waterloo Maple Inc. All rights reserved.`; end module

(**)

dudy:= Sol:-value(Uy(y,t)):

(**)

dudy(0,3);

[y = 0., t = 3., Uy(y, t) = 0.682213955850606429e-1]

(**)

Sol:-plot(Uy(y,t), y= 0, t= 0..5);

(**)

dudy0:= tau-> eval(Uy(y,t), dudy(0,tau));

proc (tau) options operator, arrow; eval(Uy(y, t), dudy(0, tau)) end proc

(**)

plot(dudy0, 0..5);

(**)

dudy0(3.);

0.682213955850606429e-1

(**)

 

 

 

Download dudy0.mw

For your first question, enter the command

Digits:= 15;

or whatever number of digits you want.

 

For your second question, make the loop like this:

for i from 1 to 5 do
     A := i;
     B := 2*i;
     C := A+B;
     print('C' = C)
end do:

Whether you use colons or semicolons inside the loop makes no difference; what matters is what you use on the end do.

Here's how to do simplify with side relations (modulo what Markiyan said), plus a few other techniques to simplify both the input and the output of this expression.

 

(**)

mul((1+x[k]+y[k])*(1-x[k])*(1-y[k]), k= 1..4);

(1+x[1]+y[1])*(1-x[1])*(1-y[1])*(1+x[2]+y[2])*(1-x[2])*(1-y[2])*(1+x[3]+y[3])*(1-x[3])*(1-y[3])*(1+x[4]+y[4])*(1-x[4])*(1-y[4])

(**)

simplify(%, {seq}([x[k]^2, y[k]^2][], k= 1..4));

x[1]*x[2]*x[3]*x[4]*y[1]*y[2]*y[3]*y[4]-x[1]*x[2]*x[3]*y[1]*y[2]*y[3]-x[1]*x[2]*x[4]*y[1]*y[2]*y[4]-x[1]*x[3]*x[4]*y[1]*y[3]*y[4]-x[2]*x[3]*x[4]*y[2]*y[3]*y[4]+x[1]*x[2]*y[1]*y[2]+x[1]*x[3]*y[1]*y[3]+x[1]*x[4]*y[1]*y[4]+x[2]*x[3]*y[2]*y[3]+x[2]*x[4]*y[2]*y[4]+x[3]*x[4]*y[3]*y[4]-x[1]*y[1]-x[2]*y[2]-x[3]*y[3]-x[4]*y[4]+1

(**)

factor(%);

(x[4]*y[4]-1)*(x[3]*y[3]-1)*(x[2]*y[2]-1)*(x[1]*y[1]-1)

(**)

 

 

Download siderels.mw

When solving for the real part of y, the y2 cancels miraculously, leaving something that can be plotted with plot3d.

 

(**)

E:= y^2+x^2=1;

x^2+y^2 = 1

(**)

E1:= eval(E, [y= y1+I*y2, x= x1+I*x2]);

(x1+I*x2)^2+(y1+I*y2)^2 = 1

(**)

E2:= solve(E1, y1);

-I*y2+(-(2*I)*x1*x2-x1^2+x2^2+1)^(1/2), -I*y2-(-(2*I)*x1*x2-x1^2+x2^2+1)^(1/2)

(**)

E3a:= evalc(Re(E2[1])); E3b:= evalc(Re(E2[2]));

(1/2)*(2*((-x1^2+x2^2+1)^2+4*x1^2*x2^2)^(1/2)-2*x1^2+2*x2^2+2)^(1/2)

-(1/2)*(2*((-x1^2+x2^2+1)^2+4*x1^2*x2^2)^(1/2)-2*x1^2+2*x2^2+2)^(1/2)

(**)

P1:= plot3d(E3a, x1= -1.5..1.5, x2= -1.5..1.5):

(**)

P2:= plot3d(E3b, x1= -1.5..1.5, x2= -1.5..1.5):

(**)

plots[display]([P1,P2], scaling= constrained, title= "Real part");

(**)

 

 

Download real_part.mw

By default, Maple does not check for discontinuities in a plot. So it picks x-values too close to the vertical asymptotes, which produces very large y-values. You can get Maple to check for discontinuities by including the discont option to plot.

plot([(2*x+3)/(x-4), (4*x+3)/(x-2)], discont);

First 357 358 359 360 361 362 363 Last Page 359 of 395