Carl Love

Carl Love

28035 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

What you are calling the "gradient" of a vector field is called the Jacobian in English. You simply need to change Gradient to Jacobian.

In English, as far as I know, one only computes the gradient of functions from R^n to R^1. The Jacobian is a generalization of the gradient to functions from R^n to R^m. The distinction is only linguistic and a matter of convention; there's no practical difference between first derivatives, gradients, and Jacobians.

When I used what seemed to me to be more-realistic values for the constants, Maple wasn't as compliant as it was for Kitonum's values of the constants.

Numeric integration over a hyper-rectangle is more robust than integration over regions with variable inner limits of integration. I can see immediately that your region of integration is a cylinder of height t over the unit disk. So a conversion to cylindrical coordinates will make the region of integration a hyper-rectangle.


 

restart:

Digits:= 15: #15 is usually the best (fastest, most robust) value.

#Your integrand, just like you entered it:
f:= exp(-1/2*(xop^2+yop^2))*exp(-((x-xop-2*Pe*(t-tp))^2+(y-yop)^2+z^2)/(4*(t-tp)))/((4*Pi^(3/2))*(t-tp)^(3/2));

(1/4)*exp(-(1/2)*xop^2-(1/2)*yop^2)*exp(-((x-xop-2*Pe*(t-tp))^2+(y-yop)^2+z^2)/(4*t-4*tp))/(Pi^(3/2)*(t-tp)^(3/2))

(1)

#Express as a triple integral with nested, inert (capital-I) Int, without limits of integration:
I3D:= Int(Int(Int(f, yop), xop), tp);

Int(Int(Int((1/4)*exp(-(1/2)*xop^2-(1/2)*yop^2)*exp(-((x-xop-2*Pe*(t-tp))^2+(y-yop)^2+z^2)/(4*t-4*tp))/(Pi^(3/2)*(t-tp)^(3/2)), yop), xop), tp)

(2)

#Convert to cylindrical coordinates:
Icyl0:= IntegrationTools:-Change(I3D, {yop= r*sin(theta), xop= r*cos(theta)});

Int(Int(Int((1/4)*exp(-(1/4)*(2*r^2*cos(theta)^2*t-2*r^2*cos(theta)^2*tp+2*r^2*sin(theta)^2*t-2*r^2*sin(theta)^2*tp+r^2*cos(theta)^2+4*cos(theta)*Pe*r*t-4*cos(theta)*Pe*r*tp+r^2*sin(theta)^2+4*Pe^2*t^2-8*Pe^2*t*tp+4*Pe^2*tp^2-2*cos(theta)*r*x-2*sin(theta)*r*y-4*Pe*t*x+4*Pe*tp*x+x^2+y^2+z^2)/(t-tp))*abs(cos(theta)^2*r+r*sin(theta)^2)/(Pi*(t-tp)*((t-tp)*Pi)^(1/2)), r), theta), tp)

(3)

#Extract the integrand:
Fcyl1:= op([1,1,1], Icyl0);

(1/4)*exp(-(1/4)*(2*r^2*cos(theta)^2*t-2*r^2*cos(theta)^2*tp+2*r^2*sin(theta)^2*t-2*r^2*sin(theta)^2*tp+r^2*cos(theta)^2+4*cos(theta)*Pe*r*t-4*cos(theta)*Pe*r*tp+r^2*sin(theta)^2+4*Pe^2*t^2-8*Pe^2*t*tp+4*Pe^2*tp^2-2*cos(theta)*r*x-2*sin(theta)*r*y-4*Pe*t*x+4*Pe*tp*x+x^2+y^2+z^2)/(t-tp))*abs(cos(theta)^2*r+r*sin(theta)^2)/(Pi*(t-tp)*((t-tp)*Pi)^(1/2))

(4)

#Simplify the integrand:
Fcyl2:= simplify(Fcyl1) assuming r >= 0, r <= 1;

