acer

26372 Reputation

29 Badges

16 years, 322 days

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

What are your ranges for alpha and eta?

Read binary file

I see that link in the Recent Questions subpanel (when I scroll down on my phone, or at the right on my computer).

I think that it didn't disappear from the Recent Questions subpanel/tab. But it wasn't listed in the recent/all set of freshly changed postings. I edited it, to make it appear there too.

@Axel Vogt Thanks, that's interesting. Above I used Maple 2022.1.

Those other examples like fsolve(IIf, 7..8) work OK for me (and quickly) in my Maple 2021.2 for Linux. I'll try and find a moment to check Windows.

@wswain You need to assign an initial value to the iterations variable. If you don't then the line,

  ++iterations;

involves a recursive assignment.

Eg,

Trial := proc() 
    
    local y2, test1, iterations:=0; 
    y2 := Vector(5, fill = 0.); 
    y2[5] := 1; 
    test1 := 1.0; 
    
    while 1e-5 < test1 do
        test1 /= 2.0; 
        ++iterations; 
    end do; 
    
    test1; 
    iterations; 
        
    end proc:

Trial();

        17

@Mac Dude Like many quadrature algorithms, Romberg's method's accuracy is related to smooth higher derivatives. It's not well suited to this problem, and accuracy problems here should not be a surprise.

I suspect that the OP asked about numeric quadrature because his choice of syntax for symbolic integration incurs a time-consuming computation. But, as I hoped to have shown, the symbolic integration can produce an accurate result quite quickly, using a more suitable syntax.

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

1 2 3 4 5 6 7 Last Page 1 of 490