Carl Love

Carl Love

28055 Reputation

25 Badges

12 years, 353 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

Before answering your main question, I want to point out that angles need to be in radians in Maple. You can enter them in degrees and then convert to radians, like this

Angles:= [30, 60, 90];
Angles:= Angles *~ Pi/180:

or

Angles:= map(convert, Angles *~ degrees, radians):

Second, in my previous post, I said that you'd probably want to use more than 4 digits in actual practice, but you kept it at 4. You should change all occurences of evalf[4] to evalf[15] or simply evalf (for reasons that I won't go into here, 15 digits usually gives the optimal balance of speed and accuracy).

Third, you should probably use different names for the two Constructor procedures: the one that constructs the Ts and the one that constructs the Ths.

Now for your main questions. You have two lists of eight matrices each, Ts and Ths, and a single matrix Q, and you want to compute a list of matrices Qb such that Qb[k] = Ts[k]^(-1).Q.Ths[k] for k from 1 to 8. Here are two ways do that:

Qb:= [seq](Ts[k]^(-1).Q.Ths[k], k= 1..8):

or

Qb:= zip((X,Y)-> X^(-1).Q.Y, Ts, Ths):

Finally, you want to construct a 6x6 matrix from 3x3 matrices A, B, and D, with pattern

A B
B D.

The command for this is

M:= < < A | B > , < B | D > >:

Let me know if that is all clear to you!

Sorry that I didn't see this earlier. This system is very bad at notifying authors when there are responses to their answers.

So, I promised to help you again if you found two of the five errors. You found four, and correctly fixed three. The one that you incorrectly fixed is the definition of h. It should be

h:= (b-a)/2/n;   (or, equivalently, h:= (b-a)/(2*n);.)

The last error is that subtle one. It's another fencepost error. What is the sequence of coefficients of the evaluation points? If n = 2, then the coefficients are

1, 4, 2, 4, 1

and if n = 3, then the coefficients are

1, 4, 2, 4, 2, 4, 1.

Do you agree? So, there is always one more 4 than there are 2s. But your loop does the 4s and 2s in pairs so there are always the same number of them. You are missing the last 4. It occurs at x = b - h. So, your initialization of aa should be

aa:= f(a)+4*f(b-h)+f(b);

 

You should first test your procedures against Student:-Calculus1:-ApproximateInt (like I described in the previous post) before checking for convergence.

I came up with two ways to evaluate the the integral of Ei in this problem: an Ei identity and integration by parts. The identity is Ei(1,x) = -Ei(-x).

restart;

DQ:= -((D@@2)(y))(r)-2*y(r)/r+y(r) = ((4/9-exp(-r)/r)*2)*r*exp(-r):

Y:= rhs(dsolve(DQ, y(r)));

exp(-r)*r*_C2+(2*exp(-r)*Ei(1, -2*r)*r+exp(r))*_C1-(8/9)*((1/4)*(Int((2*Ei(1, -2*r)*r+exp(2*r))*(4*r*exp(r)-9)*exp(-3*r), r))*r*exp(3*r)+((r+1/2+r^2)*exp(r)-1/2-(3/2)*r)*r*Ei(1, -2*r)+((1/2)*r^2+(1/2)*r+1/4)*exp(3*r)-(3/4)*exp(2*r)*(r+1/3))*exp(-4*r)

subs(Ei(1,-2*r)= -Ei(2*r), %);

exp(-r)*r*_C2+(-2*exp(-r)*Ei(2*r)*r+exp(r))*_C1-(8/9)*((1/4)*(Int((-2*Ei(2*r)*r+exp(2*r))*(4*r*exp(r)-9)*exp(-3*r), r))*r*exp(3*r)-((r+1/2+r^2)*exp(r)-1/2-(3/2)*r)*r*Ei(2*r)+((1/2)*r^2+(1/2)*r+1/4)*exp(3*r)-(3/4)*exp(2*r)*(r+1/3))*exp(-4*r)

value(%);

exp(-r)*r*_C2+(-2*exp(-r)*Ei(2*r)*r+exp(r))*_C1-(8/9)*((1/4)*(3*exp(-r)+(1/2)*Ei(2*r)*(8*r^2*exp(-2*r)+8*r*exp(-2*r)+4*exp(-2*r)-12*r*exp(-3*r)-4*exp(-3*r))-4*r-2*ln(2*r)-2*Ei(1, r))*r*exp(3*r)-((r+1/2+r^2)*exp(r)-1/2-(3/2)*r)*r*Ei(2*r)+((1/2)*r^2+(1/2)*r+1/4)*exp(3*r)-(3/4)*exp(2*r)*(r+1/3))*exp(-4*r)

