Carl Love

Carl Love

26344 Reputation

25 Badges

11 years, 134 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

The transformation that you propose has already been implemented in the reverse direction, the right side of your equation becoming the left side. In other words, the following happens automatically:

diff(inttrans:-fourier(U(x,t), x, s), t);  lprint(%);
fourier(diff(U(x,t),t),x,s)

I'm not writing that to say that what you want can't be done; I'm just saying that we need to take it into account when implementing your proposal. If it's not taken into account, your transformed result (i.e., diff@fourier) will automatically revert to its original (i.e., fourier@diff). It is for this reason that my implementation below uses the inert Diff rather than the usual diff.

map2(
    inttrans:-addtable,
    fourier,
    [diff, Diff](U(x,t),t),
(* => *) Diff(inttrans:-fourier(U(_x,t), _x, s), t),
    x, s, {U(x,t), t}
): 
#There's no response if this command is accepted.

Test it:

inttrans:-fourier(diff(U(x,t),t), x, s);  lprint(%);
Diff(fourier(U(_x,t),_x,s),t)

I changed x to _x in the result to show that it's not really the same x. This is akin to @sand15's xi. Since they're bound variables, I don't think it makes any practical difference.

Unfortunately, this only implements your transform for the specific function name U. I may be able to figure out something more general later..

There are eight groups of three elimination variables that can be used, all of which yield the desired polynomial. The process shown in the following worksheet finds all of them quickly and automatically.

restart:

Peng-Robinson_eos

In the following computation, the variables in each numerator are never separated, so we can reduce the number of variables by considering each numerator as a single variable.

eq4:= P = RT/(v-b) - aT/(v*(v+b)+b*(v-b));

P = RT/(v-b)-aT/(v*(v+b)+b*(v-b))

elims:= {(A,B,Z) =~ P/RT *~ (aT/RT, b, v)};

{A = P*aT/RT^2, B = P*b/RT, Z = P*v/RT}

So there are potentially five variables that we can eliminate:

Vs:= indets(rhs~(elims));

{P, RT, aT, b, v}

But since there are only three elimination equations, we need to pick three variables. The possible groups of three are

V3:= [combinat:-choose(Vs, nops(elims))[]];

[{P, RT, aT}, {P, RT, b}, {P, RT, v}, {P, aT, b}, {P, aT, v}, {P, b, v}, {RT, aT, b}, {RT, aT, v}, {RT, b, v}, {aT, b, v}]

Use solve on each group of three variables:

sols:= V3 =~ map2({solve}, elims, V3);

[{P, RT, aT} = {}, {P, RT, b} = {{P = Z^2*aT/(A*v^2), RT = Z*aT/(A*v), b = B*v/Z}}, {P, RT, v} = {{P = B^2*aT/(A*b^2), RT = B*aT/(A*b), v = Z*b/B}}, {P, aT, b} = {{P = Z*RT/v, aT = A*RT*v/Z, b = B*v/Z}}, {P, aT, v} = {{P = B*RT/b, aT = A*RT*b/B, v = Z*b/B}}, {P, b, v} = {{P = A*RT^2/aT, b = B*aT/(A*RT), v = Z*aT/(A*RT)}}, {RT, aT, b} = {{RT = P*v/Z, aT = A*P*v^2/Z^2, b = B*v/Z}}, {RT, aT, v} = {{RT = P*b/B, aT = A*P*b^2/B^2, v = Z*b/B}}, {RT, b, v} = {{RT = RootOf(A*_Z^2-P*aT), b = B*RootOf(A*_Z^2-P*aT)/P, v = Z*RootOf(A*_Z^2-P*aT)/P}}, {aT, b, v} = {{aT = A*RT^2/P, b = B*RT/P, v = Z*RT/P}}]

Obviously, we can't use the one whose solution set is empty. The one with RootOf could also be problematic, so we'll exclude it also, at least for now.

ok_sols:= ((lhs=op@rhs)~@remove)(
    (r-> has(r, RootOf) or nops(r) <> 1)@rhs,
    sols
);

[{P, RT, b} = {P = Z^2*aT/(A*v^2), RT = Z*aT/(A*v), b = B*v/Z}, {P, RT, v} = {P = B^2*aT/(A*b^2), RT = B*aT/(A*b), v = Z*b/B}, {P, aT, b} = {P = Z*RT/v, aT = A*RT*v/Z, b = B*v/Z}, {P, aT, v} = {P = B*RT/b, aT = A*RT*b/B, v = Z*b/B}, {P, b, v} = {P = A*RT^2/aT, b = B*aT/(A*RT), v = Z*aT/(A*RT)}, {RT, aT, b} = {RT = P*v/Z, aT = A*P*v^2/Z^2, b = B*v/Z}, {RT, aT, v} = {RT = P*b/B, aT = A*P*b^2/B^2, v = Z*b/B}, {aT, b, v} = {aT = A*RT^2/P, b = B*RT/P, v = Z*RT/P}]

