Rouben Rostamian

Solution...

Answered: Rouben Rostamian 7711

December 14 2018

1 10

The hill is We assume that the conditions are such that

the particle continues pressing against the hill with a positive force

during the motion, that is, it does not jump away or separate

from the surface.

Here is the particle's position vector:

>	r := <x(t),f(x(t))>;

Vector(2, {(1) = x(t), (2) = f(x(t))})

and the velocity vector:

>	v := diff(r,t);

Vector(2, {(1) = diff(x(t), t), (2) = (D(f))(x(t))*(diff(x(t), t))})

The unit tangent vector to the path:

>	v / sqrt(v^+ . v): T := simplify(%) assuming D(x)(t)>0;

$Vector(2, {(1) = 1/sqrt(1+(D(f))(x(t))^2), (2) = (D(f))(x(t))/sqrt(1+(D(f))(x(t))^2)})$

The upward unit normal vector to the path:

>	N := < -T[2], T[1] >;

$Vector(2, {(1) = -(D(f))(x(t))/sqrt(1+(D(f))(x(t))^2), (2) = 1/sqrt(1+(D(f))(x(t))^2)})$

Forces acting on the particle are:
1. The force of gravity (that is, its weight)
2. The reaction force of the surface in the normal direction, of (yet unknown) magnitude
3. The friction force exerted by the surface in the tangential direction, of magnitude

>	F := <0,-mg> + nN - munT;

$Vector(2, {(1) = -n*(D(f))(x(t))/sqrt(1+(D(f))(x(t))^2)-mu*n/sqrt(1+(D(f))(x(t))^2), (2) = -m*g+n/sqrt(1+(D(f))(x(t))^2)-mu*n*(D(f))(x(t))/sqrt(1+(D(f))(x(t))^2)})$

Newton's law of motion:

>	m*diff(v,t) =~ F: DE := convert(%, list);

[m*(diff(diff(x(t), t), t)) = -n*(D(f))(x(t))/(1+(D(f))(x(t))^2)^(1/2)-mu*n/(1+(D(f))(x(t))^2)^(1/2), m*(((D@@2)(f))(x(t))*(diff(x(t), t))^2+(D(f))(x(t))*(diff(diff(x(t), t), t))) = -m*g+n/(1+(D(f))(x(t))^2)^(1/2)-mu*n*(D(f))(x(t))/(1+(D(f))(x(t))^2)^(1/2)]

Eliminate the unknown between the equations:

>	solve(DE, {diff(x(t),t,t), n}): select(has, %, diff(x(t),t,t)): de := %[];

diff(diff(x(t), t), t) = -((D(f))(x(t))+mu)*(((D@@2)(f))(x(t))*(diff(x(t), t))^2+g)/(1+(D(f))(x(t))^2)

Here I am taking the hill to be a parabola for the sake of illustration. Change as needed.

>	myde := eval(de, {f = (x -> x^2), m=1, g=1, mu=1/10});

diff(diff(x(t), t), t) = -(2*x(t)+1/10)*(2*(diff(x(t), t))^2+1)/(1+4*x(t)^2)

>	ic := x(0)=0, D(x)(0)=1;

>	tmax := 2;

>	dsol := dsolve({myde, ic}, numeric, output=listprocedure, events=[[diff(x(t),t), halt]], range=0..tmax);

Warning, cannot evaluate the solution further right of 1.3669169, event #1 triggered a halt

The graph of versus

>	plots:-odeplot(dsol, [t,x(t)]);

Animation:

>	plots:-odeplot(dsol, [x(t), x(t)^2], scaling=constrained, frames=100, labels=[x,y]);

Download sliding.mw

Here's how...

Answered: Rouben Rostamian 7711

December 12 2018

2 5

