Rouben Rostamian

Convert to D...

Answered: Rouben Rostamian 7711

September 24 2022

1 6

mtaylor does not like expressions with "diff" in them, but it works well when converted to "D". So all you need is:

convert(lhs(Eq1), D):
mtaylor(%, [x,t], 2);

Geodesics between two points on a torus...

Answered: Rouben Rostamian 7711

September 21 2022

3 7

I am posting this as an "Answer" since the question answered here is quite distinct from the question that was originally asked and answered in this thread.

Geodesics connecting two points on a torus

>

restart;

>	with(plots):

>	with(VariationalCalculus):

The parametric equation of a torus with major and minor radii of and .

>	< (a+bcos(v))cos(u), (a+bcos(v))sin(u), b*sin(v) >; x := unapply(%, u, v):

Vector(3, {(1) = (a+b*cos(v))*cos(u), (2) = (a+b*cos(v))*sin(u), (3) = b*sin(v)})

Calculate the coefficients , , and of the first fundamental form:

>	xu := diff(x(u,v),u): xv := diff(x(u,v),v): E := unapply(simplify(xu^+ . xu), u, v); F := unapply(simplify(xu^+ . xv), u, v); G := unapply(simplify(xv^+ . xv), u, v);

Look for geodesics of the form . These account for all geodesics with the exception

of those corresponding to constant. The latter are plain circles so we don't need Maple

to calculate them.

From differential geometry, the length of any curve of the form

is given by L = int(Upsilon(u), u = u__1 .. u__2) where:

>	Upsilon := sqrt(E(u, v(u)) + 2F(u, v(u))diff(v(u),u) + G(u,v(u))*(diff(v(u),u))^2);

((a+b*cos(v(u)))^2+b^2*(diff(v(u), u))^2)^(1/2)

The minimum length is obtained with the help from the Euler-Lagrange equation

which leads to a second order nonlinear differential equation in :

>	EulerLagrange(Upsilon, u, v(u)): remove(type, %, equation)[]: isolate(%, diff(v(u), u, u)): de := simplify(%);

diff(diff(v(u), u), u) = -(2*b^2*(diff(v(u), u))^2+(a+b*cos(v(u)))^2)*sin(v(u))/((a+b*cos(v(u)))*b)

Pick numbers for the torus's radii:

>	a, b := 5, 2;

>	torus := plot3d(x(u,v), u=-Pi..Pi, v=-Pi..Pi, scaling=constrained, color="khaki", style=surface);

Example 1

Calculate a geodesic from to , :

>	bc1 := v(0)=0, v(Pi/2)=Pi/2;

>	dsol1 := dsolve({de, bc1}, numeric, output=operator);

>

display(
        torus,
        tubeplot([seq(eval(x(u, v(u)), dsol1))], u=0..Pi/2, color=red, radius=0.1),
        pointplot3d(eval([x(0,v(0)),x(Pi/2,v(Pi/2))], dsol1), symbol=solidsphere, symbolsize=20, color="Orange"),
lightmodel=light4, thickness=5, orientation=[-5,75,0], size=[700,400],
style=surface, axes=none);

>

Example 2

>	bc2 := v(0)=0, v(Pi/2)=Pi;

>	dsol2 := dsolve({de, bc2}, numeric, output=operator);

>

display(
        torus,
        tubeplot([seq(eval(x(u, v(u)), dsol2))], u=0..Pi/2, color=red, radius=0.1),
        pointplot3d(eval([x(0,v(0)),x(Pi/2,v(Pi/2))], dsol2), symbol=solidsphere, symbolsize=20, color="Orange"),
lightmodel=light4, thickness=5, orientation=[-65,65,0], size=[700,400],
style=surface, axes=none);

>

Example 3

The end points in this example are identical to those of Example 1 but the new geodesic,
drawn in cyan, is different. In fact, there are infinitely many geodesics that connect

those two points.

>	bc3 := v(0)=0, v(Pi/2)=5*Pi/2;

>	dsol3 := dsolve({de, bc3}, numeric, output=operator);

>

