acer

32343 Reputation

29 Badges

19 years, 327 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@Rouben Rostamian  There are a few things I find alarming about the way fsolve "works" with units. As Axel showed, the following works (Maple 2021.2, for me) with Units:-Simple loaded and in play.

   f:= r -> eqn2(6*Unit('degree'), r*Unit('m'/'sec'))/Unit(N);
   fsolve(f, 7 .. 8);

However, that takes about 11sec on my Intel Core i5-7400 @3.00GHz. That's not good.

One problem is that each call to int(...,numeric) involves a complicated integrand with a mess of subexpressions with unsimplified units all round. Those all get combined and simplified at each call. And dimensional checks are also done repeatedly.

There's also some memoization of the units analysis, so that exact same fsolve attempt returns quickly when done a second time. (That's also scary to me -- what if increases in working precision were ignored by that!)

So here's a revision in which I've avoided loading Units (and thus also Units:-Simple) altogether. A judiciously placed and terse invocation of a few things like evalf(Int(...)) instead of int(...,numeric), change of variables to clear the units out of the range of integration, combine(...units) , etc, can still get the job done.

On my machine it's about 250 times faster to compute it this way.

But I didn't edit any of the definitions prior to that of eqn2, and even in a revision of eqn2 I only switched int to Int. The "unit handling" was only added under the procedures passed to fsolve.

restart;

kernelopts('version');

`Maple 2022.1, X86 64 LINUX, May 26 2022, Build ID 1619613`

with(LinearAlgebra):

#with(Units):

N__bld := 4:

D__rtr := 30*Unit('ft'):
R__tip := D__rtr/2:
chd := 8*Unit('inch'):
R__hng := 1*Unit('ft'):

rho__air:=ThermophysicalData:-Atmosphere(10,density,useunits):

a__air:=ThermophysicalData:-Atmosphere(10, speedofsound, useunits):

M__tip := 0.62:

omega:=(M__tip*a__air)/R__tip:      

Cl__alpha := 0.1*Unit(1/'deg'):
C__L := alpha -> alpha*Cl__alpha:

I converted NACA0012_CD_alpha.xlsx to NACA0012_CD_alpha.csv, and applied

the following commands to determine a sixth degree polynomial fit.

I have pasted the resulting polynomial here, so you don't need the raw data file.

(*
data := Import("/tmp/NACA0012_CD_alpha.csv"):
Data := convert(data, 'Matrix'):
alp := LinearAlgebra:-Column(Data, 1) *~ Unit('deg');
cd := LinearAlgebra:-Column(Data, 2);
fit_me := Statistics:-PolynomialFit(6, Data, x);
*)
fit_me := 0.00340780072932875 + 0.0000826700291382254*x + 0.000146289994961542*x^2 - 1.63712719600353*10^(-6)*x^3 - 1.31704118799963*10^(-6)*x^4 + 7.08093026284748*10^(-9)*x^5 + 6.41167350162788*10^(-9)*x^6;

0.340780072932875e-2+0.826700291382254e-4*x+0.146289994961542e-3*x^2-0.1637127196e-5*x^3-0.1317041188e-5*x^4+0.7080930263e-8*x^5+0.6411673502e-8*x^6

#plot(fit_me, x=-15..15);

C__D := unapply(fit_me, x):

omega:

Note unit specification of the arctan.  Angles are measured in radians in Maple.
That may be different in MapleFlow. Check!

alpha__i := (r, theta, u) -> theta - arctan(u/(omega*r))*Unit('radian') - 8*Unit('deg')*r/R__tip + 6*Unit('deg'):

lf__s := (r, theta, u) -> 1/2*rho__air*chd*(omega^2*r^2 + u^2)*C__L(alpha__i(r, theta, u)):

dr__s := (r, theta, u) -> 1/2*rho__air*chd*(omega^2*r^2 + u^2)*C__D(alpha__i(r, theta, u)):

dr__s(r, theta, u):

lf__s(r, theta, u)*cos(alpha__i(r, theta, u)) - dr__s(r, theta, u)*sin(alpha__i(r, theta, u)):
th__s := unapply(%, r, theta, u):