simplify(%);

exp(-r)*r*_C2-2*exp(-r)*Ei(2*r)*_C1*r+exp(r)*_C1+(4/9)*exp(-r)*ln(r)*r+(4/9)*exp(-r)*ln(2)*r+(4/9)*exp(-r)*Ei(1, r)*r+(4/9)*r^2*exp(-r)-(4/9)*exp(-r)*r-(2/9)*exp(-r)+(2/9)*exp(-2*r)

YY:= %:

Verify that it satifies the original equation.

eval((rhs-lhs)(DQ), [(D@@2)(y)(r)= diff(YY,r$2), y(r)= YY]);

2*(4/9-exp(-r)/r)*r*exp(-r)+(4/9)*exp(-r)/r+(8/9)*exp(-r)-(4/9)*(exp(-r))^2/r+2*(exp(-r)*r*_C2-2*exp(-r)*Ei(2*r)*_C1*r+exp(r)*_C1+(4/9)*exp(-r)*ln(r)*r+(4/9)*exp(-r)*ln(2)*r+(4/9)*exp(-r)*Ei(1, r)*r+(4/9)*r^2*exp(-r)-(4/9)*exp(-r)*r-(2/9)*exp(-r)+(2/9)*exp(-2*r))/r-2*exp(-r)*_C2-(8/9)*exp(-r)*ln(r)-(8/9)*exp(-r)*ln(2)-(8/9)*exp(-r)*Ei(1, r)+4*exp(-r)*Ei(2*r)*_C1-2*exp(-r)*exp(2*r)*_C1/r+(2/3)*exp(-2*r)+(4/3)*(exp(-r))^2-(16/9)*exp(-r)*r

simplify(%);

0

The identity is not valid everywhere. We can do it without the identity by using integration by parts.

 

Extract the unevaluated integral:

J:= indets(Y, Int(algebraic$2))[];

Int((2*Ei(1, -2*r)*r+exp(2*r))*(4*r*exp(r)-9)*exp(-3*r), r)

with(IntegrationTools):

simplify(Expand(expand(J)));

8*(Int(exp(-2*r)*Ei(1, -2*r)*r^2, r))-18*(Int(exp(-3*r)*Ei(1, -2*r)*r, r))+4*(Int(r, r))-9*(Int(exp(-r), r))

JJ:= value(%);

8*(int(exp(-2*r)*Ei(1, -2*r)*r^2, r))-18*(int(exp(-3*r)*Ei(1, -2*r)*r, r))+2*r^2+9*exp(-r)

jj:= Parts(JJ, Ei(1,-2*r), applytoall);

-2*(2*r^2+2*r+1)*exp(-2*r)*Ei(1, -2*r)-4*r-2*ln(r)+2*(3*r+1)*exp(-3*r)*Ei(1, -2*r)-6/exp(r)-2*Ei(1, r)+9*exp(-r)

YY:= simplify(eval(Y, J= jj));

exp(-r)*r*_C2+2*exp(-r)*Ei(1, -2*r)*_C1*r+exp(r)*_C1+(4/9)*exp(-r)*ln(r)*r+(4/9)*exp(-r)*Ei(1, r)*r+(4/9)*r^2*exp(-r)-(4/9)*exp(-r)*r-(2/9)*exp(-r)+(2/9)*exp(-2*r)

Verify that it satifies the original equation.

eval((rhs-lhs)(DQ), [(D@@2)(y)(r)= diff(YY, r$2), y(r)= YY]);

2*(4/9-exp(-r)/r)*r*exp(-r)+(4/9)*exp(-r)/r+(8/9)*exp(-r)-(4/9)*(exp(-r))^2/r-2*exp(-r)*_C2-(8/9)*exp(-r)*ln(r)-(8/9)*exp(-r)*Ei(1, r)+2*(exp(-r)*r*_C2+2*exp(-r)*Ei(1, -2*r)*_C1*r+exp(r)*_C1+(4/9)*exp(-r)*ln(r)*r+(4/9)*exp(-r)*Ei(1, r)*r+(4/9)*r^2*exp(-r)-(4/9)*exp(-r)*r-(2/9)*exp(-r)+(2/9)*exp(-2*r))/r-2*exp(-r)*exp(2*r)*_C1/r-4*exp(-r)*Ei(1, -2*r)*_C1+(2/3)*exp(-2*r)+(4/3)*(exp(-r))^2-(16/9)*exp(-r)*r

