## 6554 Reputation

16 years, 53 days

## Maple worksheet...

Disclaimer: This is not an answer!

I am posting this as an "Answer" to set it apart from the twisty discussion thread that we have had so far.

I have converted the OP's equations and data into a self-contained Maple worksheet.  At the end of the worksheet I have inserted the questions that the OP is seeking answers for.  My attempts to solve the equations were unsuccessful.  See what you can do with this.

 > restart;
 > kernelopts('version'); > 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; > 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; > eqn2 := (theta, u) -> Thust(theta, u) - UD(u); >

Questions

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

 > eqn2(6*Unit('degree'), 5*Unit('m'/'sec')); > eqn2(6*Unit('degree'), 10*Unit('m'/'sec')); 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);

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);

Observation: The following does work!
I plot the 3D graph of z = eqn2(theta,u) along with the plane z=0.

The intersection curve of those two surfaces is the graph of the function U defined above.

Note: I don't know how to make surfdata work with units,
so I divide the value of  Thust by Unit(N) to remove the unit.

 > [seq([seq(         [theta, u, Thust( theta*Unit('degree'), u*Unit('ft'/'sec') )/Unit(N)],                         theta=0..15, 3)], u = 0..100, 10)]: plots:-surfdata(%): plot3d(0, theta=0..15, u=0..100): plots:-display(%,%%, style=surface); >

## Implementing recursion...

What you have shown appears to be an incomplete fragment (what is "Thust", for instance?) and therefore I cannot tell you how to fix that.  If Maple's fsolve cannot solve your equation, I doubt that your own Newton's root-finder will do any better.

Nonetheless, since you asked for a recursive implementation of Newton's root-finding method, here it is, for whatever it's worth.

 > restart;

Implements the Newton root-finding iteration.

Required arguments:

f:        a function

x0:        starting point of the root-finding iteration

Optional arguments:
eps: tolerance -- iteration stops when consecutive x values are less than eps
default eps = 1e-5
nmax: iteration is aborted if numer of iterations exceeds nmax
default nmax = 10

 > NewtIt := proc(f, x0, {eps:=1e-5, nmax:=10})         local x, xnew;         x := evalf(x0);         xnew := x - f(x)/D(f)(x);         if is(abs(xnew - x) < eps) then                 return xnew;         elif nmax = 0 then                 error "number of iterations exceeded nmax";         else                 return `procname`(f, xnew, :-eps=eps, :-nmax=nmax-1);         end if; end proc:
 > NewtIt(x->x^2 - 4, 10); > NewtIt(x->x^2 - 4, 10, eps=0.01); > NewtIt(x->x^2 - 4, 10, eps=0.01, nmax=2);

Error, (in NewtIt) number of iterations exceeded nmax

## Perhaps this is what you want?...

 > restart;
 > ftaylor := mtaylor(g(x,y), [x, y], 3); > g := (x,y) -> exp(-x^2-y^2); > ftaylor; Or maybe

 > restart;
 > ftaylor := n -> mtaylor(g(x,y), [x, y], n); > g := (x,y) -> exp(-x^2-y^2); > ftaylor(3); > ftaylor(5); ## Try this...

 > restart;
 > to_list := u -> `if`( type(u,`*`), convert(u, list), [u]):

 > to_list(a^2*b^5); > to_list(a^2); ## Shorter code...

 > restart;

Are the rows of the matrix unique?

 > uniqueRows := proc(A::Matrix)         local rows, n1, n2;         rows := convert(A, listlist);         n1 := nops(rows);         rows := convert(rows, set);         n2 := nops(rows);         return evalb(n1=n2); end proc:
 > A := Matrix([[1, 2, 3], [3, 2, 1], [4, 6, 5]]); > uniqueRows(A); > A := Matrix([[1, 2, 3], [1,2,3], [4, 6, 5]]); > uniqueRows(A); Determine if A and B are the same modulo row permutations

 > RowEqual := proc(A::Matrix, B::Matrix)         local rowsA, rowsB;         rowsA := sort(convert(A, listlist));         rowsB := sort(convert(B, listlist));         return ArrayTools:-IsEqual(rowsA, rowsB); end proc:
 > A := Matrix([[1, 2, 3], [3, 2, 1], [4, 6, 5]]); > B := Matrix([[3, 2, 1], [1, 2, 3], [4, 6, 5]]); > RowEqual(A, B); > B := Matrix([[5, 2, 1], [1, 2, 3], [4, 6, 5]]); > RowEqual(A, B); ## Not a solution...

I don't know why you refer to singular solutions. I don't see anything "singular" there.

Neither of the two solutions produced by your construction is a solutions.  You are "misleading" Maple into looking at the ODE as one of the d'Alembert type.  Maple blindly obliges and produces the wrong results.

Here is how one (illegally) casts your equation into a d'Alembert type.

1. Square both sides of the equation to get (dy/dx)^2 = 1 + x + y;
2. Rearrange into y + 1 = - x + (dy/dx)^2
3. Let z = 1 + y.  Then z = - x + (dz/dx)^2

