Rouben Rostamian

MaplePrimes Activity


These are replies submitted by Rouben Rostamian

@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.

 

@acer Thanks for the explanation.  I see the fault of my original approach and I sorry for having called it a "bug".

@acer Thanks for your detailed explanation. It's clear to me now why my approach wasn't working. It's the map[2] command that I needed. Thanks again.

@Kitonum The others have pointed out what I was doing wrong and how to fix it.  That said, I very much like your approach because that we we may do:

subsindets(sin(a+b), sin(anything),t->2*sin(op(1,t)/2)*cos(op(1,t)/2));

and obtain 2 sin(a/2 + b/2) cos(a/2 + b/2)..

That brings up a question:  Is is possible to apply two transformation rules in a single call to subsindets?  Specifically, can we call subsindents just once in order to apply the two rules defined in half_angle_rule in my original post to the expression sin(x)+cos(x)?

@Carl Love Thanks for observation. Now I see what I was doing wrong.

@jetboo Your messy notation may be hiding a simple structure.  Perhaps you can

work on simplifying the notation. For one thing, your equations don't seem

to be coupled, so why are you attempting to solve them as a system?

 

For now, let's look at your fourth equation, that is

diff(vhat(vars), f(x)) = a_4*v + u*diff(f(x), x).

In the interest of a simpler notation, let us write this in mathematically

equivalent but more understandable form:
v'(f(x)) = g(x).

where v' is the derivative of v, and v'(f(x)) means the derivative of v

evaluated at f(x).  The right-hand side, g(x), is any known function of x.
It makes no difference if it has f'(x) in it or not.

You need to decide whether you want to calculate v(f(x)) or v(x).

 

1. Let's say you want to calculate v(f(x)).  Call it w(x), that is, let

w(x) = v(f(x)).  Then by the chain rule we have

w'(x) = v'(f(x)) f'(x) = f'(x) g(x). Therefore w(x) is the integral

of f'(x) g(x) and we are done.

 

2. Let's say you want to calculate v(x).  Let f(x) = q.

Assume f is invertible. Then x = "f^(-1)"(q). Then your equation takes
the form v'(q) = g("f^(-1)"(q))  or, with a change of notation, v'(x) = g("f^(-1)"(x)).
Therefore, v(x) is the integral of g("f^(-1)"(x)), and we are done.

 

Both examples run fine on Linux.  I tried them in GUI and also in a character-based terminal.

Since Mac is really a unix clone under the hood, your error may stem from the stacksize limit in your unix shell.  I have no experience with Macs, so I don't know what sort of a shell is its defaul.  If it is something like the classical csh or tcsh, type "limit" at a terminal to see what your stacksize is. Mine says "stacksize 8192 kbytes".  You may raise the stacksize if needed.  If your shell is something like sh or bash, the command is "ulimit".

@ogunmiloro You have a linear system of 11 differential equations in 11 unknowns. For the equilibrium to be stable, all eigenvalues of the coefficient matrix should be on the left-hand side of the complex plane. The Routh-Hurwitz criterion (look it up in Wikipedia) provides necessary and sufficient conditions for that. 

In the attached worksheet I have extended my previous calculations to analyze the locations of the eigenvalues through the Routh-Hurwitz criterion. In the end, we see that the stability of the system is determined by two easy-to-verify explicit inequalities that should hold among the system's parameters. We also see that the inequalities may or may not hold, depending of the relative magnitudes of those  parameters.

Have a look at those inequalities and see whether you have enough information about your model to justify their validity.  In the worst case, you will need to supply numerical values for the parameters to verify that the inequalities hold.

This site refuses me to display the contents of my worksheet.  You will need to download it and view it in Maple.

Download worksheet: mw2.mw

First 14 15 16 17 18 19 20 Last Page 16 of 91