dr__s(r, theta, u)*cos(alpha__i(r, theta, u)) + lf__s(r, theta, u)*sin(alpha__i(r, theta, u)):
trq__s := unapply(%, r, theta, u):

#Thust := (theta, u) -> N__bld *
#        int(th__s(r, theta, u), r = R__hng .. R__tip, numeric);

Thust := proc(theta, u)
        #print(theta, u);
        (N__bld * int(th__s(r, theta, u), r = R__hng .. R__tip, numeric));
end proc:

UD := u -> 2*rho__air*u^2*Pi*R__tip^2;

proc (u) options operator, arrow; 2*rho__air*u^2*Pi*R__tip^2 end proc

eqn2 := (theta, u) -> Thust(theta, u) - UD(u);

proc (theta, u) options operator, arrow; Thust(theta, u)-UD(u) end proc

Questions

 

eqn2(theta,u) may evaluated at arbitrary theta and u values:

combine(eqn2(6*Unit('degree'), 5*Unit('m'/'sec')),units);

9351.119026*Units:-Unit(N)

combine(eqn2(6*Unit('degree'), 10*Unit('m'/'sec')),units);

-9476.904304*Units:-Unit(N)

We see that eqn2(6,u) changes sign as u goes from 5 to 10, so

so it is zero for some u between 5 and 10.
Question 1: How do we find that u?

The following (which I have commented out) does not work:

# fsolve('eqn2'(6*Unit('degree'), u*Unit('m'/'sec')), u = 5 .. 10);

f:= r -> eqn2(6*Unit('degree'), r*Unit('m'/'sec'))/Unit(N):
CodeTools:-Usage( fsolve(r->combine(f(r),units), 7 .. 8) );

memory used=1.34GiB, alloc change=76.00MiB, cpu time=12.00s, real time=10.56s, gc time=2.19s

7.749740189

ThustB := proc(theta, u)
        #print(theta, u);
        (N__bld * Int(th__s(r, theta, u), r = R__hng .. R__tip));
end proc:
eqn2B := (theta, u) -> ThustB(theta, u) - UD(u);

proc (theta, u) options operator, arrow; ThustB(theta, u)-UD(u) end proc

II:=simplify(combine(IntegrationTools:-Change(eqn2B(6*Unit('degree'), u*Unit('m'/'sec')),
                                              r=rr*Unit(ft), rr), units)/Unit(N)):
IIf:=subs(__dummy=II,proc(u) local res;
  Digits:=15;
  res:=evalf(__dummy);
  evalf[10](res);
end proc):

CodeTools:-Usage( fsolve(IIf, 7..8) );

memory used=4.08MiB, alloc change=0 bytes, cpu time=40.00ms, real time=40.00ms, gc time=0ns

7.749740188

Question 2: Suppose that we figure out how to answer the previous question.

How do we generalize it to make it work for unspecified values of the first argument,

that is, find U so that:

#U := theta -> fsolve(eqn2(theta*Unit('degree'), u*Unit('m'/'sec')), u = 0..100);

III:=simplify(combine(IntegrationTools:-Change(eqn2B(theta*Unit('degree'), u*Unit('m'/'sec')),
                                              r=rr*Unit(ft), rr), units)/Unit(N)):
IIIf:=theta->subs(__dummy=eval(III,:-theta=theta),proc(u) local res;
  Digits:=15;
  res:=evalf(__dummy);
  evalf[10](res);
end proc):

U:=theta->fsolve(IIIf(theta), 10^(-Digits+1)..100):

CodeTools:-Usage( U(6) );

memory used=5.85MiB, alloc change=0 bytes, cpu time=59.00ms, real time=59.00ms, gc time=0ns

7.749740188

CodeTools:-Usage( U(8.3) );

memory used=5.23MiB, alloc change=0 bytes, cpu time=48.00ms, real time=48.00ms, gc time=0ns

9.618205391

 

Download rotor_equation_acc.mw

If you upload and attach here your corrupted worksheet then we could try and repair it. You can use the green up-arrow in the Mapleprimes editor for that.

It doesn't help to duplicate your message in several other old Question threads.