All eight of the remaining combinations yield the desired residual polynomial:

for V in lhs~(ok_sols) do
    sol[V[]]:= eliminate({eq4,elims[]}, V)[2][]
od;

B^3-3*B^2*Z+B*Z^2+Z^3-A*B+A*Z+B^2-2*B*Z-Z^2

B^3-3*B^2*Z+B*Z^2+Z^3-A*B+A*Z+B^2-2*B*Z-Z^2

B^3-3*B^2*Z+B*Z^2+Z^3-A*B+A*Z+B^2-2*B*Z-Z^2

B^3-3*B^2*Z+B*Z^2+Z^3-A*B+A*Z+B^2-2*B*Z-Z^2

B^3-3*B^2*Z+B*Z^2+Z^3-A*B+A*Z+B^2-2*B*Z-Z^2

B^3-3*B^2*Z+B*Z^2+Z^3-A*B+A*Z+B^2-2*B*Z-Z^2

B^3-3*B^2*Z+B*Z^2+Z^3-A*B+A*Z+B^2-2*B*Z-Z^2

B^3-3*B^2*Z+B*Z^2+Z^3-A*B+A*Z+B^2-2*B*Z-Z^2

(collect(sol[aT,b,v], Z) = 0) &where subs(RT= R*T, elims);

`&where`(Z^3+(B-1)*Z^2+(-3*B^2+A-2*B)*Z+B^3-A*B+B^2 = 0, {A = P*aT/(R^2*T^2), B = P*b/(R*T), Z = P*v/(R*T)})

Verify that the substitutions return the original eq4:

solve(eval(op(1,%), op(2,%)), P);

0, 0, -(R*T*b^2-2*R*T*b*v-R*T*v^2-aT*b+aT*v)/(b^3-3*b^2*v+b*v^2+v^3)

convert(%[3], parfrac, v);

-aT/(-b^2+2*b*v+v^2)+R*T/(v-b)

 

Download Peng_Robinson.mw

