Rouben Rostamian

Numerical solution...

Answered: Rouben Rostamian 8259

May 11 2023

2 1

The second form of your differential equation has explicitly defined numerical coefficients except for V which is left undefined. So we apply dsolve() to solve the differential equation numerically while leaving V as a parameters that can be assigned values retroactively.

>

restart;

>	de := diff(f(t), t) = (.21.2894606222-.21(1.158611162/((1+0.4836862955e-1t)(1/(1+0.4836862955e-1t))^.1986136226))-.2894606222f(t))/V;

diff(f(t), t) = (0.6078673066e-1-.2433083440/((1+0.4836862955e-1*t)*(1/(1+0.4836862955e-1*t))^.1986136226)-.2894606222*f(t))/V

>	ic := f(0) = .21;

>	dsol := dsolve({de,ic}, numeric, parameters=[V], output=operator);

>	F := eval(f, dsol);

>	dsol(parameters=[1]); # set V = 1

>

F(4);

>	plot(F(t), t=0..20);

>	dsol(parameters=[2]); # set V = 2

>

F(4);

>	plot(F(t), t=0..20);

>

Download mw.mw

Try this...

Answered: Rouben Rostamian 8259

April 21 2023

1 0

for i from 1 to 24 do
	print(i);
	if irem(i,5) = 0 then
		print("press enter to continue");
		readline();
	end if;
end do;

Correction...

Answered: Rouben Rostamian 8259

April 16 2023

3 0

That's because in the definition of F you have

z=-2*cos(q)-2.5 .. -2.5-2*cos(q)

That should be

z=-2*cos(q)-2.5 ..  2.5-2*cos(q)

Simplified version...

Answered: Rouben Rostamian 8259

April 09 2023

3 1

C_R has pointed out the source of the error. You are attempting to calculate "le centre" of a paraboloid. But a paraboloid has no "centre".

It seems to me that you are making things much more complicated than need be, and plotting planes and quadratic surfaces through implicitplot3d is a poor choice. These can be much better plotted as parametric surfaces through plot3d. Here, for instance, is the paraboloid you wish to plot:

L := 4:  K := 1.3*L:
display(
	plot3d(< 2*s*cos(t), 3*s*sin(t), s^2 >,
		s=0..2, t=-Pi..Pi, color="Orange"),
	plot3d([s,t,0], s=-L..L, t=-L..L, transparency=0.6, color="DarkGray"),
	plot3d([0,s,t], s=-L..L, t=-L..L, transparency=0.6, color="DarkGray"),
	plot3d([t,0,s], s=-L..L, t=-L..L, transparency=0.6, color="DarkGray"),
	arrow([0,0,0], [K,0,0], color="Red"),
	arrow([0,0,0], [0,K,0], color="Green"),
	arrow([0,0,0], [0,0,K], color="Blue")
);

Download worksheet: quadratics.mw

Multispan beam...

Answered: Rouben Rostamian 8259

April 01 2023

3 0

This website is broken again and won't allow showing the contents of a worksheet.

This worksheet should answer your question.

Download multispan-beam.mw

Streamlines...

Answered: Rouben Rostamian 8259

March 27 2023

2 0

Did you really mean that velocity field? That velocity field is time-dependent. Plotting time-dependent streamlines calls for an animation; a static plot won't do. Usually when people talk about streamlines, they mean the streamlines of a stationary (that is, time-independent) velocity field.

Here I plot the streamlines of a stationary velocity field. This may be adapted to produce an animation of the time-dependent case if necessary, but that would call for a closer familiarity with Maple.

>

restart;

>	with(DEtools):

Velocity field:

>	v := (x,y) -> <1+x^2+y^2, x+y>;

Streamlines are the solutions of this system of differential equations

>	de := [seq](diff(<x(t),y(t)>,t) =~ v(x(t),y(t)));

Here is the corresponding vector field:

>	DEplot(de, [x(t),y(t)], t=0..1, x=-3..3, y=-3..3);

Specify initial conditions along the left edge:

>	ic := seq([x(0)=-3, y(0)=h], h=-3..3, 0.2);

Plot streamlines

>	DEplot(de, [x(t),y(t)], t=0..3, [ic], x=-3..3, y=-3..3, linecolor=black, thickness=1);

