acer

26592 Reputation

29 Badges

17 years, 0 days

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@Carl Love I thought those produce 1x3 Matrices. The OP asked for row Vectors, so I rejected those ways.

(I too concocted those ways, at first. And I thought that I'd checked them, in Maple 2022.1, ... but perhaps I missed something.)

@MalakMMK Did you look where (and how) I used either a simple iterative formula or a few iterated steps?

That is, using u[i], v[i], and u[i+1] instead of, say, x[n], y[n], and x[n+1].

Have you tried adjusting the formulas yourself?

Is this work for an academic course or project?

@MalakMMK I don't understand what you're trying to say.

@mmcdara Sure. Tabulate is just an easy-to-use front end, for some common situations.

Tabulate is built from those lower level DocumentTools widgets&things (Layout & Components), of which you've shown an example. Naturally, one can build more flexible or targeted applications from those, for specific situations.

@tomleslie The OP stated, "Note I want to monitor points during the loop running not after it finishes.", and what you've shown doesn't provide that.

@vs140580 Please add your very close followup query here, and not in a separate Question thread.

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?

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