(1/4)*exp(-(1/4)*(4*cos(theta)*Pe*r*t-4*cos(theta)*Pe*r*tp+4*Pe^2*t^2-8*Pe^2*t*tp+4*Pe^2*tp^2-2*cos(theta)*r*x-2*sin(theta)*r*y-4*Pe*t*x+4*Pe*tp*x+2*r^2*t-2*r^2*tp+r^2+x^2+y^2+z^2)/(t-tp))*r/(Pi^(3/2)*(t-tp)^(3/2))

(5)

#Express as the final symbolic intregral (with limits of integration):
Icyl2:= Int(Fcyl2, [theta= 0..2*Pi, r= 0..1, tp= 0..t]);

Int((1/4)*exp(-(1/4)*(4*cos(theta)*Pe*r*t-4*cos(theta)*Pe*r*tp+4*Pe^2*t^2-8*Pe^2*t*tp+4*Pe^2*tp^2-2*cos(theta)*r*x-2*sin(theta)*r*y-4*Pe*t*x+4*Pe*tp*x+2*r^2*t-2*r^2*tp+r^2+x^2+y^2+z^2)/(t-tp))*r/(Pi^(3/2)*(t-tp)^(3/2)), [theta = 0 .. 2*Pi, r = 0 .. 1, tp = 0 .. t])

(6)

#Supply numeric values for the 5 constants:
Icyl3:= eval(Icyl2, [t= 1/2, x= 1/2, y= 1/2, z= 1/2, Pe= 1/2]);

Int((1/4)*exp(-(1/4)*(-2*cos(theta)*r*tp+1/2+tp^2-r*sin(theta)+2*r^2-2*r^2*tp)/(1/2-tp))*r/(Pi^(3/2)*(1/2-tp)^(3/2)), [theta = 0 .. 2*Pi, r = 0 .. 1, tp = 0 .. 1/2])

(7)

#Perform numeric integration:
evalf(%);

.174309703986304

(8)

 

Download NumIntCylCoord.mw

I can confirm that value to 9 digits by using the original coordinates, performing the innermost integral symbollically, and then the two outer integrals numerically. I don't think that this method is as robust numerically as the conversion to cylindrical coordinates.

A conditional solution will either be a piecewise branch with an inequality condition (if you designated which names were to be considered parameters), or it'll contain inequalities (if you didn't). The presented case is the latter. To separate the solutions which contain inequalities from those that don't:

S:= solve({2*x+y+z = 3, x^2+y^2+z^2 = 3}, [x, y, z], real, parametric);
selectremove(x-> membertype({`<`, `<=`}, x), S);

To just get the unconditional solutions, change selectremove to remove.

 

For what it's worth:

I changed your initial condition so that it didn't reuse the independent variable x. So, I made it u(a, 1/a) = 0. This eliminated the error message, so I guess that that is the preferred form for initial conditions; however, I let the pdsolve command run 11 hours without getting any result, then I killed it. Like a zombie, it started running again. I eventually needed to kill my entire Maple session. This is in Maple 2016.

If you include a trivial inequality in the system, solve seems spurred to return only real solutions without the need to specify any other options, load any packages, or use any assumptions. This is even true if the inequality has no relation to the other equations in the system, as in

solve({x^2+y^2+z^2=3, x+y+z=3, w > 0}, {x,y,z});

     {x = 1, y = 1, z = 1}

I obtained the above result both with and without RealDomain.

Here's the procedure npar:

npar:= proc(p::procedure, $)
local L:= [op(1, eval(p))];
   nops(L) - `if`(` $` in L, 1, 0)
end proc;

This counts the parameters, not the arguments. What's the distinction? Parameters are an inherent property of a procedure; they're part of how it's written. Arguments are what's used when the procedure is called. The number of arguments to a given procedure may vary. If I declare

f:= x-> x^2;

then f has one parameter, x. If I call it using f(1,2), then it has two arguments for this call. This will not generate an error; by default it's allowed to pass extra unused arguments to a procedure. If you want to insist that f be called with exactly one argument, then you need to write it like this:

f:= proc(x,$) x^2 end proc:

Now that $ will appear as ` $` in op(1, eval(f)), which is why my procedure npar checks for it and subtracts 1 if it's found.

Inside procedure lich, you use f as if it were itself a procedure. That's fine, but then you need to pass a procedure to f, not an expression. So, you just need to change the call to lich to