display(
        torus,
        tubeplot([seq(eval(x(u, v(u)), dsol1))], u=0..Pi/2, color=red, radius=0.1),
        tubeplot([seq(eval(x(u, v(u)), dsol3))], u=0..Pi/2, color=cyan, radius=0.1),
        pointplot3d(eval([x(0,v(0)),x(Pi/2,v(Pi/2))], dsol3), symbol=solidsphere, symbolsize=20, color="Orange"),
lightmodel=light4, thickness=5, orientation=[-65,45,0], size=[700,400],
style=surface, axes=none, viewpoint=circleright);

>

Download geodesics-between-two-points.mw

Finding geodesics...

Answered: Rouben Rostamian 7711

September 16 2022

6 5

You will find the geodesic equations in the wikipedia web page that you have referenced to. Admittedly it is rather cryptic, involving Christoffel symbols in an implicit summation. Once we unpack the notation, the geodesic equations for a surface turn out to be a pair of coupled, second order, nonlinear ODEs. These equations are usually derived in an introductory course on differential geometry. To get to those equations, one needs to learn about a surface's first fundamental form and how the Christoffel symbols are expressed in terms of the coefficients of the first fundamental form.

The Christoffel() proc in the attached worksheet calculates the coefficients of the first fundamental form of any surface and then calculates and returns the corresponding Christoffel symbols. To understand the details, you need to know some differential geometry.

A second proc, Geodesic_DEs(), applies the equation in the Wikipedia page to calculate and return the surface's geodesic equations. Once you have those equations, you can apply Maple's dsolve() to calculate and plot the geodesics. The worksheet includes an illustrative example.

Maple's Differential Geometry package has a GeodesicEquations() function. Odds are that what I have calculated here can be done more quickly by applying that function but don't know how. Perhaps someone who knows may show us.

Answers to the several questions that you asked would involve some elementary differential geometry, so I cannot summarize them here in a few lines. I can, however, give short answers to two of those questions:

What does "to walk forward" mean? When you move on a surface, at each location you have a velocity vector, an acceleration vector, and a normal vector (i.e., vector perpendicular to the surface). You are "moving forward", or equivalently, "traveling on a geodesic", provided that those three vectors are co-planar.
Geodesic is the shortest distance between two points in space, but no such end points are shown in the website's animation. A geodesic typically extends indefinitely in both directions. There are no end points. A curve is geodesic if for any pair of points A, B on the curve, the part of the curve that lies between A and B is the shortest path among all curves that lie on the surface and connect A to B.

>

restart;

>	with(plots):

Receives a parametrization of a surface and returns the corresponding

Christoffel symbols , each of which is a procedure.

>

Christoffel := proc(x::procedure)
        local u, v, xu, xv, E, F, G, Eu, Ev, Fu, Fv, Gu, Gv, den, C;
        xu := diff(x(u,v),u);
        xv := diff(x(u,v),v);
        # E, F, G are the coefficients of the surface's first fundamental form
        E := xu^+ . xu;
        F := xu^+ . xv;
        G := xv^+ . xv;
        Eu := diff(E, u);
        Ev := diff(E, v);
        Fu := diff(F, u);
        Fv := diff(F, v);
        Gu := diff(G, u);
        Gv := diff(G, v);
        den := 2*(E*G - F^2);
        C[1][1,1] := unapply(simplify((G*Eu - 2*F*Fu + F*Ev)/den), u, v);
        C[1][1,2] := unapply(simplify((G*Ev - F*Gu)/den), u, v);
        C[1][2,2] := unapply(simplify(-(F*Gv - 2*G*Fv + G*Gu)/den), u, v);
        C[2][1,1] := unapply(simplify(-(E*Ev - 2*E*Fu + F*Eu)/den), u, v);
        C[2][1,2] := unapply(simplify((E*Gu - F*Ev)/den), u, v);
        C[2][2,2] := unapply(simplify((E*Gv - 2*F*Fv + F*Gu)/den), u, v);
        C[1][2,1] := C[1][1,2];
        C[2][2,1] := C[2][1,2];
        return C;
end proc:

Example: A sphere

>	x := (u,v) -> < cos(v)cos(u), cos(v)sin(u), sin(v) >;

Here are the 8 Christoffel symbols of our parametrization of the sphere:

>	Gamma := Christoffel(x):

>	Gamma[1][1,1](u,v), Gamma[1][1,2](u,v), Gamma[1][2,1](u,v), Gamma[1][2,2](u,v);