simplify(%);

0

 


Download Ei_Integration.mw

Check out command ?RootFinding,Analytic.

Here's an example. Here I use 2x2 matrices to save space. I'm sure that you'll be capable to generalize this. The essence is the use of unapply to turn one matrix with a parameter into a constructor for other matrices. (The ?unapply command is much more general than just being used for matrices.)

restart;

Make a Matrix that depends on a unspecified angle:

M:= < < cos(theta),sin(theta) > | < -sin(theta),cos(theta) > >;

M := Matrix(2, 2, {(1, 1) = cos(theta), (1, 2) = -sin(theta), (2, 1) = sin(theta), (2, 2) = cos(theta)})

Make that Matrix into a procedure that constructs more Matrices:

Constructor:= unapply(M, theta):

Make a list (not a vector) of 8 angles:

Angles:= [seq](Pi*k/4, k= 0..7);

[0, (1/4)*Pi, (1/2)*Pi, (3/4)*Pi, Pi, (5/4)*Pi, (3/2)*Pi, (7/4)*Pi]

Apply the Conctructor to the list of angles to create a list of Matrices:

Mlist:= Constructor ~ (Angles);

Mlist := [Matrix(2, 2, {(1, 1) = 1, (1, 2) = 0, (2, 1) = 0, (2, 2) = 1}), Matrix(2, 2, {(1, 1) = (1/2)*sqrt(2), (1, 2) = -(1/2)*sqrt(2), (2, 1) = (1/2)*sqrt(2), (2, 2) = (1/2)*sqrt(2)}), Matrix(2, 2, {(1, 1) = 0, (1, 2) = -1, (2, 1) = 1, (2, 2) = 0}), Matrix(2, 2, {(1, 1) = -(1/2)*sqrt(2), (1, 2) = -(1/2)*sqrt(2), (2, 1) = (1/2)*sqrt(2), (2, 2) = -(1/2)*sqrt(2)}), Matrix(2, 2, {(1, 1) = -1, (1, 2) = 0, (2, 1) = 0, (2, 2) = -1}), Matrix(2, 2, {(1, 1) = -(1/2)*sqrt(2), (1, 2) = (1/2)*sqrt(2), (2, 1) = -(1/2)*sqrt(2), (2, 2) = -(1/2)*sqrt(2)}), Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = -1, (2, 2) = 0}), Matrix(2, 2, {(1, 1) = (1/2)*sqrt(2), (1, 2) = (1/2)*sqrt(2), (2, 1) = -(1/2)*sqrt(2), (2, 2) = (1/2)*sqrt(2)})]

