Rouben Rostamian

MaplePrimes Activity


These are replies submitted by Rouben Rostamian

@ijuptilk Here is the implementation of that modified algorithm in Maple.  Its results agree with what you have shown.

restart;

with(LinearAlgebra):

doit := proc(A::Matrix)
        local Q, R, s, id;
        s := A[-1,-1];
        id := IdentityMatrix(Dimension(A));
        Q, R := QRDecomposition(A-s*id, fullspan);
        return R.Q + s*id;
end proc:

A[0] := Matrix(  [[3.0, 1.0, 4.0], [1.0, 2.0, 2.0], [0., 13.0, 2]] );

Matrix(3, 3, {(1, 1) = 3.0000000000, (1, 2) = 1.0000000000, (1, 3) = 4.0000000000, (2, 1) = 1.0000000000, (2, 2) = 2.0000000000, (2, 3) = 2.0000000000, (3, 1) = 0., (3, 2) = 13.0000000000, (3, 3) = 2})

for i from 1 to 6 do
        A[i] := doit(A[i-1]);
end do;

Matrix(3, 3, {(1, 1) = 3.5000000000, (1, 2) = -4.2635347558, (1, 3) = .26883338042, (2, 1) = -9.2059763198, (2, 2) = 1.5766961652, (2, 3) = 9.1965598768, (3, 1) = 0., (3, 2) = -1.4100418410, (3, 3) = 1.9233038348})

Matrix(3, 3, {(1, 1) = 3.8726710918, (1, 2) = 5.8456590691, (1, 3) = 11.4347283508, (2, 1) = 4.4236981716, (2, 2) = 1.7461077570, (2, 3) = -1.8760931255, (3, 1) = 0., (3, 2) = .17939094164, (3, 3) = 1.3812211512})

Matrix(3, 3, {(1, 1) = 6.6491606142, (1, 2) = -3.0031275241, (1, 3) = -3.8697004604, (2, 1) = -4.2847775607, (2, 2) = -.63360994982, (2, 3) = 10.9574256085, (3, 1) = 0., (3, 2) = -0.14483593242e-1, (3, 3) = .98444933561})

Matrix(3, 3, {(1, 1) = 7.5051964360, (1, 2) = 3.1077455977, (1, 3) = 9.7110573795, (2, 1) = 1.8714215464, (2, 2) = -1.4597440348, (2, 3) = -6.4160015539, (3, 1) = 0., (3, 2) = 0.1396082e-3, (3, 3) = .95454759878})

Matrix(3, 3, {(1, 1) = 8.1438626575, (1, 2) = .36380697772, (1, 3) = -7.5750590186, (2, 1) = -.87218400156, (2, 2) = -2.0980217062, (2, 3) = 8.8369002971, (3, 1) = 0., (3, 2) = -0.1708443559e-7, (3, 3) = .95415904870})

Matrix(%id = 36893628787730574380)

Approximate eigenvalues

sort(Diagonal(A[i-1]));

Vector(3, {(1) = -2.0102632935, (2) = .95415900373, (3) = 8.0561042898})

Precise eigenvalues

Eigenvalues(A[0]);

Vector(3, {(1) = -2.0669460973+0.*I, (2) = .95415900373+0.*I, (3) = 8.1127870936+0.*I})

I assume that by "deflate" you mean something like this:

[seq](convert(Row(A[6],i), list), i=1..3);

[[HFloat(8.056104289776464), HFloat(1.5956054456867756), HFloat(8.584132685780313)], [HFloat(0.35961451551635476), HFloat(-2.0102632935089915), HFloat(-7.860343023905163)], [HFloat(0.0), HFloat(2.5728710926568884e-16), HFloat(0.9541590037325286)]]

Download mw2.mw

Running out of memory can cause those symptoms.  Perhaps you are running several other memory-consuming applications at the same time?  Try closing your browser and other memory hogs.

@Tokoro That's a good question.  I hadn't thought about it.

My calculations show that if we split the disk into three sectors of arbitrary sizes, then the maximum total volume is achieved when one of the sectors is of zero angle, and thus the problem reduces to the two cones case.

On the other hand, if we split the disk into three sectors so that two of the sectors are of the same size, then the maximum total volume is achieved when one of the sectors is of zero angle, and the other two are exactly half of the disk each.

@FredK You may have any number of equations and functions within a single proc.  It's difficult to provide a more concrete answer in the absence of more detailed information.

The expression
    - x(t) - 1/6*x(t)^3 + 1/120*x(t)^5
looks very much like the  first few terms of the series expansion of the sine function but the signs are wrong.  I am guessing that it is meant to be
    - x(t) + 1/6*x(t)^3 - 1/120*x(t)^5

Same goes for the y(t).

How would you solve this problem with paper and pencil?  If you present an outline that shows that you understand the math, then I will be happy to show you how to code that math into Maple.

 

There are very many ways of doing that, some more sophisticated that others, and it's not all about knowing how to use Maple; tt also depends on what sort of mathematics you know.

Can you describe how you would construct the triangle if you were doing this a hundred years ago, before there were computers? Your (detailed) answer to that question will help us pick and present to you the solution that.best matches your background.

@mehdibgh In my original response to your question, I showed how to express the curve in the form x = my_x(y), where my_x(y) is calculated numerically.

Let's say you want to calculate the area bounded between your curve and the line x=0.5. We do:

# find limits of integration
a := fsolve(my_x(y) = 0.5, y = 0 .. 0.2);
b := fsolve(my_x(y) = 0.5, y = 1.6 .. 2.0);
area := int(0.5  - my_x(y), y = a .. b, numeric);