This is a d'Alembert ODE (see ?Solving d'Alembert ODEs) with f(s)=-1 and g(s)=s^2.

Where is the error?  It is in Step 1.  By squaring the two sides of the ODE, we are introducing spurious solutions.  Maple's solution would be correct (for the wrong reasons) if you replace the right-hand side of your ODE by its negative.

To avoid the issue, don't impose the dalembert requirement.  Do dsolve(ode_1) to let Maple do its own thing and produce the correct, albeit implicit, solution.

## Maybe this...

This comes close to what you have sketched.

```p := plots:-display(
plot(floor(x), x=1..8, discont=[color="Red"], color="Green"),
plot(floor(x), x=-8..-1/2, discont=[color="Red"], color="Green"),
symbol=solidcircle, symbolsize=10, thickness=3):
plottools:-transform((x,y)->[1/x, 1/y])(p);``` ## Here is how...

The array d contains the desired coordinates, one point per line.
 > restart;
 > c := plots:-spacecurve([cos(t), sin(t), t], t=0..2*Pi); > plottools:-getdata(c): d := op(3,%); ## Your f,fdiff=... command is equivalent t...

Your f,fdiff=... command is equivalent to the simpler:

`f, diffs := GenerateMatrix(sys222, var1);`

However...

A matrix representation of a system of equations is useful if the system is linear. Your system is nonlinear. What do you expect to get out of a matrix representation of a nonlinear system?

## The "symbolic" option...

Regarding the symbolic option to simplify(), the help page says:

The result of such an operation is in general not valid over the whole complex plane and can lead to incorrect results if you assume the expressions represent analytical functions.

Thus, the symbolic option should be used only when you really know what you are doing.  Preferably it should not be used at all.  In your worksheet you are comparing two simplifications, one with and one without the symbolic option.  It's not surprising that the result are inconsistent.

## Probably this is what you want...

```restart;
doit := proc(oper)
local x := 3, y := 4;
return oper(x,y);
end proc:
doit(`+`);
7
doit(`-`);
-1
doit(`*`);
12
doit(`/`);
3/4
```

## The error message says it all...

The error message says it all; it doesn't know what to do with d.

In the IBC you have

`C(0, t) = 1 + d(1 - cos(varpi*t))`

What is d?  If it's a function, then you need to specify it.  If it's a multiplier constant, then you need to provide its value and also don't forget to insert a multiplication sign between d and the adjacent parenthesis.

Regarding the differential equation for theta:

It has A and A in it, but there are no A and A in your expected solution.  Is that an oversight?

By the way, that differential equation is quite trivial to solve by hand.  It's not worth bothering with Maple with it. If you try the manual solution, you will find that a solution exists provided that  - v  - A + A = v.

Regarding the differential equation for phi:

Express the differential equation in terms of Heaviside() rather than pieceswise(), as in:

```eq2 := diff(phi(t),t) = -v + v*Heaviside(t-T) + v*Heaviside(t-T);
dsolve({eq2,phi(0)=-v});```

That yields the unique solution of that initial value problem, and therefore determines the value of phi(1).  You cannot arbitrarily demand a different value for phi(1).

## The semi-cubic path...

As dharr has correctly pointed out, a nonzero initial velocity obviates the objection that I had expressed in my prevrious message.  Following through with that insight, we end up with a semi-cubic parabolic curve y^3 = - a*x^2 for the particle's path.  The adjustment, relative to dharr's solution, is that the x compoment of the initial vecocity should be zero rather than v0.

 > restart;

The sum of kinetic and potential energies:

 > E := 1/2*m*(D(x)(t)^2 + D(y)(t)^2) + m*g*y(t) (1)

Conservation of energy

 > eq := E = eval(E, t=0); (2)

Suppose that the particle starts off at the origin, that is, , and slides down a curve (to be determined)
under the gravitational pull.  Assume that the vertical component
of the velocity is a constant for some .  Therefore .
Additionally, assume that the particle begins falling vertically,

that is, Plugging this information into the equation of

conservation of energy we obtain a differential equation for :

 > de := eval(eq, {x(0)=0, y(0)=0, D(x)(0)=0, y = (t -> -V*t)}); (3)

Solve the differential equation

 > dsol := dsolve({de, x(0)=0}); (4)

Therefore the particle's path is given by the parametric curve:

 > curve := x = rhs(dsol), y = -V*t; (5)

Eliminate to obtain the equation of the path as a relationship

between and  . We end up with a semi-cubic of the form :

 > x = subs(t=-y/V, rhs(dsol)): %^2: isolate(%, y^3); (6)

Here is the graph of the path, with a certain choice of parameters:

 > subs(curve, g=1, V=1, [x,y,t=0..1]); plot(%);  ## Do it with ImageMagick...

I don't think you can do that in Maple.  But you can do it through post-processing with open and free software ImageMagick assuming that you are comfortable with using command-line tools.

Save your Maple graphics in a file, let's say p1.png (or GIF), and then convert it to p2.png:

`convert p1.png -transparent white p2.png`

Beware that this converts all while pixels, (not just the background!) to transparent.

 1 2 3 4 5 6 7 Last Page 1 of 47
﻿