(Optional) Evaluate to floating point (in actual usage, you'd probably want more than 4 digits).

Mlist:= evalf[4](Mlist);

Mlist := [Matrix(2, 2, {(1, 1) = 1., (1, 2) = 0., (2, 1) = 0., (2, 2) = 1.}), Matrix(2, 2, {(1, 1) = .7070, (1, 2) = -.7070, (2, 1) = .7070, (2, 2) = .7070}), Matrix(2, 2, {(1, 1) = 0., (1, 2) = -1., (2, 1) = 1., (2, 2) = 0.}), Matrix(2, 2, {(1, 1) = -.7070, (1, 2) = -.7070, (2, 1) = .7070, (2, 2) = -.7070}), Matrix(2, 2, {(1, 1) = -1., (1, 2) = 0., (2, 1) = 0., (2, 2) = -1.}), Matrix(2, 2, {(1, 1) = -.7070, (1, 2) = .7070, (2, 1) = -.7070, (2, 2) = -.7070}), Matrix(2, 2, {(1, 1) = 0., (1, 2) = 1., (2, 1) = -1., (2, 2) = 0.}), Matrix(2, 2, {(1, 1) = .7070, (1, 2) = .7070, (2, 1) = -.7070, (2, 2) = .7070})]

Individual Matrices in the list are specified by index, 1 to 8.

Mlist[4];

Matrix(2, 2, {(1, 1) = -.7070, (1, 2) = -.7070, (2, 1) = .7070, (2, 2) = -.7070})

 

Download Matrix_function.mw

Yes, the "problem" is caused by the lack of end proc. You can avoid the problem two ways:

  1. When you want a blank line, use Shift-Enter (or Shift-Return) instead of Enter (or Return). Don't hit Enter until you're ready for the procedure to be processed ("parsed" is the lingo). If you do hit Enter prematurely, just ignore any error message and put the cursor back where you want it. This takes a bit of getting used to, which brings us to my preferred method...
  2. From the Insert Menu, create a Code Edit Region. This gives you a box, a subwindow, that you can use like Microsoft Notepad or other primitive text editor. Here, Enter just gives you a new line. When you want the code to be parsed (and Executed, if it's more than just a procedure definition), right-click in the box and select Execute Code. You can resize the box by left-clicking in it and then dragging one of the blue square dots on the sides or corners. But the box doesn't need to be big enough to hold your code; it will scroll if needed.

Put a semicolon after the end proc.

The proc line should not end with a semicolon.

There is also a mistake in your type declaration, but the mistake will unfortunately not generate an error message. If you want a parameter to be restricted to having both type A and type B, then the type is And(A,B). So with that the proc line becomes

Grade:= proc(mark::And(nonnegative,numeric))

You should look up a user named Rupunzel on this forum. She got help about this very same homework problem in two separate threads a few days ago. There's a Users tab at the top of the page.

For me, in worksheet text mode, the alt codes generate different characters than the characters accessible programmatically by Maple. (I haven't explored the possibility of getting more characters with MathML or XML.) But, I see that you got a thorn for character 254, which is the same as Maple's standard font's 254. So, I can just plot the standard font for you:

plots:-display(
    plots:-pointplot([[0,0],[16,16]], color= white),
    plots:-textplot(
        [seq([.5+irem(k,16), .5+iquo(k,16), convert([k],bytes)
             ,color= black, font= [times, roman, 24]
             ]
             ,k= 0..255
             )
        ,seq([irem(k,16), .75+iquo(k,16), sprintf("%d",k)
             ,color= blue, font= [times, roman, 9], align= {above, right}
             ], k= 0..255
             )
        ]
    ),
    labels= ["low nybble", "high nybble"],
    labeldirections= [horizontal,vertical],
    axes= boxed,
    axis= [gridlines= [$(1..15)]],
    tickmarks= [(1/2 +~ [seq](k, k= 0..15) =~ [$(0..9),A,B,C,D,E,F]) $ 2],
    axesfont= [times, italic, 16],
    labelfont= [HELVETICA, bolditalic, 16],
    title= "Character codes 0-255 in the regular font",
    titlefont= [HELVETICA, BOLDOBLIQUE, 24]
);

When you plot this in Maple, you'll need to enlarge the plot with your mouse so that the characters fit neatly in their boxes.

As posted, you have overlapping ranges for x, so I think that you meant for the upper bound of x in the first line to start -2*Pi instead of 2*Pi, and my code below reflects that change. Note that the constant Pi is capitalized in Maple. That being said, what I get from your ranges is a nonagon, not a hexagon.

Plotting a surface over a complicated region in the plane involves the same steps as setting up a double integral over that region. You've already done the hardest part---defining the ranges. I chose an f, and here's the rest of the steps.

rt3:= sqrt(3):  C:= Pi/3/rt3:

R:= [[x= -4*C..-2*C, y= -rt3*(x-4*C)..rt3*(x-4*C)]
    ,[x= -2*C..2*C, y= -2*Pi/3..2*Pi/3]
    ,[x= 2*C..4*C, y= -rt3*(4*C-x)..rt3*(4*C-x)]
    ]:

f:= (x,y)-> x^2+y^2:

plots:-display((r-> plot3d(f(x,y), r[])) ~ (R));

 

Download poly.mw

 

After executing the statement

gn:= g;

gn is identical to g, so any changes made to gn will automatically be made to g also. To avoid this, i.e. to create a separate independent matrix, use copy:

gn:= copy(g);

But I can't tell for sure if this is what you need because you never use gn after constructing it. So the only way that the code as posted makes sense is if the changes are made to g. Perhaps you didn't upload the whole code.

Also, the double if statement at the end can be replaced by

if not (j in s[i] xor n in s[i]) then g[j,n]:= 0 end if;

To avoid unintended conversion of monomials, you can check if the expression is of type `+`.

Example:

p:= randpoly([x,y,z]);

-62*x^2*z^3+97*x*y^3*z-73*y*z^4-56*x*y*z^2+87*x*y

evalb(p::`+`);

true

`if`(p::`+`, nops, 1)(p);

5

L:= `if`(p::`+`, [op], `[]`)(p);

[-62*x^2*z^3, 97*x*y^3*z, -73*y*z^4, -56*x*y*z^2, 87*x*y]

tL:= indets ~ (L);

[{x, z}, {x, y, z}, {y, z}, {x, y, z}, {x, y}]

nL:= nops ~ (tL);

[2, 3, 2, 3, 2]

[seq]({seq}(v=degree(L[k],v), v= tL[k]), k= 1..nops(L));

[{x = 2, z = 3}, {x = 1, y = 3, z = 1}, {y = 1, z = 4}, {x = 1, y = 1, z = 2}, {x = 1, y = 1}]

 

Download poly.mw

You simply need to do

eval(x(t), dsn1(1));

Let me know how that goes.

 

Building upon Markiyan's excellent observation about erf (which I wish I could vote up, but can't because it's a comment), we see that the Asker's original problem can be done thus (which I call going one step backwards before going forward):


restart;

f:= exp(-(x+sigma)^2/2/sigma^2);

exp(-(1/2)*(x+sigma)^2/sigma^2)

diff(int(f,x), x$(n+1));

(1/2)*2^(1/2)*sigma*((1/2)*2^(1/2)*x/sigma+(1/2)*2^(1/2))^(-n)*MeijerG([[0, 1/2], []], [[0], [(1/2)*n+1/2, (1/2)*n]], ((1/2)*2^(1/2)*x/sigma+(1/2)*2^(1/2))^2)*2^(n+1)*((1/2)*2^(1/2)/sigma)^(n+1)

And also we see that Maple's algorithm for arbitrary order differentiation could be very easily improved. It would be trivial for it to check this method after other methods have failed.

Download errdiff.mw

You've made the classic so-called fencepost error. Think of the trapezoids as the gaps between posts of a fence, and think of the evaluation points as the fenceposts. There are n of gaps, hence n+1 posts. In your CompTrap, you have n+2 evaluation points: a, b, and n done in the loop. The loop needs to stop at n-1.

Also, in your test comparison integral, you are integrating from 8 to 1. You should go from 1 to 8.

You can compare your procedure with Maple's own trapezoid rule procedure like this: The call

CompTrap(f, a, b, n);

should produce the same result (*footnote) as

Student:-Calculus1:-ApproximateInt(f(x), x= a..b, method= trapezoid, partition= n);

This is a good way to test whether your procedure is correctly implementing the rule. 

It is also good to compare against the integral done by Maple's more sophisticated methods, as you did. This is a way of testing whether the rule itself is good, not testing your implementation of it. You can also study the rate of convergence of the rule.

A matter of style in your procedure: Be careful of how you mix "pure" exact symbolic computation with floating-point numeric computation. Your procedure has 0.5 near the end. You should probably change that to exact (1/2).

Regarding your Simpson's Rule procedure: In addition to the fencepost error, there five other errors in that, three algorithmic and two syntactic. Two of the algorithmic errors are super obvious and the third is subtle; both syntactic errors are essentially the same. If you find any two of the five, I'll point out another two. Deal? Also, why did you change the parameter n from the trapezoid procedure to j in the Simpson's procedure? Maple doesn't care which you use, but the change may have clouded your thinking and led to one of the errors.

You can test/compare your procedure with Maple's own implementation of Simpson's rule by using the same Student command as above, but changing trapezoid to simpson.

(*footnote): After you apply evalf, there may be a small deviation due to rounding error. This may affect the last digit.

You were close. That's actually fairly impressive for one's first day with Maple (and the day is young in my time zone).

An unevaluated function does need to have at least a name, something to put immediately in front of the (a, b, c, d). You are right that no explict expression or -> is needed. I think from your example that you want to name it B. Do this:

restart;
mtaylor(B(a,b,c,d), [a=e,b=f,c=g,d=h], 2);

Does that give you what you want? I made the order 2 because 1 is too trivial an example for you to see what Maple is doing. Please feel free to followup here on MaplePrimes. Have fun exploring Maple.

First 384 385 386 387 388 389 390 Last Page 386 of 395