Download mw.mw

Rotating helix...

Answered: Rouben Rostamian 8259

March 26 2023

2 3

This worksheet provides an alternative approach to what I had provided earlier in RevolSurfaceQ.mw and adds animation as you have requested.

>

restart;

>	with(plots):

rotates the vector about the vector through an angle ,

The length of the vector is immaterial as long as it is nonzero.

>

R := proc(u::Vector(3))
        local pname, n, phi, nlen, t;
        pname := op(procname);
        if not type(procname, 'indexed') or nops([pname])<>2 then
                error "Two indices needed";
        end if;
        n := pname[1];
        phi := pname[2];
        if not type(n, 'Vector(3)') then
                error "expected Vector(3) for the first argument of R[n,phi]";
        end if;
        nlen := simplify(sqrt(n^+ . n));
        if nlen = 0 then
                error "the vector n in R[n,phi] should have a nonzero length";
        end if;
        n /= nlen;
        t := (n^+ . u)*n;
        t + (u - t)*cos(phi) + LinearAlgebra:-CrossProduct(n,u)*sin(phi);
end proc:

Example

Let's rotate the helix

>	H := < cos(s), sin(s), s/4 >;

Vector(3, {(1) = cos(s), (2) = sin(s), (3) = (1/4)*s})

about the line

>	L := t -> <t,0,t>;

that is, about the vector

>	n := <1, 0, 1>;

Vector(3, {(1) = 1, (2) = 0, (3) = 1})

through an angle . We get

>	S := R[<1,0,1>,q](H);

Vector(3, {(1) = ((1/4)*sqrt(2)*cos(s)+(1/16)*sqrt(2)*s)*sqrt(2)+cos(q)*(cos(s)-((1/4)*sqrt(2)*cos(s)+(1/16)*sqrt(2)*s)*sqrt(2))-(1/2)*sin(q)*sqrt(2)*sin(s), (2) = cos(q)*sin(s)+sin(q)*(-(1/8)*sqrt(2)*s+(1/2)*sqrt(2)*cos(s)), (3) = ((1/4)*sqrt(2)*cos(s)+(1/16)*sqrt(2)*s)*sqrt(2)+cos(q)*((1/4)*s-((1/4)*sqrt(2)*cos(s)+(1/16)*sqrt(2)*s)*sqrt(2))+(1/2)*sin(q)*sqrt(2)*sin(s)})

For a fixed value of and varying , the vector calculated above
represents the parametric equation of the rotated helix.

If both and vary, the vector represents the parametric equation

of the surface swept by the helix as it rotates about the vector .

Here is what the rotating helix looks like:

>	demo1 := q -> display( tubeplot([seq](R[n,q](H)), s=-4Pi..4Pi, radius=0.08, numpoints=200, color="Green"), tubeplot([seq](L(s)), s=-2..2.5, radius=0.08, color="Red")):

>	animate(demo1, [q], q=0..2*Pi, lightmodel=light4, style=surface, orientation=[-120,75,0], labels=[x,y,z]);

and here is the rotating helix and the surface it sweeps:

>	S := a -> display( plot3d(R[n,q](H), s=-4Pi..4Pi, q=0..a, color="Orange"), tubeplot([seq](R[n,a](H)), s=-4Pi..4Pi, radius=0.08, numpoints=200, color="Green"), tubeplot([seq](L(s)), s=-2..2.5, radius=0.08, color="Red") ):

>	animate(S, [a], a=0..2*Pi, style=surface, lightmodel=light4, orientation=[-120,75,0], labels=[x,y,z], frames=60, axes=framed);

Download rotate-helix-animate.mw

Animation of the ruled surface...

Answered: Rouben Rostamian 8259

March 25 2023

3 4

Here is an animation of the generatrix of your ruled surface.

And here is the worksheet the produced it.

Download ruled-surface.mw

Functions rather than indexed names...

Answered: Rouben Rostamian 8259

March 20 2023

4 0

Finite differences work more naturally in Maple if we use function form rather than indexed name for the variables. Here is an illustration.

Finite differences in Maple

Finite difference calculations work more naturally in Maple

if is expressed as a function .

>

restart;

and

>	dx := (i,j) -> (u(i+1,j) - u(i,j))/h; dy := (i,j) -> (u(i,j+1) - u(i,j))/h;