@Rouben Rostamian  Alternatively,

f := u -> int(cos(1+arctan(u,r)), r=1..2, numeric) - u^2:

fsolve(f, 0..1);

       0.4910845909

fsolve('f'(u), u=0..1);

       0.4910845909

This can be managed by avoiding evaluation of a symbolic function call such as f(u), thus f treated like a black box.

I haven't investigated yet as to the precise nature of the difficulty.

@MANUTTM I find it difficult to tell exactly what you mean, in this new and related query.

Are you trying to say that you want to compute the globally best D as a function of inputs g and Sr? If so, then is y also an input of that, or is the optimization allowed to vary y over its stated range?

Also, did you see the result of simplification of the Hessian Matrix in my revision of your worksheet above? Do you know how that impacts on your convexity query?

I have deleted three duplicate postings of this. Please stop submitting duplicates of it.

Also, if you have followup details on this query then please add them here, instead of spawning a separate Question thread.

The data of this Question has been changed, so I've deleted my Answer and my followup alternative since they pertained to that particular form of data.

The original Question referrred to only the following data:

and asked about how to produce an image like this:

@Carl Love You example works in Maple 2022.1, but not in all older versions. Below I mention some edits that allow it to work more as intended in older versions (that supported Typesetting).

That text denoting the color, eg. green in your code, should be within escaped string-quotes, ie, \"green\" .

The hex variants should be prefixed by #, rather than x. (In Maple 2022.1 there must be some preceding character(s), and it needn't specifically be x or #)

Also, use in some plotting situations, that name constructed by nprintf should be terminated by the semicolon character.

Maple 2022.1 respects the leading/trailing blanks, but it's trickier in some older versions. This is one reason I preferred (here) the Typesetting exports rather than the MathML style name-markup. However the OP has not yet told us precisely how he intends on using these colored things, nor his Maple version, and so I cannot know which is best here.

@Ronan If you are using 2D Input mode then the anonymous procedure can be delimited by brackets.

Eg, you can paste in this:

   (add@map)((proc(E) ``(lcoeff(E))*normal(E/lcoeff(E)); end proc)@add, [Q(expand(pn))])

That "end" goes with the "proc".

I have deleted a duplicate of this topic.

Put your followup queries and details on this here, not in separate Question threads.

@mehdibgh Yes, fsolve is having problems finding/accepting the multivariate roots if the ranges's boundary is very close.

restart

EQ := Matrix(4, 1, {(1, 1) = 32.1640740637930*Tau[1]-0.172224519601111e-4*Tau[2]-0.270626540730518e-3*Tau[3]+0.1570620334e-9*P[1]+0.3715450960e-14*sin(t), (2, 1) = -0.172224519601111e-4*Tau[1]+32.1667045885952*Tau[2]+0.587369829416537e-4*Tau[3]-0.1589565489e-8*P[1]+0.1004220091e-12*sin(t), (3, 1) = -0.270626540730518e-3*Tau[1]+0.587369829416537e-4*Tau[2]+32.1816411689934*Tau[3]-0.7419658527e-8*P[1]+0.5201228088e-12*sin(t), (4, 1) = 0.1570620334e-9*Tau[1]-0.1589565489e-8*Tau[2]-0.7419658527e-8*Tau[3]+601.876235436204*P[1]})

V := Matrix(1, 4, {(1, 1) = Tau[1], (1, 2) = Tau[2], (1, 3) = Tau[3], (1, 4) = P[1]})

q := 0

X := Matrix(4, 1, {(1, 1) = -0.1156532164e-15*sin(t), (2, 1) = -0.3121894613e-14*sin(t), (3, 1) = -0.1616209235e-13*sin(t), (4, 1) = -0.2074537757e-24*sin(t)})

t := 0.5617399480550e-5

fsolve({seq(EQ[r, 1], r = 1 .. 4)}, {seq(V[1, r] = (1+0.1e-8)*X[r, 1] .. q, r = 1 .. 4)})

{P[1] = -0.1165350731e-29, Tau[1] = -0.6496703177e-21, Tau[2] = -0.1753692918e-19, Tau[3] = -0.9078892919e-19}