AFAIK, pdsolve(..., numericcan only handle PDE systems with 2 differentiation variables. You have X, R, and t.

Here are two ways:

subsindets(expr, specindex(a), x-> b[op(x)])

Or, if you know that all the a's that you want have exactly 1 index:

subsindets(expr, a[anything], x-> b[op(x)])

How do you expect it to react when you give it a y-range that's out-of-bounds (into complex numbers) for the function? Narrow the y-range to -1.53..1.53.

The command FunctionAdvisor will give you a huge amount of information about a special function. Try

FunctionAdvisor(Ei[1](x)).

The notations Ei[a](x) and Ei(a, x) seem to be used interchangeably, and

 Ei(x) = Ei(0, x) = Ei[0](x) = exp(-x)/x.

Here's another way, just to show the possibilities. I often prefer index-free methods, which let me treat matrices and vectors as whole mathematical objects in their own right rather than "poking inside" them.

I construct a matrix of products of powers of X and by taking the outer product of a vector of powers of and a vector of powers of Y. Then I take the elementwise product of that with m. Then use add, which can be used without indices.

(X,Y):= (2,3):  (r,c):= op(1, m) -~ 1:
add(<X^~($0..r)>.<Y^~($0..c)>^%T *~ m);

                           
-11.518

I will readily admit that the above, done on a large scale, will increase your garbage collection time. 

I've done this with respect to your Solusi5, but it can be done for any of your plots.

Colors:= [red, blue, green]:
F:= [A, l, S](t):  tF:= [t, F[]]:
plots:-display(
    plots:-odeplot(
        Solusi5, `[]`~(t, F), t= 0..5, numpoints= 1000,
        color= Colors, legend= [A5, I5, S5]
    ),
    annotation= (typeset @ xcoordinate)(
        proc(t) option remember; uses T= Typesetting;
            (T:-mtable @ op @ (T:-mtr @ T:-mtd @ T:-mtext)~)(
                sprintf~("%a = %a", tF, evalf[5](eval(tF, Solusi5(t)))),
                mathcolor=~ [black, Colors[]], mathsize= 14
            )
        end proc
    )
);

This will initially appear to be a plot just like your version. The difference will become apparent when you drag your mouse over the curves. Instead of the usual box with the x and y coordinates, you will get a box with a column of numbers. The top line gives the value of t at the cursor; the next 3 lines give the values of A, l, and for that value of t. Each of these lines is the same color as its corresponding curve.

There are no Maple help pages for the Typesetting package operators that I used: mtablemtrmtdmtextmathsize, and mathcolor. These are MathML operators which you can find some help for on many independent webpages; then you need to convert them to Maple syntax. There are Maple help pages for typeset and annotation, which aren't part of the Typesetting package. The definition of xcoordinate is on the annotation help page. The sprintf is a standard Maple command, and very similar to the sprintf in C and other languages. The @ is the function composition operator (help page ?@). The ~ is the elementwise operator (help page ?operators,elementwise). The `[]` is the list-building operator in prefix form (a tiny amount of help at ?use). 

Let me know if you'd like this format changed in any way.

Like this:

select(hastype, A01 + _z, specfunc(tanh))

Edit: I added + _z to avoid the problem that acer mentioned.

This quickly returns a fairly simple elliptic expression:

int(1/sqrt((a-t)*(t-b)*(t-c)*(t-d)), t= b..a) assuming d<c, c<b, b<a;

You shouldn't make assumptions about a variable of integration, t in this case; it's done automatically.

The easiest way to get the nontrivial solution is

solve({eq||(1..4)}, {DC, tau, ton, T})

Since there's only one solution, assumptions can't help.

Do

eval(A1, exp= cos)

It seems like the motivation for your Question is doing some preparation for solving your BVP by the shooting method. Acer's implicitplot shows that there are likely 4 solutions. Maple does have a direct numeric solver for BVPs, but it's difficult to use when there are multiple solutions. Shooting works well and fast in this case. The entire worksheet below executes on my computer in less than 1 second, generating plots of y(x) and y'(x) for all 4 solutions.

restart:

st:= time():  #CPUtime at start

#The fsolve command below will likely fail if the dsolve values computed
#at the right boundary are not done at higher precision than that requested
#for the fsolve left-boundary solution.
Digits:= 15:

#
ode:= diff(y(x),x,x) + y(x)^2*diff(y(x),x) + cos(y(x)) = x;
bc:= (1+y)*D(y) - 3*y: #common part of both BCs
bc1:= bc(0);  # "= 0" is implied in both BCs.
bc2:= (bc - y^2)(1) + 1/2;

diff(diff(y(x), x), x)+y(x)^2*(diff(y(x), x))+cos(y(x)) = x

(1+y(0))*(D(y))(0)-3*y(0)

-y(1)^2+(1+y(1))*(D(y))(1)-3*y(1)+1/2

bc1_sol:= solve(bc1, D(y)(0));

3*y(0)/(1+y(0))

ICs:= (ic1:= y(0)= a), D(y)(0) = subs(ic1, bc1_sol);

y(0) = a, (D(y))(0) = 3*a/(1+a)

dsol_ivp:= dsolve({ode, ICs}, numeric, parameters= [a]):

#The syntax in the next procedure requires that y(1) be replaced with y(x)
#in bc2, but this will only be used for x=1.
bc2_prep:= subs([y(1)= y(x), D(y)(1)= diff(y(x),x)], bc2);

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

BC2:= proc(a)
    #Digits is a special type of global variable called an "environment
    #variable", which means that its value be reset to its previous value
    #every time this procedure exits.
    Digits:= 15;
    if a::numeric then
        dsol_ivp(parameters= [a]);
        eval(bc2_prep, dsol_ivp(1))
    else
        'BC2'(a)
    fi
end proc
:

Digits:= 10:  #precision for fsolve

#These ranges come from acer's implicitplot:
A_rngs:= [-4..-3, -2.5..-1.5, -0.99..-0.7, 0..0.3];

[-4 .. -3, -2.5 .. -1.5, -.99 .. -.7, 0 .. .3]

a_sols:= fsolve~(BC2(a), a=~ A_rngs);

[-3.535939019, -2.002235920, -.7734241018, .1155584192]

for a in a_sols do
    dsol_ivp(parameters= [a]);
    (print@plots:-odeplot)(
        dsol_ivp, `[]`~(x, [y(x), diff(y(x), x)]), x= 0..1,
        axes= boxed, view= [0..1.1, DEFAULT], gridlines,
        legend= [y, `y'`], size= [400,200], legendstyle= [location= right]
    )
od:

time()-st;  #CPUtime for this whole worksheet

.766

 

Download Shooting_method.mw

Like this:

restart;
mysphere:= (x - a)^2 + (y - b)^2 + (z - c)^2 - 10^2 = 0;
L:= [seq](seq(seq(`if`(igcd(a,b,c)=1, mysphere, NULL), a= 1..5), b= 1..5), c= 1..5);

 

The command to merge the iterated exponents is combine. For example, starting with the original W, you could do

combine(convert(W, exp));

For future Questions, please note that uploaded worksheets are much more useful to us than uploaded screenshots. In this particular case, the Question is simple enough that a photo was sufficient, but that's rarely the case.

combine can also work directly on your hyperbolic trig functions, so this produces an equivalent result:

convert(combine(W), exp);

3 4 5 6 7 8 9 Last Page 5 of 381