and

>	(dx(i,j) - dx(i-1,j))/h: simplify(%): d2x := unapply(%, i, j); (dy(i,j) - dy(i,j-1))/h: simplify(%): d2y := unapply(%, i, j);

>	(d2x(i+1,j) - d2x(i,j))/h: simplify(%): d3x := unapply(%, i, j);

The Laplacian

>	d2x(i,j) + d2y(i,j): simplify(%): Lap := unapply(%, i, j);

The x derivative of the Laplacian

>	(Lap(i+1,j) - Lap(i,j))/h: simplify(%);

Download finite-difference-expressions.mw

Making a phase portrait...

Answered: Rouben Rostamian 8259

March 19 2023

2 1

I can't tell whether you have picked the coefficients of the DEs based on some specific need or just randomly. The initial conditions and the plotting range are certainly picked randomly. A phase portrait is interesting near the equilibria. You should indentify the locations of the equilibria and then set up the plotting scene near that location.

Have a look at this modified version of your worksheet.

>

restart;

>	with(plots):

>	with(DEtools):

The vector field:

>	rhs1 := (A,l,S) -> `αAP`A - `μA`A; rhs2 := (A,l,S) -> `αlS`lS - `μl`l - `βlA`lA; rhs3 := (A,l,S) -> `αS`S - `βSl`lS - 0.25SA;

Equilibria:

>	equilibria := solve([rhs1(A,l,S), rhs2(A,l,S), rhs3(A,l,S)], [A,l,S]);

The differential equations:

>	de1 := diff(A(t), t) = rhs1(A(t),l(t),S(t)); de2 := diff(l(t), t) = rhs2(A(t),l(t),S(t)); de3 := diff(S(t), t) = rhs3(A(t),l(t),S(t));

Parameters:

>	`αAP` := 2; `μA` := 1; `αlS` := 2; `μl` := 0.2; `βlA` := 0.5; `αS` := 100; `βSl` := 1; `βSA` = 0.25;

Where are the equilibria?

>	equilibria;

Pick initial conditions near the second equilibrium:

>	ICs := seq([A(0)=0,l(0)=z1,S(0)=0.3], z1=90..110, 10), seq([A(0)=0.1,l(0)=z1,S(0)=0.15], z1=90..110, 10);

How many initial conditions?

>	nops([ICs]);

Specify colors for plotting orbits, one color per initial condition:

>	DEplot3d([de1, de2, de3], [A,l,S], t=0..3, [ICs], linecolor=["Gold", "Cyan", "Magenta", "Red", "Green", "Blue"], thickness=5, stepsize=0.01, arrows=hybrid, dirgrid=[5,5,5], animatecurves=true);

Warning, numpoints adjusted from 301 to 313 for animation

>

Download deplot3d-animate-orbits.mw

This ought to do it...

Answered: Rouben Rostamian 8259

March 19 2023

2 0

Download mw.mw

Solution...

Answered: Rouben Rostamian 8259

March 17 2023

2 2

>

restart;

>	de := (diff(H(K),K))^3+4K^4H(K)^4diff(H(K),K)+8K^4*H(K)^5=0;

The differential equation above is not well-defined. A properly defined first order

differential equation should be of the form dH/dK = something. Therefore we

begin with solving your equation for the derivative:

>	convert(de, D); sol := [ solve(%, D(H)(K)) ];

>	nops(sol);

We see that there are three possible solutions. Two of them are complex-valued.

We discard those and keep what is left:

>	DE := diff(H(K),K) = simplify(remove(has, sol, I))[] assuming K > 0;

diff(H(K), K) = -(1/3)*12^(1/3)*(K^2*H(K)^2*12^(1/3)-K^(2/3)*((4*K^4*H(K)^6+27*H(K)^4)^(1/2)*3^(1/2)-9*H(K)^2)^(2/3))*H(K)*K^(2/3)/((4*K^4*H(K)^6+27*H(K)^4)^(1/2)*3^(1/2)-9*H(K)^2)^(1/3)

This is the differential equation that you should be working with.
Let's solve it:

>	dsol := dsolve({DE, H(0)=1/2}, numeric);

>	plots:-odeplot(dsol, K=0..7);

Download mw.mw