``

Download SoalNewton2_ac.mw

You haven't stated why an optimization approach is not acceptable.

@dharr For some reason I had convinced myself that the solution of the linear problem lay outside the OP's provided ranges -- hence my least-squares approach. I might not have looked properly at the range end-points X that the OP provided.

@gravatarq I deleted a duplicate posting of this.

@mehdibgh Your description is that of an inverted range, since with your t value as 1 the values in X will be negative and you have the lower range value q as 0.

Also, you have the same kind of indexing mistake with X, which you have as a 4x1 Matrix, and not a Vector.

You could try the ranges as,

     { seq(V[1, r] = X[r, 1] .. q, r = 1 .. 4) }

but you might not find a solution.

[edit] Perhaps you might actually be satisfied with a solution to an optimization problem of minimizing some related objective (eg. least-squares). Here I raise the working precision to try and get a more accurate result:

restart

EQ := Matrix(4, 1, {(1, 1) = 32.1640740637930*Tau[1]-0.172224519601111e-4*Tau[2]-0.270626540730518e-3*Tau[3]+0.1570620334e-9*P[1]+0.3715450960e-14*sin(t), (2, 1) = -0.172224519601111e-4*Tau[1]+32.1667045885952*Tau[2]+0.587369829416537e-4*Tau[3]-0.1589565489e-8*P[1]+0.1004220091e-12*sin(t), (3, 1) = -0.270626540730518e-3*Tau[1]+0.587369829416537e-4*Tau[2]+32.1816411689934*Tau[3]-0.7419658527e-8*P[1]+0.5201228088e-12*sin(t), (4, 1) = 0.1570620334e-9*Tau[1]-0.1589565489e-8*Tau[2]-0.7419658527e-8*Tau[3]+601.876235436204*P[1]})

V := Matrix(1, 4, {(1, 1) = Tau[1], (1, 2) = Tau[2], (1, 3) = Tau[3], (1, 4) = P[1]})

q := 0

X := Matrix(4, 1, {(1, 1) = -0.1156532164e-15*sin(t), (2, 1) = -0.3121894613e-14*sin(t), (3, 1) = -0.1616209235e-13*sin(t), (4, 1) = -0.2074537757e-24*sin(t)})

t := 1

Digits := 60; Sol := Optimization:-Minimize(add(EQ[r, 1]^2, r = 1 .. 4), seq(V[1, r] = X[r, 1] .. q, r = 1 .. 4))

[0.23051887707880648450227137668244873934400e-43, [P[1] = -0.174566332911092633926811253720853832324936144346344610109133e-24, Tau[1] = -0.973188258955118595763778228079271664727641452276400379153389e-16, Tau[2] = -0.262698373446757694397996566086762382450096865275922610909883e-14, Tau[3] = -0.135999317663106703497601318807782949900124793323219378929709e-13]]

Sol[2]

[P[1] = -0.174566332911092633926811253720853832324936144346344610109133e-24, Tau[1] = -0.973188258955118595763778228079271664727641452276400379153389e-16, Tau[2] = -0.262698373446757694397996566086762382450096865275922610909883e-14, Tau[3] = -0.135999317663106703497601318807782949900124793323219378929709e-13]

evalf[5](evalf[60](eval(add(EQ[r, 1]^2, r = 1 .. 4), Sol[2])))

0.23052e-43

NULL

Download SoalNewton_acc.mw

@fkohlhepp I have deleted a duplicate of this Question thread.

If you have any additional details or related followup queries then feel free to add them here. Please don't spawn a separate Question thread for this.

And here is an animation, based upon the OP's revised code...

I let it be animated w.r.t. the parameter t passed to CircleParm for point P1.

The procedure F produces each plotted frame.