lich(x-> sin(x^2), 0, Pi, 10);

We usually don't take kindly to screenshots. In the future, please upload a worksheet using the green uparrow on the toolbar of the MaplePrimes editor.

To get longitude and latitude lines on a sphere, use plot3d with coords= spherical:

plot3d(
   [1, phi, theta], theta= 0..2*Pi, phi= -Pi/2..Pi/2,
   coords= spherical, scaling= constrained
);

You can control the number of lines (to some extent) with the grid= [m, n] option.

Here's a worksheet that summarizes all that I've said so far, and also works out some of the finer details and nuisances of implementation.


 

restart:

Digits:= 15: #15 is usually best (fastest, most robust) value.

f:= (x::numeric, a::numeric, b::numeric, c::numeric)->
   a*I + b*x*exp(c*I):

g:= (x::numeric)->
   7*I + 5*x*exp(3*I):

So, we expect the result a = 7, b = 5, c = 3. But note that f(a, -b, c-Pi) = f(a, b, c), so there are equivalent results that are just as valid.

(A,B):= (-2,2): #interval of x-values (arbitrarily chosen)

H:= (a::numeric, b::numeric, c::numeric)->
   Re(int(
     (x-> abs(f(x,a,b,c) - g(x))^2), A..B,
     numeric, digits= iquo(Digits,2)
   ))
:

#In H:
#1. Re is needed because the result is sometimes returned as x + 0.*I.
#2. We must use operator form to avoid evaluations with symbolic x.
#3. The squared norm seems to work better.
#4. Results will be computed to 7 digits, and there'll be 8 guard digits (15-7 = 8).

infolevel[Optimization]:= 5:
#Ask the minimizer to reveal some of its internal process.

CodeTools:-Usage(Optimization:-Minimize(H));
#Find (a,b,c) that minimize H, and get timing stats.

NLPSolve: calling NLP solver

NLPSolve: using method=pcg
NLPSolve: number of problem variables 3
NLPSolve: optimality tolerance set to 0.3256082241e-11
NLPSolve: iteration limit set to 50
NLPSolve: trying evalhf mode
NLPSolve: trying evalf mode

attemptsolution: number of major iterations taken 10
memory used=0.52GiB, alloc change=428.01MiB, cpu time=5.56s, real time=5.66s, gc time=234.38ms

 

[3.66534600000000022*10^(-16), Vector(3, {(1) = 7.00000000054163, (2) = -5.00000000065577, (3) = -.141592653649165})]

(1)

This says that a = 7, b = -5, c = -.14159 (3 - Pi). The value on the left is the minimal value of the norm.

 


 

Download Complex_norm_minimization.mw

If you want the number to display with the number of decimal places with which it was entered, but you don't want to explicitly specify that number of places, then use

FloatString:= (x::float)-> sprintf("%.*f", -op(2,x), x);
FloatString(0.0141);

To either convert command, append the clause assuming t1 > 0.

If your expression for drag force is df, then all that you need is

Cd:= V-> df *~ V;

If you have more trouble, please post a Maple worksheet (not a .docx).

Whether you use a colon or a semicolon to terminate a statement in a loop makes no difference. The only one that matters is the one after the final end do.

The loop depth for which output is displayed is controlled by the environment variable printlevel. By default, this is set to 1, so only the output of unnested loops is displayed. If you increase it:

printlevel:= 2;

then the output of every command in your doubly nested loops will be shown. I strongly advise that you only use this technique for debugging because, in general, the volume of output will be overwhelming.

The proper way to display the output of commands within loops or other control structures is to use print. That is, for any X that you want displayed, use print(X).

To return the entry of minimal modulus in a single line of code, use

rts[min[index](abs~(rts))];

It is also possible (in a very short single line of code) to sort the entries by modulus. This is especially useful for lists of eigenvalues.

This can also be done with simplify with side relations:

sol:= dsolve({diff(y(t),t$2) + w^2*y(t) = p/m*exp(-a*t), y(0) = 0, D(y)(0) = 0}):
simplify(sol, {a= n*w, p= m*w^2*u});

 

First 187 188 189 190 191 192 193 Last Page 189 of 395