>	Gamma[2][1,1](u,v), Gamma[2][1,2](u,v), Gamma[2][2,1](u,v), Gamma[2][2,2](u,v);

Example: A torus of major radius and minor radius

>	x := (u,v) -> < (a+bcos(v))cos(u), (a+bcos(v))sin(u), b*sin(v) >;

>	Gamma := Christoffel(x):

>	Gamma[1][1,1](u,v), Gamma[1][1,2](u,v), Gamma[1][2,1](u,v), Gamma[1][2,2](u,v);

>	Gamma[2][1,1](u,v), Gamma[2][1,2](u,v), Gamma[2][2,1](u,v), Gamma[2][2,2](u,v);

The geodesic equations

>

Geodesic_DEs := proc(x::procedure)
        local Gamma, de, i, j, k, X;
        Gamma := Christoffel(x);
        for i from 1 to 2 do
           de[i] := diff(X[i](t), t, t) + add(add(
                        Gamma[i][j,k](X[1](t),X[2](t))*diff(X[j](t),t)*diff(X[k](t),t),
                        j=1..2), k=1..2) = 0;
        end do:
        eval([de[1],de[2]], {X[1]=u, X[2]=v})[];   # change notation from X to u, v
        return (%);
end proc:

>

Geodesics on a torus

Here we calculate and plot geodesics on a torus. You may change the torus

to any other parametrized surface. Everthing else remains the same.

>	< (a+bcos(v))cos(u), (a+bcos(v))sin(u), b*sin(v) >; eval(%, {a=5, b=2}); x := unapply(%, u, v):

Vector(3, {(1) = (a+b*cos(v))*cos(u), (2) = (a+b*cos(v))*sin(u), (3) = b*sin(v)})

Vector[column](%id = 36893628224392653508)

>	the_torus := plot3d(x(u,v), u=-Pi..Pi, v=-Pi..Pi, scaling=constrained, style=wireframe);

The differential equations of geodesics on torus

>	DEs := Geodesic_DEs(x);

diff(diff(u(t), t), t)-4*sin(v(t))*(diff(v(t), t))*(diff(u(t), t))/(5+2*cos(v(t))) = 0, diff(diff(v(t), t), t)+(cos(v(t))*sin(v(t))+(5/2)*sin(v(t)))*(diff(u(t), t))^2 = 0

Pick any desired starting point and direction on the torus. The angle determines
the geodesic's initial direction.

>	unassign('theta'); ICs := u(0) = 0, v(0) = 0, D(u)(0) = cos(theta), D(v)(0) = sin(theta);

Pick a and solve the geodesic DEs

>	theta:= 0.95*Pi/2: dsol := dsolve({DEs, ICs}, numeric, output=operator);

Plot the geodesic:

>	eval(x(u(t),v(t)), dsol); the_curve := spacecurve(%, t=0..15, color="Red", thickness=5);

Vector(3, {(1) = (5+2*cos(v(t)))*cos(u(t)), (2) = (5+2*cos(v(t)))*sin(u(t)), (3) = 2*sin(v(t))})

Here is a composite drawing of the torus and the geodesic

>	display(the_torus, the_curve);

Download geodesics.mw

Here it is...

Answered: Rouben Rostamian 7711

September 10 2022

2 0

I did this demo for my students some time ago.

Download mw.mwmw.mw

PS: It's also possible to plot the springs as coils, as here and here but I prefer the wavy springs shown above for this particular example.

Patches on a sphere...

Answered: Rouben Rostamian 7711

September 03 2022

4 3

This worksheet will work provided that