You didn't specify the PDE, so I will explain this in the context of the heat equation ∂u/∂t = ∂²u/∂x². You plug u(x,t)=X(x) T(t) into the PDE, and after some work arrive at an ordinary differential equation for T(t), let's say something like T_n'(t) + a_n T_n(t) = 0, where T_n' is the derivative of T_n. You solve that ODE, and then express the solution of the PDE as u(x,t)=sum(X_n(x) T_n(t)). It follows that the derivative of u(x,t) with respect to t is sum(X_n(x) T_n'(t)). But you know T_n'(t) from the ODE, to wit, T_n'(t) = − a_n T_n(t). We conclude that ∂/∂t u = −sum(a_n X_n(t) T_n(t)).

Another answer...

Answered: Rouben Rostamian 7711

November 01 2018

1 1

As an alternative to Carl Love's implicitdiff approach, here is one with plain diff. I haven't compared the results but they ought to be equivalent.

>

restart;

>	w := exp(exp(xb)(r-1)/(1+varphiexp(xb))); # inserted "*" before (r-1)

exp(exp(x*b)*(r-1)/(1+varphi*exp(x*b)))

>	constraint := ln(r) = varphi(exp(xb))(r-1)/(1+varphiexp(x*b))-1;

ln(r) = varphi*exp(x*b)*(r-1)/(1+varphi*exp(x*b))-1

Rearrange the constraint equation:

>	my_phi := isolate(constraint, varphi);

Eliminte varphi in w. Call the result W (uppercase w):

>	W := eval(w, my_phi);

exp(exp(x*b)*(r-1)/(1-(ln(r)+1)/(-r+ln(r)+2)))

Now we may compute the derivative of W with respect to b:

>	simplify(diff(W, b));

To compute the derivative of W with respect to varphi, we note
that dW/`dφ` = dW/dr*(dr/`dφ`)

>	simplify(diff(W,r)/diff(rhs(my_phi),r));

-exp((r-2)*exp(x*b)+2*x*b)*(-r+ln(r)+2)^2*(r-1)*r^(-exp(x*b))/(ln(r)*r+1)

Animation...

Answered: Rouben Rostamian 7711

October 25 2018

5 7

Ronan gave a link to my post from almost two years ago where I illustrated the use of quaternions in tracking rotations in space. That motivated me to clean up the code a bit and then apply it to the phenomenon of the flipping T-handle which is the subject of the current thread. The code, which combines Euler's equations of motion and quaternions, is sufficiently commented to let you adapt it to other situations.

Here is the worksheet T-handle-flip.mw and here is the resulting animation:

Fourier series produces the weak solutio...

Answered: Rouben Rostamian 7711

October 23 2018

5 5

@ecterrab You are correct in stating that this problem has "no solution". To be precise, this problem has no classical solution.

A great deal of effort was devoted throughout the 20th century to make sense of exactly such boundary value problems which arise in just about every application, especially in physics. The concept of generalized derivatives and weak solutions were introduced precisely for that purpose. If we dismiss weak solutions, we will be throwing away most of the useful applications of the PDEs.

Consider, for instance, a slender (one-dimensional) bar currently at a uniform temperature of 100 degrees. We bring one end of it in contact with an ice cube. We want to determine how the bar's temperature evolves in subsequent times.

The temperature, u(x,t), satisfies the PDE u_t = u_xx, where subscripts denote differentiation. The initial condition is u(x,0)=100. The boundary condition is u(0,t)=0. Certainly these two conditions are incompatible but we don't want throw up our hands and say that this problem is not solvable. True, no classical solution exists, but a weak solution does, and that's what's relevant to this question.

Mathematica's solution that was posted earlier (and is now deleted) expresses the solution of the PDE in the Fourier series. The Fourier series solution picks up the weak solution, therefore that solution was correct in the sense that it produced the expected weak solution to the problem. In the worksheet below I show how to reproduce Mathematica's solution with Maple's pdsolve(). The process hits a few snags but eventually we arrive at the desired solution.

P.S.: To be clear, I am not here to berate the terrific job that you have done (and continue doing) with extending Maple's capabilities in many directions. All I want to point out is that the Fourier series picks up the weak solution, and that it should not be dismissed.

[Worksheet produced in Maple 2018.1]

>

restart;

>	pde := diff(u(x,y),x,x) + diff(u(x,y),y,y) = 0;

>	bc := u(0,y)=0, u(Pi,y)=sinh(Pi)*cos(y), u(x,0)=sin(x), u(x,Pi)=-sinh(x);

>	result := pdsolve({pde,bc});

Snag #1: pdsolve() returns empty. I follow Robert Lopez's lead and split up

the problem into four problems, then use superposition:

>	bc1 := u(0,y)=0, u(Pi,y)=0, u(x,0)=0, u(x,Pi)=0;

>	bc2 := u(0,y)=0, u(Pi,y)=sinh(Pi)*cos(y), u(x,0)=0, u(x,Pi)=0;

>	bc3 := u(0,y)=0, u(Pi,y)=0, u(x,0)=sin(x), u(x,Pi)=0;

>	bc4 := u(0,y)=0, u(Pi,y)=0, u(x,0)=0, u(x,Pi)=-sinh(x);

Attempt to solve the problem with bc1:

>	sol1 := pdsolve({pde,bc1});

Snag #2: The exact solution of {pde,bc1} is u(x,y)=0. Maple fails to detect that.
I don't know why. The boundary conditions are certainly compatible.

Here I intervene and set the solution myself:

>	sol1 := u(x,y) = 0;

Atempt to solve with bc2:

>	sol2_tmp := pdsolve({pde,bc2});

u(x, y) = Sum(2*n*sinh(Pi)*((-1)^n+1)*sin(y*n)*(exp(2*x*n)-1)*exp(n*(Pi-x))/(Pi*(exp(2*Pi*n)-1)*(n^2-1)), n = 1 .. infinity)

Snag #3: The first term (n=1) of the summation above evaluates to 0/0.

In fact, all terms with odd n should evaluate to zero. Here we intervene and

replace n by 2*n:

>	applyop(eval, [2,1], sol2_tmp, n=2*n): sol2 := simplify(%) assuming n::posint;

u(x, y) = Sum(8*n*sinh(Pi)*sin(2*y*n)*(exp(4*x*n)-1)*exp(2*n*(Pi-x))/(Pi*(exp(4*Pi*n)-1)*(4*n^2-1)), n = 1 .. infinity)

The boundary conditions bc3 and bc4 give no problem:

>	sol3 := pdsolve({pde,bc3});

u(x, y) = sin(x)*(exp(2*Pi-y)-exp(y))/(exp(2*Pi)-1)

>	sol4 := pdsolve({pde,bc4});

u(x, y) = Sum((-1)^n*n*(exp(2*Pi)-1)*sin(x*n)*(exp(2*y*n)-1)*exp((-y+Pi)*n-Pi)/(Pi*(n^2+1)*(exp(2*Pi*n)-1)), n = 1 .. infinity)

The original problem's weak solution is obtained as a superposition

of the previously solved four sub-problems. Ideally, Maple's pdsolve()

should obtain this solution in one shot and I am confident that at some future time it will:

>	sol := rhs(sol1) + rhs(sol2) + rhs(sol3) + rhs(sol4);

Sum(8*n*sinh(Pi)*sin(2*y*n)*(exp(4*x*n)-1)*exp(2*n*(Pi-x))/(Pi*(exp(4*Pi*n)-1)*(4*n^2-1)), n = 1 .. infinity)+sin(x)*(exp(2*Pi-y)-exp(y))/(exp(2*Pi)-1)+Sum((-1)^n*n*(exp(2*Pi)-1)*sin(x*n)*(exp(2*y*n)-1)*exp((-y+Pi)*n-Pi)/(Pi*(n^2+1)*(exp(2*Pi*n)-1)), n = 1 .. infinity)

We truncate the series:

>	my_sol := eval(sol, infinity=50);

Sum(8*n*sinh(Pi)*sin(2*y*n)*(exp(4*x*n)-1)*exp(2*n*(Pi-x))/(Pi*(exp(4*Pi*n)-1)*(4*n^2-1)), n = 1 .. 50)+sin(x)*(exp(2*Pi-y)-exp(y))/(exp(2*Pi)-1)+Sum((-1)^n*n*(exp(2*Pi)-1)*sin(x*n)*(exp(2*y*n)-1)*exp((-y+Pi)*n-Pi)/(Pi*(n^2+1)*(exp(2*Pi*n)-1)), n = 1 .. 50)

and plot the result:

>	plot3d(my_sol, x=0..Pi, y=0..Pi, style=patchcontour, contours=20);

Download pde-bc.mw

Maybe this?...

Answered: Rouben Rostamian 7711

October 17 2018

0 0

I don't quite understand what you mean by "link a worksheet written in Maple 16 into WordPress" but let me try.

WordPress can display HTML code as well as images. In Malpe's menubar select File->Print and save your worksheet as a graphics image. On my Maple 2018 on Linux the image is saved in the PostScript format. In your Maple 16 it may be something else, but it does not matter much. You may convert one image into any other image format such as PNG, which then you pass on to WordPress.

Here it is...

Answered: Rouben Rostamian 7711

October 13 2018

0 1

plot([0, sin(x)], x=0..2*Pi, axes=boxed, color=[black,blue]);

Use the Layout palette...

Answered: Rouben Rostamian 7711

October 10 2018

1 0

In the Layout palette pick the entry that consists of a green "A" on top of a red "b". Then edit to replace the "A" with "max" and the "B" with the desired subscript. Here is what we get:

This works in Maple 2018. I don't have Maple 2015 to test.

Download worksheet: mw.mw

It depends.......

Answered: Rouben Rostamian 7711

October 08 2018

0 0

The answer depends on what you assume about x.

Is x itself bounded? If yes, then cos(i*x) is bounded.

Is x imaginary and unbounded? If yes, then i*x is real, and therefore cos(i*x) is bounded because sin(i*x)^2 + cos(i*x)^2 = 1.

Is x real and unbounded? If so, then:

convert(cos(I*x), exp) assuming x::real;

shows that cos(i*x) is 1/2 exp(x) + 1/2 exp(-x), and therefore cos(i*x) is unbounded.

And here is an exercise for you: Suppose x is complex. What can you tell about the boundedness of cos((i*x)?

There is no nontrivial solution...

Answered: Rouben Rostamian 7711

September 28 2018

3 7

Need differential equations...

Answered: Rouben Rostamian 7711

September 09 2018

0 0

The purpose of phaseportrait is to plot solutions of a system of differential equations. There is no differential equation in what you have shown.

Suggestion: Maple's help page on phaseportrait is quite limited and not very instructive. To its credit, it refer the reader to DEplot's help page which does what phaseportrait is supposed to do, and is much better documented. Look up the examples in DEplot's help page and come back here if you still have questions.

Solution...

Answered: Rouben Rostamian 7711

September 02 2018

3 1

The abs(t) seems to confuse Maple. The following equivalent alternative works:

Download mw.mw

Linearization...

Answered: Rouben Rostamian 7711

August 29 2018

3 0

>

restart:

>	eq[1]:=diff(x[1](t), t)-2x[1](t)(N1N2N3-x[1](t)x[3](t))/(N1N2*N6):

>

eq[2]:=diff(x[2](t), t)+((1-N1*N2*N6/x[1](t))*(N1*N2*N3-x[2](t)*x[3](t))/((N1*N2+x[2](t)*(1-N1*N2*N6/x[1](t)))*N1*N2)-2*x[2](t)*(1-N1*N2*N6/x[1](t))*(N1*N2*N3-x[1](t)*x[3](t))/(x[1](t)*N1*N2*(N1*N2+x[2](t)*(1-N1*N2*N6/x[1](t)))))/((1-N1*N2*N6/x[1](t))/(x[2](t)*(N1*N2+x[2](t)*(1-N1*N2*N6/x[1](t))))-ln(1+(1-N1*N2*N6/x[1](t))*x[2](t)/(N1*N2))/x[2](t)^2):

>

eq[3]:=diff(x[3](t), t)-(-(-2*N1*N2*N3*ln(1+(1-N1*N2*N6/x[1](t))*x[2](t)/(N1*N2))*(N1*N2*N3-x[1](t)*x[3](t))/x[1](t)-N3*(x[3](t)*N2*N1*ln(1+(1-N1*N2*N6/x[1](t))*x[2](t)/(N1*N2))/x[2](t)+N3*N2*N1*(1-N1*N2*N6/x[1](t)-N2*N1*ln(1+(1-N1*N2*N6/x[1](t))*x[2](t)/(N1*N2))/x[2](t))/x[2](t))-N5*(x[3](t)^2*N2*N1*ln(1+(1-N1*N2*N6/x[1](t))*x[2](t)/(N1*N2))/x[2](t)+2*x[3](t)*N3*N2*N1*(1-N1*N2*N6/x[1](t)-N2*N1*ln(1+(1-N1*N2*N6/x[1](t))*x[2](t)/(N1*N2))/x[2](t))/x[2](t)+N3^2*N2*N1*((1/2)*(1-N1*N2*N6/x[1](t))^2-N2*N1*(1-N1*N2*N6/x[1](t)-N2*N1*ln(1+(1-N1*N2*N6/x[1](t))*x[2](t)/(N1*N2))/x[2](t))/x[2](t))/x[2](t))-N2*N1*ln(1+(1-N1*N2*N6/x[1](t))*x[2](t)/(N1*N2))/(x[2](t)*F1)-Kin*x[3](t)^2-Kex*N1*N2*(x[3](t)+N3*(1-N1*N2*N6/x[1](t)))^2/(N1*N2+x[2](t)*(1-N1*N2*N6/x[1](t))))*x[1](t)*(N1*N2*N6*x[2](t)-N1*N2*x[1](t)-x[2](t)*x[1](t))*F1/(N1*N2)-Aol^2*Dol*F1*LR*N5*N3^2*N1^2*N2^2*N6^2+2*Aol^2*Dol*F1*LR*N5*N3^2*N1*N2*N6*x[1](t)+2*Aol^2*Dol*F1*LR*N5*N3*N1*N2*N6*x[1](t)*x[3](t)-Aol^2*Dol*F1*LR*N5*N3^2*x[1](t)^2-2*Aol^2*Dol*F1*LR*N5*N3*x[1](t)^2*x[3](t)-Aol^2*Dol*F1*LR*N5*x[1](t)^2*x[3](t)^2-2*Aol*F1*LR*N3^2*N1*N2*x[1](t)+2*Aol*F1*LR*N3*x[1](t)^2*x[3](t)-LR*x[1](t)^2+F1*N4*N1*N2*N6^2*x[3](t)^2*x[2](t)-F1*N4*N1*N2*N6*x[1](t)*x[3](t)^2-F1*N4*N6*x[1](t)*x[3](t)^2*x[2](t)+N1*N2*N6^2*x[2](t)-N1*N2*N6*x[1](t)-N6*x[1](t)*x[2](t))/(Aol*F1*LR*x[1](t)^2+F1*N6*x[1](t)*x[2](t)+F1*ln(-(N1*N2*N6*x[2](t)-N1*N2*x[1](t)-x[2](t)*x[1](t))/(N1*N2*x[1](t)))*x[1](t)^2*x[2](t)-F1*N1*N2*N6^2*x[2](t)+F1*N1*N2*N6*x[1](t)+F1*N1*N2*ln(-(N1*N2*N6*x[2](t)-N1*N2*x[1](t)-x[2](t)*x[1](t))/(N1*N2*x[1](t)))*x[1](t)^2-F1*N1*N2*N6*ln(-(N1*N2*N6*x[2](t)-N1*N2*x[1](t)-x[2](t)*x[1](t))/(N1*N2*x[1](t)))*x[1](t)*x[2](t)):

A differential equation is expected to be something of the form dx/dt = F(x(t)). That is, one defines the function F first, and then defines the differential equation afterwards in terms of F. In your equations the function F is hidden, so we first untangle your notation to set things straight. This should not have been necessary if you had done things correctly in the beginning by defining the F functions first.

>	isolate(eq[1], diff); eval(rhs(%), [x[1](t)=x[1], x[2](t)=x[2], x[3](t)=x[3]]); F[1] := unapply(%, [x[1], x[2], x[3]]); # here is F Eq[1] := diff(x[1](t),t) = F[1](x[1](t), x[2](t), x[3](t)); # here is the DE

(1)

Similarly, we untangle the notation for eq[2] and eq[3]:

>	isolate(eq[2], diff): eval(rhs(%), [x[1](t)=x[1], x[2](t)=x[2], x[3](t)=x[3]]): F[2] := unapply(%, [x[1], x[2], x[3]]): # here is F Eq[2] := diff(x[1](t),t) = F[2](x[1](t), x[2](t), x[3](t)): # here is the DE

>

>	isolate(eq[3], diff): eval(rhs(%), [x[1](t)=x[1], x[2](t)=x[2], x[3](t)=x[3]]): F[3] := unapply(%, [x[1], x[2], x[3]]): # here is F Eq[3] := diff(x[1](t),t) = F[3](x[1](t), x[2](t), x[3](t)): # here is the DE

>

Find the equilibria:

>	solve({F[1](x[1],x[2],x[3])=0, F[2](x[1],x[2],x[3])=0, F[3](x[1],x[2],x[3])=0}, {x[1],x[2],x[3]});

${x[1] = N1*N2*N6+N1*N2*exp(RootOf(6*LR*N6+2*F1*Kex*N3^2*(exp(_Z))^3+F1*N3^2*N5*(exp(_Z))^3-2*F1*Kex*N3^2*(exp(_Z))^2+2*F1*Kin*N3^2*(exp(_Z))^2+2*F1*N3^2*N6*(exp(_Z))^2-2*F1*Kin*N3^2*exp(_Z)-F1*N3^2*N5*exp(_Z)-2*F1*N3^2*N6*exp(_Z)+2*F1*Kin*N3^2*N6*exp(_Z)+2*F1*N3^2*N4*N6*exp(_Z)+2*F1*Kex*N3^2*N6*(exp(_Z))^2-2*LR+2*Aol^2*Dol*F1*LR*N3^2*N5*(exp(_Z))^3-2*Aol^2*Dol*F1*LR*N3^2*N5*(exp(_Z))^2+2*(exp(_Z))^3*_Z-4*(exp(_Z))^2*_Z+2*_Z*exp(_Z)+2*N6*(exp(_Z))^3-6*LR*(exp(_Z))^2-4*N6^2*exp(_Z)-4*N6*(exp(_Z))^2+6*LR*exp(_Z)+2*N6*exp(_Z)+2*LR*(exp(_Z))^3+2*N6^3*exp(_Z)+4*N6^2*(exp(_Z))^2+2*Aol^2*Dol*F1*LR*N3^2*N5*N6*(exp(_Z))^2+2*LR*N6^3-6*LR*N6^2+2*F1*N3^2*(exp(_Z))^3-4*F1*N3^2*(exp(_Z))^2+2*F1*N3^2*exp(_Z)+6*LR*N6^2*exp(_Z)+6*LR*N6*(exp(_Z))^2-12*LR*N6*exp(_Z)+2*N6^2*exp(_Z)*_Z+4*N6*(exp(_Z))^2*_Z-4*N6*exp(_Z)*_Z))-N1*N2, x[2] = N1*N2*N6+N1*N2*exp(RootOf(6*LR*N6+2*F1*Kex*N3^2*(exp(_Z))^3+F1*N3^2*N5*(exp(_Z))^3-2*F1*Kex*N3^2*(exp(_Z))^2+2*F1*Kin*N3^2*(exp(_Z))^2+2*F1*N3^2*N6*(exp(_Z))^2-2*F1*Kin*N3^2*exp(_Z)-F1*N3^2*N5*exp(_Z)-2*F1*N3^2*N6*exp(_Z)+2*F1*Kin*N3^2*N6*exp(_Z)+2*F1*N3^2*N4*N6*exp(_Z)+2*F1*Kex*N3^2*N6*(exp(_Z))^2-2*LR+2*Aol^2*Dol*F1*LR*N3^2*N5*(exp(_Z))^3-2*Aol^2*Dol*F1*LR*N3^2*N5*(exp(_Z))^2+2*(exp(_Z))^3*_Z-4*(exp(_Z))^2*_Z+2*_Z*exp(_Z)+2*N6*(exp(_Z))^3-6*LR*(exp(_Z))^2-4*N6^2*exp(_Z)-4*N6*(exp(_Z))^2+6*LR*exp(_Z)+2*N6*exp(_Z)+2*LR*(exp(_Z))^3+2*N6^3*exp(_Z)+4*N6^2*(exp(_Z))^2+2*Aol^2*Dol*F1*LR*N3^2*N5*N6*(exp(_Z))^2+2*LR*N6^3-6*LR*N6^2+2*F1*N3^2*(exp(_Z))^3-4*F1*N3^2*(exp(_Z))^2+2*F1*N3^2*exp(_Z)+6*LR*N6^2*exp(_Z)+6*LR*N6*(exp(_Z))^2-12*LR*N6*exp(_Z)+2*N6^2*exp(_Z)*_Z+4*N6*(exp(_Z))^2*_Z-4*N6*exp(_Z)*_Z))-N1*N2, x[3] = N3/(N6+exp(RootOf(6*LR*N6+2*F1*Kex*N3^2*(exp(_Z))^3+F1*N3^2*N5*(exp(_Z))^3-2*F1*Kex*N3^2*(exp(_Z))^2+2*F1*Kin*N3^2*(exp(_Z))^2+2*F1*N3^2*N6*(exp(_Z))^2-2*F1*Kin*N3^2*exp(_Z)-F1*N3^2*N5*exp(_Z)-2*F1*N3^2*N6*exp(_Z)+2*F1*Kin*N3^2*N6*exp(_Z)+2*F1*N3^2*N4*N6*exp(_Z)+2*F1*Kex*N3^2*N6*(exp(_Z))^2-2*LR+2*Aol^2*Dol*F1*LR*N3^2*N5*(exp(_Z))^3-2*Aol^2*Dol*F1*LR*N3^2*N5*(exp(_Z))^2+2*(exp(_Z))^3*_Z-4*(exp(_Z))^2*_Z+2*_Z*exp(_Z)+2*N6*(exp(_Z))^3-6*LR*(exp(_Z))^2-4*N6^2*exp(_Z)-4*N6*(exp(_Z))^2+6*LR*exp(_Z)+2*N6*exp(_Z)+2*LR*(exp(_Z))^3+2*N6^3*exp(_Z)+4*N6^2*(exp(_Z))^2+2*Aol^2*Dol*F1*LR*N3^2*N5*N6*(exp(_Z))^2+2*LR*N6^3-6*LR*N6^2+2*F1*N3^2*(exp(_Z))^3-4*F1*N3^2*(exp(_Z))^2+2*F1*N3^2*exp(_Z)+6*LR*N6^2*exp(_Z)+6*LR*N6*(exp(_Z))^2-12*LR*N6*exp(_Z)+2*N6^2*exp(_Z)*_Z+4*N6*(exp(_Z))^2*_Z-4*N6*exp(_Z)*_Z))-1)}$

(2)

That's a mess. We back up and do things manually as much as possible. We have:

>	F[1](x[1],x[2],x[3]) = 0; solve(%, x[1]); tmp1 := %[2]; # ignore the solution x[1]=0; pick the other one and call it tmp1

(3)

Similarly, solve F[2] =0 for x[2]:

>	simplify(F[2](tmp1,x[2],x[3])); tmp2 := solve(%, x[2]); # solve for x[2] and call the solution tmp2

(-N6*x[3]+N3)*(N1*N2*N3-x[2]*x[3])*x[2]^2/(N2*N1*((N1*N2*N3+x[2]*(-N6*x[3]+N3))*ln(((N1*N2+x[2])*N3-N6*x[2]*x[3])/(N1*N2*N3))-x[2]*(-N6*x[3]+N3)))

(4)

It remians to solve the equation F[3]=0. But it does not have a solution in terms of elementary function:

>	simplify(F[3](tmp1, tmp2, x[3]));

(5)

So we seek a numerical solution. Toward that end, we make a list of the problem's parameters:

>	indets({F[1](x[1],x[2],x[3]), F[2](x[1],x[2],x[3]), F[3](x[1],x[2],x[3])}, name);

(6)

Pick numerical values for the parameters:

>	params:= F1=1.493520000, LR=0.3178477690, Aol=0.7843287882, Dol= 0.5128205128e-1, N4= 3.670047218, N5= 7.340094435, N1= 1.054254618, N2= 0.5146253739e-1, N6=0.650636, N3=2.64407, Kin=42.48, Kex=0.68

(7)

Let's see what F[3] looks like as a function of x[3]:

>	plot(eval(F[3](tmp1, tmp2, x[3]), {params}), x[3]=0..200);

We see that there is a root near x[3]=150. We apply fsolve() to obtain the precise value. We call it X[3]:

>	eval(F[3](tmp1, tmp2, x[3]), {params}): X[3] := fsolve(%, x[3]);

(8)

With that in hand, we may evaluate the other two coordinates of the equilibrium:

>	X[1] := eval(tmp1, {params, x[3]=X[3]});

(9)

>	X[2] := eval(tmp2, {params, x[3]=X[3]});

(10)

Sanity check: Let's verify that we have made no mistakes:

>	eval([F[1](X[1],X[2],X[3]), F[2](X[1],X[2],X[3]), F[3](X[1],X[2],X[3])], {params});

(11)

That's not bad. The third value is not as small as it can be, but if necessary, we may increase the Digits value in order to obtain greater accuracy.

Aside: With the hindsight obtained from our manual calculation above, we may ask Maple to find the equilibrium at once through a single application of the fsolve() command:

>	eval({F[1](x[1],x[2],x[3])=0, F[2](x[1],x[2],x[3])=0, F[3](x[1],x[2],x[3])=0}, {params}): fsolve(%, {x[1]=0..1e-2,x[2]=0..1e-2,x[3]=140..160});

(12)

which is identical to what we obtained earlier.

>

Now we linearize the system about the equilibrium. We begin by evaluating F[1], F[2] and F[3] at the parameter values:

>	G := eval([F[1](x[1],x[2],x[3]), F[2](x[1],x[2],x[3]), F[3](x[1],x[2],x[3])], {params}):

Here is the Jacobian matrix at any x[1], x[2], x[3]:

>	J := < diff(G[1],x[1]), diff(G[1],x[2]), diff(G[1],x[3]); diff(G[2],x[1]), diff(G[2],x[2]), diff(G[2],x[3]); diff(G[3],x[1]), diff(G[3],x[2]), diff(G[3],x[3]) >:

Evaluate the Jacobian at the equilibrium to obtain the linearized matrix:

>	A := eval(J, {x[1]=X[1], x[2]=X[2], x[3]=X[3]});

_rtable[18446883888842582486]

(13)

Find the eigenvalues of the linearized matrix

>	LinearAlgebra:-Eigenvalues(A);

_rtable[18446883888835980694]

(14)

We have two positive eigenvalues, indicating that the equilibrium is unstable.

>

Download mw2.mw

Solve a differential equation...

Answered: Rouben Rostamian 7711

August 20 2018

0 5

Let me write L instead of LL for greater legibility. By differentiating your expression for L we get:

Solve this as a differential equation:
dsol := dsolve({de, ic}, numeric);

Then plot it:

plots:-odeplot(dsol), phi=0..2*Pi);

To read of the value of the solution at phi=4, do:

dsol(4);

Here is one way...

Answered: Rouben Rostamian 7711

August 20 2018

1 0

iR := rand(1..nR):
iC := rand(1..nC):
C(iR(),iC()); # pick a random element of C
C(iR(),iC()); # pick another random element of C

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

Solution...

Here's how...

Another answer...

Animation...

Fourier series produces the weak solutio...

Maybe this?...

Here it is...

Use the Layout palette...

It depends.......

There is no nontrivial solution...

Need differential equations...

Solution...

Linearization...

Solve a differential equation...

Here is one way...

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