The answer is area = 0.6770481296.  The centroid can be calculated in the same way.

@mehdibgh Integration can be done easily in the case of the specific functions that you have provided.

restart;

kernelopts(version);

`Maple 2021.2, X86 64 LINUX, Nov 23 2021, Build ID 1576349`

eq := 1/2 - sqrt(x*exp(-x^2-y^2)) = (x-1)*exp(-x^2-y^2);

1/2-(x*exp(-x^2-y^2))^(1/2) = (x-1)*exp(-x^2-y^2)

plots:-implicitplot(eq, x=0..2, y=-2..2);

Explicit solution: y as a function x

sol := solve(eq, y);

(-x^2-ln((1/2)*(2*x-1+(3*x^2-2*x)^(1/2))/(x^2-2*x+1)))^(1/2), -(-x^2-ln((1/2)*(2*x-1+(3*x^2-2*x)^(1/2))/(x^2-2*x+1)))^(1/2), (-x^2-ln(-(1/2)*(-2*x+1+(3*x^2-2*x)^(1/2))/(x^2-2*x+1)))^(1/2), -(-x^2-ln(-(1/2)*(-2*x+1+(3*x^2-2*x)^(1/2))/(x^2-2*x+1)))^(1/2)

plot([sol[3], sol[4]], x=0.7..1.5, color=[red,blue]);

Find the x domain:

a,b := fsolve(sol[3], x=0.5..1.5, maxsols=2);

.7350072152, 1.400869250

Calculate the area of the top half, then double the result

area := 2*int(sol[3], x=a..b);

.6557733696

 

 

 

Download mw.mw

 

@Navid555 I wrote "if u is zero, then its derivative is automatically zero, so your second condition is redundant and should be removed".  You seem to disagree.

I have nothing more to say.

 

 

I am sure that you know that if a function is identically zero on an interval, then its derivative is also identically zero on that interval.

Your initial conditions u(x, 0) = 0, D[1](u)(x, 0) = 0 say that u and its derivative are zero on the interval (0,4). But if u is zero, then its derivative is automatically zero, so your second condition is redundant and should be removed.

After removing D[1](u)(x, 0) = 0, your initial boundary value problem will have incomplete initial conditions. That's why Maple is unable to produce a solution.

Go back to the application that gives rise to the PDE and IBC and see what the missing condition is.

@acer That's an excellent approach, and I particularly like it since it is set up to be quite straightforward and transparent, though nontrivial.  I am going to be busy at work all day, but I will take a closer look tomorrow.  I am thinking of adding a to_sine::truefalse option which gives preference to expressing the result in terms of sines rather than cosines whenever possible.  For instance, 1 - cos(x)^2  will become sin(x)^2.

@acer That's an incredibly complex code which does incredibly complex things!  It does the expected:

> 2*sin(x)*cos(x);
                          2 sin(x) cos(x)
> H(%);
                               sin(2 x)

and also the unexpected:

 > 4*sin(x)^2*cos(x);
                                       2
                               4 sin(x)  cos(x)

> H(%);
                               cos(x) - cos(3 x)

I would have expected 2*sin(2*x)*sin(x) as the answer but the combine step in the definition of H intervenes and, for better or worse, takes it one step further.

@Carl Love Your combine/half_angle code is very clever and it rightly belongs to Maple's default collection of convert(...) functions.  Until then, I will keep it as an add-on in my initialization file.

As to the double_angle conversion, applying combine is not ideal. For instance

expr := 4*sin(x)*cos(x)*sin(y)*cos(y);
                     expr := 4 sin(x) cos(x) sin(y) cos(y)
combine(expr);
                   1/2 cos(-2 y + 2 x) - 1/2 cos(2 y + 2 x)

I would have preferred the answer of sin(2*x)*sin(2*y). I don't know how to construct an extension to convert that will do that.

 

@Carl Love Thank you very much for that example. Using subsindets is new to me and I found your example very instructive.

Here is a generalized version of your construction. I intend to place in my Maple initialization file.

restart;

`convert/half_angle` := proc(expr)
        subsindets(
                expr,
                {sin, cos, tan, cot}(anything),
                f ->
                        local x := op(f)/2;
                        if op(0,f) = sin then
                                2*sin(x)*cos(x)
                        elif op(0,f) = cos then
                                2*cos(x)^2 - 1
                        elif op(0,f) = tan then
                                2*tan(x)/(1-tan(x)^2)
                        else # cot
                                (cot(x)^2-1)/(2*cot(x))
                        end if
        );
end proc:

 

Example uses: 

convert(sin(2*a)*cos(4*b), half_angle);

2*sin(a)*cos(a)*(2*cos(2*b)^2-1)

convert(sin(8*a+4*b), half_angle);

2*sin(4*a+2*b)*cos(4*a+2*b)

convert(tan(x), half_angle);

2*tan((1/2)*x)/(1-tan((1/2)*x)^2)

convert(sin(cos(2*x)), half_angle);

2*sin(cos(x)^2-1/2)*cos(cos(x)^2-1/2)

convert(<sin(x), cos(x)>, half_angle);

Vector(2, {(1) = 2*sin((1/2)*x)*cos((1/2)*x), (2) = 2*cos((1/2)*x)^2-1})
 

Download convert-to-half-angle.mw

I also played around with the idea of making something similar for the double-angle identities:

sin(x)*cos(x) = sin(2*x) / 2;
cos(x)^2 = (1 + cos(2*x)) / 2;
sin(x)^2 = (1 - cos(2*x)) / 2;

but did not get very far.  If you have suggestions about how to go about this, I will be eager to hear them.

 

First 22 23 24 25 26 27 28 Last Page 24 of 99