the north pole is inside the patch;
the patch is "star-shaped" (that's a technical term) with respect to the pole, that is, every meridian (with 0 <= theta <= Pi) intersects the patch's boundary only once .

These conditions are satisfied in your case.

>

restart;

>	with(plots):

Tranlsates cartesian coordinates (given as a vector) to spherical coordinates (returned as a list)

>	cartesian_to_spherical := proc(v::Vector(3)) local x, y, z, r, theta, phi; x, y, z := seq(v); r := sqrt(x^2+y^2+z^2); theta := arccos(z/r); phi := arctan(y/r, x/r); return simplify([r, theta, phi]); end proc:

>	S := (theta, phi) -> < sin(theta)cos(phi), sin(theta)sin(phi), cos(theta) >;

(1)

Shrink the plotted sphere just a bit so that the superimposed patches stand out:

>	sphere_plot := plot3d(0.999*S(theta,phi), theta=0..Pi, phi=-Pi..Pi, color="Orange", style=patch);

C1 is given as a parametric curve in Cartesian coordinates. The name of the parameter, ,
is misleading since it is unrelated to the spherical coordinate of the same name.

I changed the parameter name to to reduce confusion.

>

C1 := u -> <
        0.3071867916 - 0.3748494088*sin(u) + 0.3092507622*cos(u),
        0.3858134286*cos(u) - 0.5068582065 + 0.1636706976*sin(u),
        0.3224530297*sin(u) + 0.6143735837 + 0.1636706976*cos(u)>:
umin, umax := 0.8870700404, 1.362663780:
# convert the curve to spherical coordinates
cartesian_to_spherical(C1(u)):
my_theta, my_phi := %[2], %[3]:
curve_plot1 := spacecurve(C1(u), u=umin..umax, color="Cyan", thickness=3):
patch_plot1 := plot3d(S(theta,my_phi), theta=0..my_theta, u=umin..umax, color="Red"):

Repeat for curves 2 and 3:

>

C2 := u -> <
        0.5464684013 - 0.3515280247*sin(u),
        0.4434990987*cos(u),
        0.2704061730*sin(u) + 0.7104089214 >:
umin, umax := 0.7339, 1.8831:
# convert the curve to spherical coordinates
simplify(cartesian_to_spherical(C2(u))):
my_theta, my_phi := %[2], %[3]:
curve_plot2 := spacecurve(C2(u), u=umin..umax, color="Cyan", thickness=3):
patch_plot2 := plot3d(S(theta,my_phi), theta=0..my_theta, u=umin..umax, color="Green"):

>

C3 := u -> <
        0.2150819669 - 0.3346119696*sin(u),
        0.3437498130*cos(u),
        0.07873222801*sin(u) + 0.9140983608 >:
umin, umax := -0.2908, 2.4602:
# convert the curve to spherical coordinates
cartesian_to_spherical(C3(u)):
my_theta, my_phi := %[2], %[3]:
curve_plot3 := spacecurve(C3(u), u=umin..umax, color="Cyan", thickness=3):
patch_plot3 := plot3d(S(theta,my_phi), theta=0..my_theta, u=umin..umax, color="Blue"):

>	display(sphere_plot, curve_plot1, patch_plot1, curve_plot2, patch_plot2, curve_plot3, patch_plot3, style=surface);

Download Patch_on_a_Sphere_rr.mw

Adjustments...

Answered: Rouben Rostamian 7711

August 25 2022

5 4

The equations of motion are of the 2nd order in each of the variables Ext, theta, and phi. So you will need six initial conditions but you have provided only five. You need to specify an initial value for D(Ext).

That alone doesn't help, however. There must be some error in your derivation of the equations of motion. I didn''t have the patience to examine the ad hoc details of your calculations. I just re-did the derivation in a systematic and generic way.

In the attached worksheet I have changed parameter values and initial conditions to achieve a good-looking result. You may restore them to your original selections if you wish.

>

restart;

>	with(VariationalCalculus):

>	with(plots):

Bob's position vector

>	R := (e, theta, phi) -> (a + e) * < sin(theta)cos(phi), sin(theta)sin(phi), cos(theta) >;

Bob's position vector at time t:

>	r := R(e(t), theta(t), phi(t));

Vector(3, {(1) = (a+e(t))*sin(theta(t))*cos(phi(t)), (2) = (a+e(t))*sin(theta(t))*sin(phi(t)), (3) = (a+e(t))*cos(theta(t))})

Bob's velocity vector

>	v := diff(r, t);

The kinetic and potential energies:

>	KE := 1/2m simplify(v^+ . v); PE := -mgr[3] + 1/2ke(t)^2;

(1/2)*m*(-(cos(theta(t))-1)*(cos(theta(t))+1)*(a+e(t))^2*(diff(phi(t), t))^2+(a+e(t))^2*(diff(theta(t), t))^2+(diff(e(t), t))^2)

The Lagrangian:

>	L := KE - PE;

(1/2)*m*(-(cos(theta(t))-1)*(cos(theta(t))+1)*(a+e(t))^2*(diff(phi(t), t))^2+(a+e(t))^2*(diff(theta(t), t))^2+(diff(e(t), t))^2)+m*g*(a+e(t))*cos(theta(t))-(1/2)*k*e(t)^2

The Euler-Lagrange equations:

>	EulerLagrange(L, t, [e(t), theta(t), phi(t)]): EL := remove(type, %, equation);

${(1/2)*m*(-2*(cos(theta(t))-1)*(cos(theta(t))+1)*(a+e(t))*(diff(phi(t), t))^2+2*(a+e(t))*(diff(theta(t), t))^2)+m*g*cos(theta(t))-k*e(t)-m*(diff(diff(e(t), t), t)), (1/2)*m*(sin(theta(t))*(cos(theta(t))+1)*(a+e(t))^2*(diff(phi(t), t))^2+(cos(theta(t))-1)*sin(theta(t))*(a+e(t))^2*(diff(phi(t), t))^2)-m*g*(a+e(t))*sin(theta(t))-2*m*(a+e(t))*(diff(theta(t), t))*(diff(e(t), t))-m*(a+e(t))^2*(diff(diff(theta(t), t), t)), -m*(diff(theta(t), t))*sin(theta(t))*(cos(theta(t))+1)*(a+e(t))^2*(diff(phi(t), t))-m*(cos(theta(t))-1)*(diff(theta(t), t))*sin(theta(t))*(a+e(t))^2*(diff(phi(t), t))+2*m*(cos(theta(t))-1)*(cos(theta(t))+1)*(a+e(t))*(diff(phi(t), t))*(diff(e(t), t))+m*(cos(theta(t))-1)*(cos(theta(t))+1)*(a+e(t))^2*(diff(diff(phi(t), t), t))}$

>	ics := { e(0) = 0, D(e)(0) = 0, theta(0) = Pi/3, D(theta)(0) = 0, phi(0) = Pi/6, D(phi)(0) = 1 };

>	g := 1; m := 1; k := 5; a := 1; tmax := 30;

>	dsol := dsolve(EL union ics, numeric, output=operator);

>	odeplot(dsol, [t, e(t)], t=0..tmax);

>	orbit := spacecurve(eval([r[1],r[2],-r[3]], dsol), t=0..tmax);

>

frame := proc(t)
        local r;
        eval(R(e(t), theta(t), phi(t)), dsol);
        r := convert(%, list);
        display(
                plottools:-line([0,0,0], r, thickness=3, color=red),
                pointplot3d([0,0,0], symbol=solidbox, symbolsize=20, color=black),
                pointplot3d(r, symbol=solidsphere, symbolsize=40, color=red));
        # reflect about the horizontal plane so that the z
        # axis points upward
        plottools:-reflect(%, [[0,0,0], [1,0,0], [0,1,0]]);
end proc:

>	animate(frame, [t], t=0..tmax, background=orbit, scaling=constrained, frames=100);

Download pendulum-elastic-3d.mw

Try this...

Answered: Rouben Rostamian 7711

August 24 2022

3 2

You didn't provide samples of your expressions, so I made up one. See if this does what you want.

Download mw.mw

Calculating the flux...

Answered: Rouben Rostamian 7711

August 23 2022

2 6

You have asked how to calculate the partial derivative of u with respect to y. Here is how.

>

restart;

>

N:=10;

>	PDE:=diff(u(y, t), t) = diff(u(y, t), y, y);

>	ICBC:= {u(y,0) = 0, u(0,t) = cos(t), u(N, t) = 0};

>	dy := N/100; # space step dt := dy^2; # time step

>	sol1:=pdsolve(PDE,ICBC,numeric,spacestep=dy, timestep=dt):

>	sol1:-plot3d(t=0..Pi);

>	vals := sol1:-value(u(y,t), output=listprocedure);

This proc captures the value of the solution :

>	U := eval(u(y,t), vals);

Here is the solution evaluated at :

>

U(1,0.5);

Here is the flux, , evaluated at :

>	-D[1](U)(1,0.5);

The plot of the flux (it takes some time, be patient):

>	plot3d(-D[1](U)(y,t), y=0..N, t=0..Pi);

Download numerical-flux.mw

Maple worksheet...

Answered: Rouben Rostamian 7711

August 15 2022

3 7

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

Automatically loading the Units[Simple] subpackage

>	N__bld := 4;

>	D__rtr := 30Unit('ft'); R__tip := D__rtr/2; chd := 8Unit('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.1Unit(1/'deg'); C__L := alpha -> alphaCl__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);

proc (x) options operator, arrow; 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 end proc

>

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/(omegar))Unit('radian') - 8Unit('deg')r/R__tip + 6*Unit('deg');

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

>	dr__s := (r, theta, u) -> 1/2rho__airchd(omega^2r^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 -> 2rho__airu^2PiR__tip^2;

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

>

Questions

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

>	eqn2(6Unit('degree'), 5Unit('m'/'sec'));

>	eqn2(6Unit('degree'), 10Unit('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'(6Unit('degree'), uUnit('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(thetaUnit('degree'), uUnit('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( thetaUnit('degree'), uUnit('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);

>

Download rotor-equation.mw

Implementing recursion...

Answered: Rouben Rostamian 7711

August 13 2022

1 1

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

Download newton-iteration.mw

Perhaps this is what you want?...

Answered: Rouben Rostamian 7711

August 12 2022

0 1

>

restart;

>	ftaylor := mtaylor(g(x,y), [x, y], 3);

g(0, 0)+(D[1](g))(0, 0)*x+(D[2](g))(0, 0)*y+(1/2)*(D[1, 1](g))(0, 0)*x^2+(D[1, 2](g))(0, 0)*x*y+(1/2)*(D[2, 2](g))(0, 0)*y^2

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

Download mw.mw

Try this...

Answered: Rouben Rostamian 7711

July 22 2022

1 1

Shorter code...

Answered: Rouben Rostamian 7711

July 17 2022

2 0

TomLeslie has already answered your question. Here is a shorter alternative.

>

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

Matrix(3, 3, {(1, 1) = 1, (1, 2) = 2, (1, 3) = 3, (2, 1) = 3, (2, 2) = 2, (2, 3) = 1, (3, 1) = 4, (3, 2) = 6, (3, 3) = 5})

>	uniqueRows(A);

>	A := Matrix([[1, 2, 3], [1,2,3], [4, 6, 5]]);

Matrix(3, 3, {(1, 1) = 1, (1, 2) = 2, (1, 3) = 3, (2, 1) = 1, (2, 2) = 2, (2, 3) = 3, (3, 1) = 4, (3, 2) = 6, (3, 3) = 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]]);

Matrix(3, 3, {(1, 1) = 1, (1, 2) = 2, (1, 3) = 3, (2, 1) = 3, (2, 2) = 2, (2, 3) = 1, (3, 1) = 4, (3, 2) = 6, (3, 3) = 5})

>	B := Matrix([[3, 2, 1], [1, 2, 3], [4, 6, 5]]);

Matrix(3, 3, {(1, 1) = 3, (1, 2) = 2, (1, 3) = 1, (2, 1) = 1, (2, 2) = 2, (2, 3) = 3, (3, 1) = 4, (3, 2) = 6, (3, 3) = 5})

>	RowEqual(A, B);

>	B := Matrix([[5, 2, 1], [1, 2, 3], [4, 6, 5]]);

Matrix(3, 3, {(1, 1) = 5, (1, 2) = 2, (1, 3) = 1, (2, 1) = 1, (2, 2) = 2, (2, 3) = 3, (3, 1) = 4, (3, 2) = 6, (3, 3) = 5})

>	RowEqual(A, B);

Download mw.mw

Not a solution...

Answered: Rouben Rostamian 7711

July 08 2022

2 2

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.

Square both sides of the equation to get (dy/dx)^2 = 1 + x + y;
Rearrange into y + 1 = - x + (dy/dx)^2
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...

Answered: Rouben Rostamian 7711

July 05 2022

5 0

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

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

7711 Reputation

22 Badges

MaplePrimes Activity

These are answers submitted by Rouben Rostamian

Convert to D...

Geodesics between two points on a torus...

Finding geodesics...

Here it is...

Patches on a sphere...

Adjustments...

Try this...

Calculating the flux...

Maple worksheet...

Implementing recursion...

Perhaps this is what you want?...

Try this...

Shorter code...

Not a solution...

Maybe this...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again