Carl Love

Carl Love

28070 Reputation

25 Badges

13 years, 30 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

In this case, it means that you made a syntax error. The y in your difeq should be y(x). That is how you indicate to Maple that y is the dependent variable and x is the independent variable.

Of course you get errors: You're using Matlab code in Maple. The syntax is completely different. Maple syntax is much more compact and intuitive. Here is how to do it in Maple:

a0:= 80;  a1:= -5;  b1:= -5*(3)^0.5;  w:= Pi/12;  k:= 0.2;
U0:= [65,70,95]:
T:= 0..100:
c0:= u0-a0-(k^2*a1-k*w*b1)/(k^2+w^2);
c1:= (k^2*a1-k*w*b1)/(k^2+w^2);
d1:= (k*w*a1+k^2*b1)/(k^2+w^2);
u:= a0+c0*exp(-k*t)+c1*cos(w*t)+d1*sin(w*t);
plot(
     [seq(u, u0= U0)], t= T, legend= [seq('u'(0)= u0, u0= U0)],
     color= [red, black, green], linestyle= [solid, dash, solid]
);

Using the S and t from your original Question, the idea of Tom Leslie can be implemented like this:

S3:= zip((x,y)-> [x[1],x[2]*y], S, CurveFitting:-ArrayInterpolation(rtable(t), S[..,1]));

By the way, "serie" is not an English word. The word is "series" whether it's singular or plural.

Your understanding of the theory is a bit off. The function has a minimum if the matrix of second derivatives (called the Hessian) is positive definite. This is the same as the eigenvalues of the Hessian being all positive. The determinant of a matrix is the product of its eigenvalues, but it alone will not reveal the signs of the eigenvalues.

The worksheet below answers your problem #1.

 

f:= 4*x^2+y^2+2*z^2+2*a*x*y-4*x*z+2*y*z:
V:= [x,y,z]:

G:= VectorCalculus:-Gradient(f, V):

solve(convert(G,list), V);

[[x = 0, y = 0, z = 0]]

(1)

So (0,0,0) is a critical point regardless of a.

H:= VectorCalculus:-Hessian(f, V);

H := Matrix(3, 3, {(1, 1) = 8, (1, 2) = 2*a, (1, 3) = -4, (2, 1) = 2*a, (2, 2) = 2, (2, 3) = 2, (3, 1) = -4, (3, 2) = 2, (3, 3) = 4})

(2)

The critical point will be a relative minimum if all three eigenvalues of the Hessian are positive.

E:= LinearAlgebra:-Eigenvalues(H):

plot(convert(E,list), a= -10..10, legend= [E1,E2,E3], gridlines= false);

 

There is apparently a small region where all three are positive.

solve(E[2], a);

0, -2

(3)

So (0,0,0) is a relative minimum for a in (-2,0). What happens at a = -2 and a = 0 is unclear. For a outside [-2,0], (0,0,0) is a saddle point.

 

Download 2nd_derivative_test.mw

 

@tomleslie A second alternative:

M:= Matrix(5,`+`);

seq(`*`(L4[1..k]), k= 1..nops([L4]));

There is a scaling problem because the values in X are so close together. You can rescale simply by subtracting 616.3 from each element of X:

g:= expand(subs(x= x+616.3, Statistics:-Fit(a+b*x+c*x^2, X -~ 616.3, Y, x));

Now, you just get a warning message saying that there are zero degrees of freedom, which just means that there are the same number of data points as parameters.

Similar results can be achieved with

CurveFitting:-PolynomialInterpolation(X,Y,x);

If you have a some guesses as to the general form of the curve and you are looking for a curve that approximates the data but does not necessarily pass through every point, use Statistics:-NonlinearFit or Statistics:-LinearFit. It's hard to say more without seeing your data.

This is a very difficult problem, and in Maple several intricate steps are required to solve it. I am amazed that it was assigned as a homework problem. Here's an overview of the steps:

1. Get a parameterized numeric dsolve solution with the initial velocity and and air density as parameters.
2. For any given density and initial velocity, extract the numeric x- and y-coordinate functions from the dsolve solution.
3. Use fsolve on the x-coordinate function to determine the time required for the ball to get to the home-run distance.
4. Put that time into the y-coordinate function to determine whether the ball clears the home-run fence.
5. Repeat steps 2-4 with different initial velocities to determine the minimal initial velocity needed to clear the home-run fence. (Unfortunately, fsolve is too fussy to perform this step. I had to resort to a crude bisection root-finding method.)
6. Repeat steps 2-5 with a different air density.

By Boyle's Law, pressure is proportional to density, so 10% lower air pressure corresponds to 10% lower air density.


restart:

#Parameters (v0 and rho are specified after dsolve)
m:= 0.145:      #mass of baseball
#v0:=           #initial velocity
g:= 9.81:       #gravity
#rho:=          #density of air
C_d:= 0.35:     #drag coefficient
r:= 0.037:      #radius
y0:= 1.5:       #y(0): height from which ball is hit
yf:= 2.5:       #height of home-run fence
xf:= 114:       #home-run distance
theta:= Pi/4.:  #angle of elevation of the hit

v_x:= diff(x(t),t):
v_y:= diff(y(t),t):
eqx:= diff(v_x,t) = -((C_d)*rho*Pi*(r^2)*(v_x)*sqrt((v_x)^2 +(v_y)^2))/(2*m);
eqy:= diff(v_y,t) = -((C_d)*rho*Pi*(r^2)*(v_y)*sqrt((v_x)^2 +(v_y)^2))/(2*m)-g;
ics:= x(0) = 0, y(0) = y0, D(x)(0) = v0*cos(theta), D(y)(0) = v0*sin(theta):
Sol:= dsolve({eqx,eqy,ics}, numeric, parameters= [v0, rho], output= listprocedure):

#Procedure to compute the minimal initial velocity needed to clear the fence given
#a certain density of air.
HomerunVelocity:= proc(rho)
local
     v0,
     ClearTheFence:= proc(v0)
     local X,Y,r;
          if not v0::numeric then return 'procname'(args) end if;
          Sol(parameters= [:-v0= v0, :-rho= rho]);
          #Extract numeric procedures for X and Y
          (X,Y):= eval([x,y](t), Sol)[]:
          #What is the height above the fence when the distance is xf?
          Y(fsolve(X(t)=xf)) - yf
     end proc
;
     #fsolve has a lot of trouble trying to solve ClearTheFence(v0)=0.
     #Thus, we resort to the crude bisection method.
     evalf[6](
          Student:-NumericalAnalysis:-Bisection(
               ClearTheFence(v0), [30,100], tolerance= 1e-6
          )
     )
end proc:  

diff(diff(x(t), t), t) = -0.5190669380e-2*rho*(diff(x(t), t))*((diff(x(t), t))^2+(diff(y(t), t))^2)^(1/2)

diff(diff(y(t), t), t) = -0.5190669380e-2*rho*(diff(y(t), t))*((diff(x(t), t))^2+(diff(y(t), t))^2)^(1/2)-9.81

v_sealevel:= HomerunVelocity(1.2);

46.6597

#Try with air 10% less dense.
v_denver:= HomerunVelocity(1.2*.9);

44.9244

Sol(parameters= [v0= v_sealevel, rho= 1.2]):

P1:= plots:-odeplot(Sol, [[x(t),y(t)]], t= 0..6, color= red):

Sol(parameters= [v0= v_denver, rho= 1.2*0.9]):

P2:= plots:-odeplot(Sol, [[x(t),y(t)]], t= 0..6, color= green):

plots:-display([P1,P2], view= [0..xf, DEFAULT]);

 

 

 

Download homerun.mw

 

Try using a restart command. I am able to execute your Explore command without any problem.

Use an epsilon (I used 1 below) to change all your strict inequalities to nonstrict. Then pass them to Optimization:-LPSolve with any objective function (I used 0 below). 

e:= 1:
Optimization:-LPSolve(
     0,
     {x1 >= 0+e, x1+x2 >= 300+e, x2 <= 1000-e,
     x1+x2 <= 700-e, x2 >= 0+e, x1 <= 1000-e}
);

I think that my approach to making a color bar legend is sufficienty different from those found in the links that Alejandro provided that it is worth posting; however, it suffers the same drawback as all the others---that Maple's array plotting function allocates equal space to the color bar and to the main plot. The major difference is how I use plottools:-transform to extract the z-range from the plot.

First, make your plot and assign it to a variable. It can be any 3D plot produced with shading= zhue; it does not need to be produced with plot3d.

plot3d(x+y, x= 0..1, y= 0..1, shading= zhue);
P:= %:

(Plot omitted for brevity.) The following appliable module extracts the minimum and maximum z-coordinates from a 3D plot.

minMax:= module()
local min, Max;
export
     filter:= proc(x,y,z)
          if z::numeric then
               if z < min then min:= z end if;
               if z > Max then Max:= z end if
          end if;
          [x,y,z]
     end proc,
     ModuleApply:= proc(P::specfunc(anything, PLOT3D))
          min:= infinity;  Max:= -infinity;
          plottools:-transform(filter)(P);
          `..`(min,Max)
     end proc
;
end module;

Now use the plot's data to make the color bar:

plots:-tubeplot(
     [0,0,z], z= minMax(P),
     shading= zhue, scaling= constrained, radius= .2, axes= frame,
     tubepoints= 3, style= patchnogrid, orientation= [90,90],
     axis[1,2]= [tickmarks= 0], lightmodel= none
);
colorbar:= %:

Now use an array plot to put the two next to each other:

plots:-display(< P | colorbar >);

(Unfortunately, MaplePrimes cannot handle array plots.)

 

 

 

The following will produce a very smooth curve that fits the data rather closely.

m:= min(y_data): M:= max(y_data):
Y:= [seq(m..M, (M-m)/100)]:
X:= CurveFitting:-ArrayInterpolation(y_data, x_data, Y, degree= 3, method= spline):
plot(zip(`[]`, Y, X));

Since you said that you don't quite understand the math behind it, here it is in more detail.

F:= (x,y)-> x*(x+y)*exp(y-x):

We first note that the function is differentiable for all x and y.

Compute all first partial derivatives.
dFdx:= simplify(D[1](F)(x,y));

dFdy:= simplify(D[2](F)(x,y));

Set first partial derivatives to 0 and solve simultaneously. In the language of analytic geometry, this corresponds to finding where the tangent plane is horizontal. In the language of vector calculus, we would say that these points are where the gradient is 0. These points are called critical points.

CritPts:= [solve({dFdx=0, dFdy=0}, {x,y})];

Now we need to classify the critical points as minima, maxima, or saddle points. Algebraically, the easiest way to do this is to use the second derivatives. Since we have a computer algebra system, that's the way that we'll use.

Compute the matrix of second partial derivatives, which is called the Hessian.

Hess:= Matrix(2, (i,j)-> simplify(D[i,j](F)(x,y)));

We evaluate the Hessian at each critical point and then compute the eigenvalues. If all the eigenvalues are positive, the point is a minimum. If all the eigenvalues are all negative, it's a maximum. If some eigenvalues are positive, some are negative, and none are zero, then it's a saddle point. If any eigenvalues are zero, the test is inconclusive.

for cp in CritPts do convert(evalf(LinearAlgebra:-Eigenvalues(eval(Hess, cp))), list) od;
                  [2.414213562, -0.414213562]
                  [0.3543123712, 0.0516934784]

So (0,0) is a saddle point and (1/2, -3/2) is a minimum.

 

(Man, I'm sick and tired of posting whole worksheets by cutting and pasting each paragraph separately! Please fix this so that I can upload whole worksheets! Sometimes it works for me, sometimes not.)

The command Maplets:-Elements:-TextField has an option tooltip= "...". That's what you want. See ?Maplets,Elements,TextField

First 273 274 275 276 277 278 279 Last Page 275 of 395