Here it is...

Answered: Rouben Rostamian 8259

March 11 2023

6 5

>

restart;

>	with(plots):

We wish to determine the equation of the surface obtained by rotating
the line : about the line , and plot that surface.

We observe that the rotation axis, may be parameterized as
or equivalently
, where :

>	u := < 1, 0, 1 >; L := t -> <0,0,-1> + u*t;

Vector(3, {(1) = 1, (2) = 0, (3) = 1})

Introduce the mutually orthogonal unit vectors and , both of which
are also orthogonal to :

>	v := <0,1,0>;

Vector(3, {(1) = 0, (2) = 1, (3) = 0})

>	LinearAlgebra:-CrossProduct(u,v): w := %/sqrt(% . %);

Vector(3, {(1) = -(1/2)*sqrt(2), (2) = 0, (3) = (1/2)*sqrt(2)})

Also observe that the line may be parametrized as

>	C := s -> <0,s,s>;

For any point on the line there is a point on the line so that the line segment

is orthogonal to , that is, is the nearest point on the line to the point .

Let , and We determine by asserting that

is perpendicular to :

>	P := <0, s, s>; Q := <t, 0, t-1>; (Q - P)^+ . u = 0; isolate(%, t); Q := eval(Q, %);

Vector(3, {(1) = 0, (2) = s, (3) = s})

Vector(3, {(1) = t, (2) = 0, (3) = t-1})

Vector[column](%id = 36893628651414056100)

As the line rotates about , the point traces a circle about the point .

The circle's radius is the length of :

>	r := sqrt(simplify((Q - P)^+ . (Q-P)));

(1/2)*(6*s^2+4*s+2)^(1/2)

The entire surface of revolution consists of a family of such circles. We conclude that

the parametric equation of the surface of revolution is

>	S := Q + r(vcos(q) + w*sin(q));

$Vector(3, {(1) = 1/2+(1/2)*s-(1/4)*sqrt(6*s^2+4*s+2)*sin(q)*sqrt(2), (2) = (1/2)*sqrt(6*s^2+4*s+2)*cos(q), (3) = -1/2+(1/2)*s+(1/4)*sqrt(6*s^2+4*s+2)*sin(q)*sqrt(2)})$

Plotting

>	display( tubeplot([seq](L(t)), t=-1..2, radius=0.06, color="Red"), # the line L tubeplot([seq](C(s)), s=-2..2, radius=0.06, color="Green"), # the line C plot3d(S, s=-2..2, q=-Pi..Pi, color="Gold"), scaling=constrained, style=surface, lightmodel=light2);

Download surface-of-revolution.mw

Spurious solutions...

Answered: Rouben Rostamian 8259

February 23 2023

0 1

You are right, the equation has no (real) solutions.

If you were solving that equation by hand, you would square both sides and get

x^2 - 10*x + 1 = -8*x^2 + 9*x - 1

and then solve the quadratic and obtain x = 1/9, 2. Then you would substitute these into the original equation to check whether they solve the original equation. If it happens that they don't, they are referred to as spurious solutions and should be thrown out.

In the present case, x = 1/9 and 2 are both spurious. Ideally Maple should check for that but it doesn't. Just be aware of the issue and check it yourself.

Linear or nonlinear?...

Answered: Rouben Rostamian 8259

February 23 2023

2 4

Your equ4 is
Z1T = ZL*(I1T + I2)/I1T

which is equivalent to

Z1T*I1T = ZL*(I1T + I2)

The left-hand side of that equation is the product of two of your unknowns:

vars_unknown := [V1, I1T, V1T, Z1, Z1T]

and therefore your system of equations is nonlinear. It's no wonder that
LinearAlgebra:-GenerateMatrix produces nonsense.

It seems that you believe that the system of equations should be linear. If so, then you must have entered some of the equations incorrectly. Check!

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

8259 Reputation

23 Badges

MaplePrimes Activity

These are answers submitted by Rouben Rostamian

Numerical solution...

Try this...

Correction...

Simplified version...

Multispan beam...

Streamlines...

Rotating helix...

Animation of the ruled surface...

Functions rather than indexed names...

Making a phase portrait...

This ought to do it...

Solution...

Here it is...

Spurious solutions...

Linear or nonlinear?...

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