(It could also be done with Explore, optionally with additional parameters passed for F's other arguments -- being the values passed to CircleParm for P2,P3,P4.)

restart;

with(plots): with(plottools):with(LinearAlgebra): unprotect(D);

_EnvHorizontalName := 'x':
_EnvVerticalName := 'y':

Vdot := proc (U, V)local i: add(U[i]*V[i], i = 1 .. 2) end proc:
dist := proc (M, N) sqrt(Vdot(expand(M-N), expand(M-N))) end proc:
MinDistPoint := proc (A::[algebraic, algebraic], B::[algebraic, algebraic], P::[algebraic, algebraic])
local D, R, V;
description "Point on line AB at minimum distance from P";
D := `<,>`(A-B); R := `<,>`(`<|>`(0, -1), `<|>`(1, 0));
V := `<,>`(`<,>`(P) . D, -R . `<,>`(A) . `<,>`(B)); `~`[`/`]([V . D, R . V . D], ` $`, D . D)
end proc:

EqBIS := proc(P, U, V) local a, eq1, M1, t, PU, PV, bissec1; a := (P - U)/LinearAlgebra:-Norm(P - U, 2) + (P - V)/LinearAlgebra:-Norm(P - V, 2); M1 := P + a*t; eq1 := op(eliminate({x = M1[1], y = M1[2]}, t)); return op(eq1[2]); end proc:

Cen := proc(M, N, R) local eq1, eq2, sol; eq1 := EqBIS(M, N, R) = 0; eq2 := EqBIS(N, M, R) = 0; sol := simplify(solve({eq1, eq2}, {x, y})); return [subs(sol, x), subs(sol, y)]; end proc:

CircleParm := t -> [(-t^2 + 1)/(t^2 + 1), 2*t/(t^2 + 1)]:

 

F := proc(parm1,parm2,parm3,parm4)
local P1,P2,P3,P4,r1,r2,r3,r4,p1,p2,p3,p4,Cir1,Cir2,Cir3,Cir4,
      H,C4,C1,C2,C3,J,K,Pts,Ii,Poly,SegP1P4,SegP2P3,quadri,rect,tex;
try
P1 := Vector(CircleParm(parm1)):
P2 := Vector(CircleParm(parm2)):
P3 := Vector(CircleParm(parm3)):
P4 := Vector(CircleParm(parm4)):
p1 := convert(P1, list);
p2 := convert(P2, list);
p3 := convert(P3, list);
p4 := convert(P4, list);
C4 := Vector(Cen(P1, P2, P3)):H := convert(C4, list):
C1 := Vector(Cen(P4, P2, P3));
J := convert(C1, list);
C2 := Vector(Cen(P4, P1, P3));
Ii := convert(C2, list);
C3 := Vector(Cen(P4, P2, P1));
K := convert(C3, list);
Pts := [P1, P2, P3, P4, C1, C2, C3, C4]:
r4 := dist(MinDistPoint(p1, p2, H), H);
r2 := dist(MinDistPoint(p1, p4, Ii), Ii);
r1 := dist(MinDistPoint(p2, p3, J), J);
r3 := dist(MinDistPoint(p2, p1, K), K);
Cir1 := circle(Ii, r2, color = black);
Cir2 := circle(K, r3, color = black);
Cir3 := circle(J, r1, color = black);
Cir4 := circle(H, r4, color = black);
Poly := polygonplot([p1, p2, p4, p3, p1], color = blue, transparency = 0.9);
SegP1P4 := plot([p1, p4], color = black);
SegP2P3 := plot([p2, p3], color = black);
quadri := plot([p1, p2, p4, p3, p1], color = black);
rect := polygonplot([H, Ii, J, K, H], color = red, transparency = 0.7);
tex := textplot([[p1[], "P1"], [p2[], "P2"], [p3[], "P3"], [p4[], "P4"], [H[], "H"], [Ii[], "I"], [J[], "J"], [K[], "K"]], align = ["above", "right"]);
display(Poly, tex, quadri, rect, SegP1P4, SegP2P3, Cir1, Cir2, Cir3, Cir4, implicitplot([x^2 + y^2 - 1], x = -2 .. 2, y = -4 .. 2, colour = [blue], scaling = constrained), pointplot(Pts, symbolsize = 16), axes = none);
catch:
  NULL;
end try;
end proc:

animate(F,[param1,5,-1/10,-19/2], param1=1/10..1);

Download jamet_anim1.mw

First 101 102 103 104 105 106 107 Last Page 103 of 592