Advanced Search

Rouben Rostamian

7711 Reputation

22 Badges

19 years, 3 days
Contact Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

Here is a phase portrait...

Rouben Rostamian 7711
June 08 2020
3 8
> 

restart;

Don't begin with the differential equations.  Define their right-hand sides first:

> 

f[0] := (x,y1,y2) -> alpha*x*(1 - x/N) - beta[1]*sqrt(x)*y1/(1 + h[1]*beta[1]*sqrt(x)) - beta[2]*sqrt(x)*y2/(1 + h[2]*beta[2]*sqrt(x)) - d*E*x;

proc (x, y1, y2) options operator, arrow; alpha*x*(1-x/N)-beta[1]*sqrt(x)*y1/(1+h[1]*beta[1]*sqrt(x))-beta[2]*sqrt(x)*y2/(1+h[2]*beta[2]*sqrt(x))-d*E*x end proc

> 

f[1] := (x,y1,y2) -> -omega[1]*y1 + Mu[1]*beta[1]*sqrt(x)*y1*(1 - y1/(beta[1]*sqrt(x)))/(1 + h[1]*beta[1]*sqrt(x));

proc (x, y1, y2) options operator, arrow; -omega[1]*y1+Mu[1]*beta[1]*sqrt(x)*y1*(1-y1/(beta[1]*sqrt(x)))/(1+h[1]*beta[1]*sqrt(x)) end proc

> 

f[2] := (x,y1,y2) -> -omega[2]*y2 + Mu[2]*beta[2]*sqrt(x)*y2*(1 - y2/(beta[2]*sqrt(x)))/(1 + h[2]*beta[2]*sqrt(x));

proc (x, y1, y2) options operator, arrow; -omega[2]*y2+Mu[2]*beta[2]*sqrt(x)*y2*(1-y2/(beta[2]*sqrt(x)))/(1+h[2]*beta[2]*sqrt(x)) end proc

Now define the differential equations:

> 

de[0] := diff(x(t),t) = f[0](x(t),y1(t),y2(t));

diff(x(t), t) = alpha*x(t)*(1-x(t)/N)-beta[1]*x(t)^(1/2)*y1(t)/(1+h[1]*beta[1]*x(t)^(1/2))-beta[2]*x(t)^(1/2)*y2(t)/(1+h[2]*beta[2]*x(t)^(1/2))-d*E*x(t)

> 

de[1] := diff(y1(t),t) = f[1](x(t),y1(t),y2(t));

diff(y1(t), t) = -omega[1]*y1(t)+Mu[1]*beta[1]*x(t)^(1/2)*y1(t)*(1-y1(t)/(beta[1]*x(t)^(1/2)))/(1+h[1]*beta[1]*x(t)^(1/2))

> 

de[2] := diff(y2(t),t) = f[2](x(t),y1(t),y2(t));

diff(y2(t), t) = -omega[2]*y2(t)+Mu[2]*beta[2]*x(t)^(1/2)*y2(t)*(1-y2(t)/(beta[2]*x(t)^(1/2)))/(1+h[2]*beta[2]*x(t)^(1/2))

Numbers come last:

> 

params := alpha = 0.75, omega[1] = 0.15, omega[2] = 0.20, N = 0.8, beta[1] = 0.65, beta[2] = 0.70, Mu[1] = 0.75, Mu[2] = 0.74, h[1] = 0.005, h[2] = 0.004, d = 0.05, E = 0.35;

alpha = .75, omega[1] = .15, omega[2] = .20, N = .8, beta[1] = .65, beta[2] = .70, Mu[1] = .75, Mu[2] = .74, h[1] = 0.5e-2, h[2] = 0.4e-2, d = 0.5e-1, E = .35

Here is the system:

> 

sys := eval({de[0], de[1],de[2]}, {params});

{diff(x(t), t) = .75*x(t)*(1-1.250000000*x(t))-.65*x(t)^(1/2)*y1(t)/(1+0.325e-2*x(t)^(1/2))-.70*x(t)^(1/2)*y2(t)/(1+0.280e-2*x(t)^(1/2))-0.175e-1*x(t), diff(y1(t), t) = -.15*y1(t)+.4875*x(t)^(1/2)*y1(t)*(1-1.538461538*y1(t)/x(t)^(1/2))/(1+0.325e-2*x(t)^(1/2)), diff(y2(t), t) = -.20*y2(t)+.5180*x(t)^(1/2)*y2(t)*(1-1.428571429*y2(t)/x(t)^(1/2))/(1+0.280e-2*x(t)^(1/2))}

Need some trial-and-error, and some insight, to find a good range for plotting.

> 

xmin, xmax := 0.1, 0.6;
y1min, y1max := 0.0, 0.6;
y2min, y2max := 0.1, 0.2;

.1, .6

0., .6

.1, .2

Select a good set of initial conditions:

> 

ic := seq([x(0)=xmax, y1(0)=y1max, y2(0)=(1-i/5)*y2min+i/5*y2max], i=0..5),
      seq([x(0)=xmax, y1(0)=0.5*y1max, y2(0)=(1-i/5)*y2min+i/5*y2max], i=0..5),
      seq([x(0)=0.5*xmax, y1(0)=y1max, y2(0)=(1-i/5)*y2min+i/5*y2max], i=0..5),
      seq([x(0)=xmax, y1(0)=0.1, y2(0)=(1-i/5)*y2min+i/5*y2max], i=0..5):

Now plot:

> 

DEtools:-DEplot3d(sys, [x(t),y1(t),y2(t)], t=0..40, stepsize=0.1,
    [ic], x=xmin..xmax, y1=y1min..y1max, y2=y2min..y2max,
    arrows=cheap, linecolor=black, thickness=1, orientation=[30,72,0]);

We see that there is a stable equilibrium near x=0.3, y1=0.2, y2=0.15.  Let's find it:

> 

eval({f[0](x,y1,y2),f[1](x,y1,y2),f[2](x,y1,y2)}, {params}):
fsolve(%, {x=0.3, y1=0.2, y2=0.15});

{x = .3697176032, y1 = .1948334234, y2 = .1549004442}

 

 

 

Download mw.mw

 

The Method of Lines...

Rouben Rostamian 7711
May 27 2020
4 2

I will show how to solve the PDE with the Method of Lines.  For no particular reason I have changed all your coefficients to "1"s, and have added boundary conditions which are missing in your worksheet.

Additionally, I should point out that you may solve the PDE directly with Maple's pdsolve. In this worksheet I solve and plot the solutions both with pdsolve, and with Method of Lines, and of course the two solutions agree.

> 

restart;

 1.  Solving the PDE with Maple's pdsolve

 

Here is how to solve the PDE with Maple's pdsolve.

> 

pde := diff(u(x,t),t) = - diff(u(x,t),x$2) - diff(u(x,t),x$4) - u(x,t)*diff(u(x,t),x);

diff(u(x, t), t) = -(diff(diff(u(x, t), x), x))-(diff(diff(diff(diff(u(x, t), x), x), x), x))-u(x, t)*(diff(u(x, t), x))

> 

ic := u(x,0) = sin(Pi*x/L);

u(x, 0) = sin(Pi*x/L)

> 

bc := u(0,t)=0, D[1,1](u)(0,t)=0, u(L,t)=0, D[1,1](u)(L,t)=0;

u(0, t) = 0, (D[1, 1](u))(0, t) = 0, u(L, t) = 0, (D[1, 1](u))(L, t) = 0

> 

pdsol := pdsolve(pde, eval({ic,bc}, L=1), numeric, spacestep=0.01);

_m140335696338400

> 

pdsol:-plot3d(t=0..0.020);

 

 

 2.  Solving the PDE with the Method of Lines

 

We divide the interval [0, L] into n >= 4equal subintervals and write h = L/n

for the length of each subinterval.  Let x__i = ih, i = 0 .. n be the coordinates

of the endpoints (the nodes) of the intervals .  Note that there are
n+1 nodes altogether, with nodes 0 and n on the boundary, and nodes
1 through n-1 in the interior.  In my comments below, I write u_x
for the derivative of u(x, t) with respect to "x."

 

The finite difference discretizations seen below are pretty standard

in numerical analysis.

> 

# discretize u_x, i=1,..,n-1
d1 := add([-1,0,1] *~ [u[i+k](t) $k=-1..1])/(2*h);

(1/2)*(-u[i-1](t)+u[i+1](t))/h

> 

# discretize u_xx, i=1,..,n-1
d2 := add([1,-2,1] *~ [u[i+k](t) $k=-1..1])/h^2;

(u[i-1](t)-2*u[i](t)+u[i+1](t))/h^2

> 

# discretize u_xx at i=0
d2_left := add([2,-5,4,-1] *~ [u[k](t) $k=0..3])/h^2;

(2*u[0](t)-5*u[1](t)+4*u[2](t)-u[3](t))/h^2

> 

# discretize u_xx at i=n
d2_right := add([-1,4,-5,2] *~ [u[n+k](t) $k=-3..0])/h^2;

(-u[n-3](t)+4*u[n-2](t)-5*u[n-1](t)+2*u[n](t))/h^2

> 

# discretize u_xxxx, i=2,..,n-2
d4 := add([1,-4,6,-4,1] *~ [u[i+k](t) $k=-2..2])/h^4;

(u[i-2](t)-4*u[i-1](t)+6*u[i](t)-4*u[i+1](t)+u[i+2](t))/h^4

> 

# the discretized PDE at any i=2..n-2
eq__i := diff(u[i](t),t) = - d2 - d4 - u[i](t)*d1;  

diff(u[i](t), t) = -(u[i-1](t)-2*u[i](t)+u[i+1](t))/h^2-(u[i-2](t)-4*u[i-1](t)+6*u[i](t)-4*u[i+1](t)+u[i+2](t))/h^4-(1/2)*u[i](t)*(-u[i-1](t)+u[i+1](t))/h

Boundary conditions at the left end are u=0 and u_xx=0.
These determine u[0] and u[1] in terms of u[2] and u[3].

> 

u[0](t)=0, d2_left = 0:
bc_left := solve({%}, {u[0](t), u[1](t)});

{u[0](t) = 0, u[1](t) = (4/5)*u[2](t)-(1/5)*u[3](t)}

Boundary conditions at the right end are u=0 and u_xx=0.
These determine u[n-1] and u[n] in terms of u[n-2] and u[n-3].

> 

u[n](t)=0, d2_right = 0:
bc_right := solve({%}, {u[n-1](t), u[n](t)});

{u[n](t) = 0, u[n-1](t) = -(1/5)*u[n-3](t)+(4/5)*u[n-2](t)}

We pick a reasonable n and generate the system of ODEs

> 

n := 15:
L := 1:
h := L/n:
T := 0.02:   # solution will be computed on the time interval 0..T
sys := { seq(eq__i, i=2..n-2) } union bc_left union bc_right:

and here are the corresponding initial conditions

> 

IC := {seq(u[i](0)=eval(rhs(ic), x=i*h), i=2..n-2)}:

At this point, the PDE has been converted into a system of n+1 ODEs

in n+1 unknowns.  (That's the purpose of the Methods of Lines.)

We call Maple's dsolve to solve that system.

> 

dsol := dsolve(sys union IC, numeric,
    output=operator, range=0..T, maxfun=0):

Here is the plot of u__5(t)

> 

plots:-odeplot(dsol, [t, u[5](t)]);

This proc plots the solution u(x, t) produced by Method of Lines:

> 

my_plot3d := proc()
  local i, j, V, m:=40, k:= T/m;
  V := Array(0..n, 0..m, datatype=float[8]):
  for i from 0 to m do
    V[0..n,i] := Array([eval('u[j](i*k)' $j=0..n, dsol)]);
  end do;
  plots:-display(GRID(0..L, 0..T, V), labels=['x','t','u'], _rest);
end proc:

Here is what the solution looks like.  Compare to that produced by Maple's pdsolve:

> 

my_plot3d();

 

Download method-of-lines-illustration.mw

 

Another solution...

Rouben Rostamian 7711
May 24 2020
2 2

Here's an alternative to what Carl wrote.  This one avoids implicitplot, and therefore produces crisper graphics.

> 

restart;
> 

with(plots):
> 

with(plottools):

Let's see what the overall picture looks like

> 

display([
  sphere([0,0,0], 4, color=red),
  sphere([-4,0,2], 7, color=blue)
], scaling=constrained, transparency=0.6, labels=[x,y,z], orientation=[180,90,-90]);

Here is a 2D cross sectional diagram

> 

display([
  circle([ 0,0], 4, color=red),
  circle([-4,2], 7, color=blue),
  line([0,0], [-4,2], color="Green", thickness=3)
], scaling=constrained);

The green line segment in the diagram above connects the centers of the

circles and makes an angle of sigma with respect to the vertical axis:

> 

sigma := arctan(4/2);

arctan(2)

If we rotate the above diagram by an angle sigmaso that the green line segment

becomes vertical, the diagram simplifies for subsequent calculations:

> 

display([
  circle([0,0], 4, color=red),
  circle([0,4], 7, color=blue),
  line([0,0], [0, sqrt(2^2+4^2)], color="Green", thickness=3),
  line([0,0], [7*sqrt(15)/8, -17/8], color=red, thickness=3),
  line([0,4], [7*sqrt(15)/8, -17/8], color=blue, thickness=3)
], scaling=constrained, labels=[x,r]);

Write down the equations of the two circles, and then find their intersection points:

> 

x^2+r^2=16, x^2+(r-4)^2=49;
rx := allvalues(solve({%}));

r^2+x^2 = 16, x^2+(r-4)^2 = 49

{r = -17/8, x = (7/8)*15^(1/2)}, {r = -17/8, x = -(7/8)*15^(1/2)}

The angles of the red and blue line segments in the diagram above, relative

to the vertical axis:

> 

alpha := Pi/2 + arctan(  17/8,7*sqrt(15)/8);
beta :=  Pi/2 + arctan(4+17/8,7*sqrt(15)/8);

(1/2)*Pi+arctan((17/105)*15^(1/2))

(1/2)*Pi+arctan((7/15)*15^(1/2))

Now we plot the object in the spherical coordinates:

> 

p1 := display([
  plot3d([4*sin(theta)*cos(phi), 4*sin(theta)*sin(phi), 4*cos(theta)], theta=0..alpha, phi=0..2*Pi),
  plot3d([7*sin(theta)*cos(phi), 7*sin(theta)*sin(phi), 4+7*cos(theta)], theta=beta..Pi, phi=0..2*Pi)
], color=["Red","Green"]);

Restore the original orientation by rotating about the y axis by sigma

> 

plottools:-rotate(p1, sigma, [[0,0,0],[0,1,0]]):
display(%, labels=[x,y,z], scaling=constrained, style=surface, orientation=[-125,75,0]);

> 

 

 

Download mw.mw

Testing for a solution...

Rouben Rostamian 7711
May 22 2020
2 0

Here is a slightly modified and extended version of your worksheet.   The edited parts are in green.
 

> 

restart:
> 

with(PDEtools):
> 

with(LinearAlgebra):
> 

alias(f=f(x,t),g=g(x,t));

f, g

> 

eq1:=diff(f,x)=-I*eta*f +I*exp(-I*t)*g;

diff(f, x) = -I*eta*f+I*exp(-I*t)*g

> 

eq2:=diff(g,x)=-I*eta*g +I*exp(I*t)*f;

diff(g, x) = -I*eta*g+I*exp(I*t)*f

> 

eq3:=diff(f,t)=(I*eta^2-I/2)*f +I*eta*exp(-I*t)*g;

diff(f, t) = (I*eta^2-(1/2)*I)*f+I*eta*exp(-I*t)*g

> 

eq4:=diff(g,t)=(-I*eta^2+I/2)*g +I*eta*exp(I*t)*f;

diff(g, t) = (-I*eta^2+(1/2)*I)*g+I*eta*exp(I*t)*f

> 

#### The solution of (2)-(5) is
> 

#eq5:=f=I*(c1*exp(A)-c2*exp(-A))*exp(-i*t/2);
 eq5:=f=I*(c1*exp(A)-c2*exp(-A))*exp(-I*t/2);  # changed i to I  

f = I*(c1*exp(A)-c2*exp(-A))*exp(-((1/2)*I)*t)

> 

#eq6:=g=(c2*exp(A)-c1*exp(-A))*exp(i*t/2);
 eq6:=g=(c2*exp(A)-c1*exp(-A))*exp(I*t/2);     # changed i to I

g = (c2*exp(A)-c1*exp(-A))*exp(((1/2)*I)*t)

> 

#### where
> 

#c1=sqrt(h-sqrt(h^2-1))/sqrt(h^2-1);c2=sqrt(h+sqrt(h^2-1))/sqrt(h^2-1);A=sqrt(h^2-1)*(x+I*h*t);
> 

# replaced the three semicolons with commas:
> 

vals := c1=sqrt(h-sqrt(h^2-1))/sqrt(h^2-1), c2=sqrt(h+sqrt(h^2-1))/sqrt(h^2-1), A=sqrt(h^2-1)*(x+I*h*t);

c1 = (h-(h^2-1)^(1/2))^(1/2)/(h^2-1)^(1/2), c2 = (h+(h^2-1)^(1/2))^(1/2)/(h^2-1)^(1/2), A = (h^2-1)^(1/2)*(x+I*h*t)

> 

#### How to verify (6) and (7) is the solution of (2)-(5)?
> 

z1 := eval({eq5, eq6}, {vals}):   # plug in the values of c1, c2, A
> 

# Test the proposed solution.  If the test is successful,

# we will get a list of zeros.
> 

z2 := pdetest(z1, [eq1,eq2,eq3,eq4]);

[exp(-((1/2)*I)*t-(h^2-1)^(1/2)*x)*(I*exp(2*(h^2-1)^(1/2)*x+I*(h^2-1)^(1/2)*h*t)*(h-(h^2-1)^(1/2))^(1/2)*(h^2-1)^(1/2)-I*exp(2*(h^2-1)^(1/2)*x+I*(h^2-1)^(1/2)*h*t)*(h+(h^2-1)^(1/2))^(1/2)-exp(2*(h^2-1)^(1/2)*x+I*(h^2-1)^(1/2)*h*t)*(h-(h^2-1)^(1/2))^(1/2)*eta+I*exp(-I*(h^2-1)^(1/2)*h*t)*(h+(h^2-1)^(1/2))^(1/2)*(h^2-1)^(1/2)+I*exp(-I*(h^2-1)^(1/2)*h*t)*(h-(h^2-1)^(1/2))^(1/2)+exp(-I*(h^2-1)^(1/2)*h*t)*(h+(h^2-1)^(1/2))^(1/2)*eta)/(h^2-1)^(1/2), exp(((1/2)*I)*t-(h^2-1)^(1/2)*x)*(I*exp(2*(h^2-1)^(1/2)*x+I*(h^2-1)^(1/2)*h*t)*(h+(h^2-1)^(1/2))^(1/2)*eta+exp(2*(h^2-1)^(1/2)*x+I*(h^2-1)^(1/2)*h*t)*(h+(h^2-1)^(1/2))^(1/2)*(h^2-1)^(1/2)-I*exp(-I*(h^2-1)^(1/2)*h*t)*(h-(h^2-1)^(1/2))^(1/2)*eta+exp(2*(h^2-1)^(1/2)*x+I*(h^2-1)^(1/2)*h*t)*(h-(h^2-1)^(1/2))^(1/2)+exp(-I*(h^2-1)^(1/2)*h*t)*(h-(h^2-1)^(1/2))^(1/2)*(h^2-1)^(1/2)-exp(-I*(h^2-1)^(1/2)*h*t)*(h+(h^2-1)^(1/2))^(1/2))/(h^2-1)^(1/2), -exp(-((1/2)*I)*t-(h^2-1)^(1/2)*x)*(I*exp(2*(h^2-1)^(1/2)*x+I*(h^2-1)^(1/2)*h*t)*(h+(h^2-1)^(1/2))^(1/2)*eta+exp(2*(h^2-1)^(1/2)*x+I*(h^2-1)^(1/2)*h*t)*(h-(h^2-1)^(1/2))^(1/2)*(h^2-1)^(1/2)*h-exp(2*(h^2-1)^(1/2)*x+I*(h^2-1)^(1/2)*h*t)*(h-(h^2-1)^(1/2))^(1/2)*eta^2-I*exp(-I*(h^2-1)^(1/2)*h*t)*(h-(h^2-1)^(1/2))^(1/2)*eta+exp(-I*(h^2-1)^(1/2)*h*t)*(h+(h^2-1)^(1/2))^(1/2)*(h^2-1)^(1/2)*h+exp(-I*(h^2-1)^(1/2)*h*t)*(h+(h^2-1)^(1/2))^(1/2)*eta^2)/(h^2-1)^(1/2), exp(((1/2)*I)*t-(h^2-1)^(1/2)*x)*(I*exp(2*(h^2-1)^(1/2)*x+I*(h^2-1)^(1/2)*h*t)*(h+(h^2-1)^(1/2))^(1/2)*(h^2-1)^(1/2)*h+I*exp(2*(h^2-1)^(1/2)*x+I*(h^2-1)^(1/2)*h*t)*(h+(h^2-1)^(1/2))^(1/2)*eta^2+I*exp(-I*(h^2-1)^(1/2)*h*t)*(h-(h^2-1)^(1/2))^(1/2)*(h^2-1)^(1/2)*h-I*exp(-I*(h^2-1)^(1/2)*h*t)*(h-(h^2-1)^(1/2))^(1/2)*eta^2+exp(2*(h^2-1)^(1/2)*x+I*(h^2-1)^(1/2)*h*t)*(h-(h^2-1)^(1/2))^(1/2)*eta-exp(-I*(h^2-1)^(1/2)*h*t)*(h+(h^2-1)^(1/2))^(1/2)*eta)/(h^2-1)^(1/2)]

> 

# Those expressions are too complicated to see whether
# or not they are zeros.  We plug in numbers for x, t, eta,
# and h, to see whether at least for these numbers we get
# zeros.

> 

evalf(eval(z2, {x=0, t=0, eta=0, h=5}));

[2.886751346*I, 2.886751346, -17.32050808, 17.32050808*I]

> 

# We see that the result is not a list of zeros.  Therefore
# the proposed "solution" is not a solution.
> 

 

 

Download verification-2.mw

 

See if this helps...

Rouben Rostamian 7711
May 17 2020
1 0
> 

restart;
> 

eq := x^2/3^2 + y^2/4^2=1;

(1/9)*x^2+(1/16)*y^2 = 1

> 

x(t) = 1/sqrt(coeff(lhs(eq), x, 2))*cos(t);
y(t) = 1/sqrt(coeff(lhs(eq), y, 2))*sin(t);

x(t) = 3*cos(t)

y(t) = 4*sin(t)

 

The general solution...

Rouben Rostamian 7711
May 07 2020
2 0

Here is how to plot a general ellipsoid with semiaxis lengths a, b, c, centered at x0, y0, z0. The variables t and s are the polar and azimuthal angles in spherical coordinates.

restart;
a, b, c, x0, y0, z0 := 2, 3, 1, 10, 11, 12:
plot3d([x0 + a*sin(t)*cos(s), y0 + b*sin(t)*sin(s), z0 + c*cos(t)],
    t=0..Pi, s=0..2*Pi, scaling=constrained, labels=[x,y,z]);

Tenth degree polynomial...

Rouben Rostamian 7711
May 06 2020
1 0

We can read off the y values of the function at 9 points x=-4, -3, ... , 4, 5 (except for x=2) on the graph that you have provided.  To that collection I added estimated values of the function at -4.25 and 5.25, making for a total of eleven (x,y) pairs.  Then we fit the data with a 10th degree polynomial (which has 11 coefficients) and solve for the coefficients.

Note that the graph does not quite look like the given graph.  One reason for that is that in the given graph the coordinates are not spaced uniformly.  For instance, the point with coordinate −5 on the y axis falls where −3 would be.

That said, I agree with Kitonum's recommendation of using splines in such problems.  We are just lucky that a tenth degree polynomial provides a reasonable answer in this particular case.  In general, fitting a single polynomial to arbitrary data will not produce good results.

> 

restart;
> 

f := x -> add(a[i]*x^i, i=0..10);

proc (x) local i; options operator, arrow; add(a[i]*x^i, i = 0 .. 10) end proc

> 

eqns := f(-4.25)= 4, f(-4) = 3, f(-3) = 1, f(-2) = 1, f(-1) = 1,
        f(0) = 0, f(1) = -1, f(3) = -5, f(4) = -5, f(5) = 1, f(5.25)=3:
> 

sol := solve({eqns}):
> 

plot(eval(f(x), sol), x=-4..5.25, gridlines);

 

Download mw.mw

 

The roller coaster...

Rouben Rostamian 7711
May 05 2020
4 2

It's easier to apply Newton's equations directly rather than appeal to indirect quantities such as energy and centripetal force.  This worksheet shows how.

> 

restart;
> 

Typesetting:-Settings(typesetdot=true):

The position vector of the mass

> 

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

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

Velocity (which is tangent to the track):

> 

v := diff(r,t);

Vector(2, {(1) = diff(x(t), t), (2) = diff(y(t), t)})

Unit normal to the track:

> 

< -v[2], v[1] >:
n := % / sqrt(%^+ . %);

Vector(2, {(1) = -(diff(y(t), t))/sqrt((diff(y(t), t))^2+(diff(x(t), t))^2), (2) = (diff(x(t), t))/sqrt((diff(y(t), t))^2+(diff(x(t), t))^2)})

Resultant force acting on the mass, consisting of the downward weight mg 

and the (yet unknown) reaction force tau(t) in the direction n:

 

> 

F := < 0, -m*g> + tau(t)*n;

Vector(2, {(1) = -tau(t)*(diff(y(t), t))/sqrt((diff(y(t), t))^2+(diff(x(t), t))^2), (2) = -m*g+tau(t)*(diff(x(t), t))/sqrt((diff(y(t), t))^2+(diff(x(t), t))^2)})

We want the projection of the resultant force in the n direction be "2 m g," 

that is, F.n = 2*mg.  We solve that equation for tau(t), and simplify:

> 

F^+ . n = 2*m*g;
solve(%, tau(t)):
simplify(%, {diff(y(t),t)^2 + diff(x(t),t)^2 = V^2}):
simplify(%) assuming V>0:
expand(%):
eval(%, V=sqrt(diff(y(t),t)^2 + diff(x(t),t)^2)):
eq := tau(t) = %;

tau(t)*(diff(y(t), t))^2/((diff(y(t), t))^2+(diff(x(t), t))^2)+(-m*g+tau(t)*(diff(x(t), t))/((diff(y(t), t))^2+(diff(x(t), t))^2)^(1/2))*(diff(x(t), t))/((diff(y(t), t))^2+(diff(x(t), t))^2)^(1/2) = 2*m*g

tau(t) = m*g*(diff(x(t), t))/((diff(y(t), t))^2+(diff(x(t), t))^2)^(1/2)+2*m*g

Write down Newton's equation of motion, and then eliminate the

reaction force tau(t)from it by applying the result of the previous

computation:

> 

m*diff(r,t,t) =~ F;
eval(%, eq):
de := simplify(%);

Vector(2, {(1) = m*(diff(diff(x(t), t), t)) = -tau(t)*(diff(y(t), t))/sqrt((diff(y(t), t))^2+(diff(x(t), t))^2), (2) = m*(diff(diff(y(t), t), t)) = -m*g+tau(t)*(diff(x(t), t))/sqrt((diff(y(t), t))^2+(diff(x(t), t))^2)})

Vector[column](%id = 18446883779531074006)

That's a system of two second order differential equations in the

two unknowns x and y.  Here are some initial conditions to go with it.

Change as desired.

> 

ic := x(0)=0, y(0)=0, D(x)(0)=1, D(y)(0)=-1;

x(0) = 0, y(0) = 0, (D(x))(0) = 1, (D(y))(0) = -1

Numerical values of the parameters:

> 

params := g=1, m=1;

g = 1, m = 1

Apply the parameter values, and then solve the system numerically:

> 

eval(de, {params}):
dsol := dsolve({%[1], %[2], ic}, numeric);

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45) else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 26, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..63, {(1) = 4, (2) = 4, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 1, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = 0.5047658755841546e-2, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..4, {(1) = .0, (2) = 1.0, (3) = .0, (4) = -1.0}, datatype = float[8], order = C_order)), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..4, {(1) = .1, (2) = .1, (3) = .1, (4) = .1}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, 1..4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0}, datatype = float[8], order = C_order), Array(1..4, 1..4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, 1..4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0}, datatype = float[8], order = C_order), Array(1..4, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = 0, (2) = 0, (3) = 0, (4) = 0}, datatype = integer[8]), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..8, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = 0, (2) = 0, (3) = 0, (4) = 0}, datatype = integer[8])]), ( 8 ) = ([Array(1..4, {(1) = .0, (2) = 1.0, (3) = .0, (4) = -1.0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = 1.0, (2) = 1.9142135623730951, (3) = -1.0, (4) = .9142135623730951}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = x(t), Y[2] = diff(x(t),t), Y[3] = y(t), Y[4] = diff(y(t),t)]`; if Y[2]^2+Y[4]^2 < 0 then YP[1] := undefined; return 0 end if; YP[2] := -(2*evalf((Y[2]^2+Y[4]^2)^(1/2))+Y[2])*Y[4]/(Y[2]^2+Y[4]^2); YP[4] := (2*evalf((Y[2]^2+Y[4]^2)^(1/2))*Y[2]-Y[4]^2)/(Y[2]^2+Y[4]^2); YP[1] := Y[2]; YP[3] := Y[4]; 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = x(t), Y[2] = diff(x(t),t), Y[3] = y(t), Y[4] = diff(y(t),t)]`; if Y[2]^2+Y[4]^2 < 0 then YP[1] := undefined; return 0 end if; YP[2] := -(2*evalf((Y[2]^2+Y[4]^2)^(1/2))+Y[2])*Y[4]/(Y[2]^2+Y[4]^2); YP[4] := (2*evalf((Y[2]^2+Y[4]^2)^(1/2))*Y[2]-Y[4]^2)/(Y[2]^2+Y[4]^2); YP[1] := Y[2]; YP[3] := Y[4]; 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 26 ) = (Array(1..0, {})), ( 25 ) = (Array(1..0, {})), ( 24 ) = (0) ] )) ] ); _y0 := Array(0..4, {(1) = 0., (2) = 0., (3) = 1., (4) = 0.}); _vmap := array( 1 .. 4, [( 1 ) = (1), ( 2 ) = (2), ( 3 ) = (3), ( 4 ) = (4) ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); _i := false; if _par <> [] then _i := `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then _i := `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) or _i end if; if _i then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch: end try else try procname(procname("right")) catch: end try end if end if end if; return elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch: end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch: end try end if end if end if; return elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error end try; if _dat[17] <> _dtbl[1][17] then _dtbl[1][17] := _dat[17]; _dtbl[1][10] := _dat[10] end if; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; if type(_EnvDSNumericSaveDigits, 'posint') then _dat[4][26] := _EnvDSNumericSaveDigits else _dat[4][26] := Digits end if; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [t, x(t), diff(x(t), t), y(t), diff(y(t), t)], (4) = []}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol := 1; _ndsol := _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error end try end proc

The solution looks like a nice roller-coaster!

> 

plots:-odeplot(dsol, [x(t), y(t)], t=0..6, scaling=constrained);

 

 
 

Download mw.mw

 

Solution...

Rouben Rostamian 7711
May 02 2020
2 6

This looks like a homework problem.  Here I will provide most of the solution but in the future please show some evidence of what you can do on your own before appealing for help.

> 

restart;

The general  system of two autonomous differential equations has the form
"w ' = A(w,f)"
"f ' = B(w,f)"
where a prime indicates derivative with respect to x.  In your case your have

> 

A := (w,f) -> f;
B := (w,f) -> a*w + b*w^2 + c*w^3;

proc (w, f) options operator, arrow; f end proc

proc (w, f) options operator, arrow; a*w+b*w^2+c*w^3 end proc

and therefore your equations are

> 

de1 := diff(w(x),x) = A(w(x),f(x));

diff(w(x), x) = f(x)

> 

de2 := diff(f(x),x) = B(w(x),f(x));

diff(f(x), x) = a*w(x)+b*w(x)^2+c*w(x)^3

The system's equilibra are obtained by solving A(w, f) = 0, B(w, f) = 0:

> 

equilib := solve({A(w,f)=0, B(w,f)=0}, {w,f}, explicit);

{f = 0, w = 0}, {f = 0, w = (1/2)*(-b+(-4*a*c+b^2)^(1/2))/c}, {f = 0, w = -(1/2)*(b+(-4*a*c+b^2)^(1/2))/c}

Linearize the system about each of the equilibria.  I do it here for one, you do the others.

 

The Jacobian matrix is

> 

J := < diff(A(w,f),w),  diff(A(w,f),f);
       diff(B(w,f),w),  diff(B(w,f),f) >;

Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 3*c*w^2+2*b*w+a, (2, 2) = 0})

Evaluate the Jacobian at equilib[3]:

> 

J2 := simplify(eval(J, equilib[3]));

Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = (1/2)*(b*sqrt(-4*a*c+b^2)-4*a*c+b^2)/c, (2, 2) = 0})

Here are the eigenvalues at equilib[3]:

> 

eigs := LinearAlgebra:-Eigenvalues(J2);

Vector(2, {(1) = (1/2)*sqrt(2)*sqrt(c*(b*sqrt(-4*a*c+b^2)-4*a*c+b^2))/c, (2) = -(1/2)*sqrt(2)*sqrt(c*(b*sqrt(-4*a*c+b^2)-4*a*c+b^2))/c})

> 

params :=
    a = (r+1)*(1+p1*s)/(2*(1-p1))-1/(n^2),
    b = (r+1)*(3-r)*(1-p1*s^2)/(8*(1-p1))-3/(2*n^4),
    c = (r+1)*(3-r)*(5-3*r)*(1+p1*s^3)/(48*(1-p1))-5/(2*n^6),
    r=-0.8, p1=0.5, s=0.6, n=1.6;

a = (r+1)*(p1*s+1)/(2-2*p1)-1/n^2, b = (r+1)*(3-r)*(-p1*s^2+1)/(8-8*p1)-(3/2)/n^4, c = (r+1)*(3-r)*(5-3*r)*(p1*s^3+1)/(48-48*p1)-(5/2)/n^6, r = -.8, p1 = .5, s = .6, n = 1.6

> 

evalf(subs(params, eigs));

Vector(2, {(1) = .4498767329, (2) = -.4498767329})

One eigenvalue is positive, the other is negative, therefore equlib[3] is a saddle.

> 

 

Now let's do the phase portrait

> 

ic := seq([w(0)=h, f(0)=0], h=-2..2, 0.2),
      seq([w(0)=0, f(0)=h], h=0.5..1, 0.1),
      seq([w(0)=0, f(0)=h], h=-1..0.5, 0.1):
> 

subs(params, {de1, de2}):
DEtools:-DEplot(%, [w(x),f(x)], x=-10..10, [ic],
    w=-2..2, f=-1..1, linecolor=black, thickness=1,
    stepsize=0.1, arrows=none);

 

Download mw.mw

 

 

Comparison or solutions...

Rouben Rostamian 7711
May 01 2020
2 0
> 

restart;
> 

de := diff(a3(z), z) = -K*a10*a20*sqrt(1-a3(z)^2/a10^2)*sqrt(1-a3(z)^2/a20^2);

diff(a3(z), z) = -K*a10*a20*(1-a3(z)^2/a10^2)^(1/2)*(1-a3(z)^2/a20^2)^(1/2)

> 

ic := a3(0) = 0;

a3(0) = 0

Arbitrary choice of parameters

> 

params := a10=1.23, a20=2.34, K=4.23;

a10 = 1.23, a20 = 2.34, K = 4.23

Numeric soluition

> 

eval({de,ic}, {params});
dsol := dsolve(%, numeric, range=0..1);

{a3(0) = 0, diff(a3(z), z) = -12.174786*(1-.6609822196*a3(z)^2)^(1/2)*(1-.1826283878*a3(z)^2)^(1/2)}

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45) else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "left" ) = 0., ( "right" ) = 1., ( "complex" ) = false ] ) _xout := _xin; _pars := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 26, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..63, {(1) = 1, (2) = 1, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 1, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 1, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = 1.0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = 0.41459938234984545e-3, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..1, {(1) = .0}, datatype = float[8], order = C_order)), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..1, {(1) = .1}, datatype = float[8], order = C_order), Array(1..1, {(1) = -1.2299999999798035}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = -.0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..6, {(1, 1) = -.0, (1, 2) = -.0, (1, 3) = -.0, (1, 4) = -.0, (1, 5) = -.0, (1, 6) = -.0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8]), Array(1..1, {(1) = -1.2299999999798037}, datatype = float[8], order = C_order), Array(1..1, {(1) = -1.2299999999798035}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0.296010311817696e-7}, datatype = float[8], order = C_order), Array(1..1, {(1) = -1.2299999999798035}, datatype = float[8], order = C_order), Array(1..1, {(1) = -0.18901959074120248e-6}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8])]), ( 8 ) = ([Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = -12.174786}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..1, {(1, 1) = .0, (2, 0) = .0, (2, 1) = .0, (3, 0) = .0, (3, 1) = .0, (4, 0) = .0, (4, 1) = .0, (5, 0) = .0, (5, 1) = .0, (6, 0) = .0, (6, 1) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = a3(z)]`; if -.1826283878*Y[1]^2 < -1 then YP[1] := undefined; return 0 end if; if -.6609822196*Y[1]^2 < -1 then YP[1] := undefined; return 0 end if; YP[1] := -12.174786*evalf((1-.6609822196*Y[1]^2)^(1/2))*evalf((1-.1826283878*Y[1]^2)^(1/2)); 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = a3(z)]`; if -.1826283878*Y[1]^2 < -1 then YP[1] := undefined; return 0 end if; if -.6609822196*Y[1]^2 < -1 then YP[1] := undefined; return 0 end if; YP[1] := -12.174786*evalf((1-.6609822196*Y[1]^2)^(1/2))*evalf((1-.1826283878*Y[1]^2)^(1/2)); 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 26 ) = (Array(1..0, {})), ( 25 ) = (Array(1..0, {})), ( 24 ) = (0) ] )), ( 3 ) = (array( 1 .. 26, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..63, {(1) = 1, (2) = 1, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 1, (9) = 0, (10) = 1, (11) = 305, (12) = 305, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 662, (19) = 30000, (20) = 5, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = 1.0, (2) = 0.10e-5, (3) = .8581218124879464, (4) = 0.500001e-14, (5) = .0, (6) = 0.41459938234984545e-3, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..1, {(1) = .0}, datatype = float[8], order = C_order)), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..1, {(1) = .1}, datatype = float[8], order = C_order), Array(1..1, {(1) = -1.2299999999798035}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = -.0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..6, {(1, 1) = -.0, (1, 2) = -.0, (1, 3) = -.0, (1, 4) = -.0, (1, 5) = -.0, (1, 6) = -.0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8]), Array(1..1, {(1) = -1.2299999999798037}, datatype = float[8], order = C_order), Array(1..1, {(1) = -1.2299999999798035}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0.296010311817696e-7}, datatype = float[8], order = C_order), Array(1..1, {(1) = -1.2299999999798035}, datatype = float[8], order = C_order), Array(1..1, {(1) = -0.18901959074120248e-6}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8])]), ( 8 ) = ([Array(1..1, {(1) = -1.2299999999798035}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = -.0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..1, {(1, 1) = .5904539721647559, (2, 0) = .5904539721647559, (2, 1) = -1.2299999999798035, (3, 0) = -1.2299999999798035, (3, 1) = .7157363544252189, (4, 0) = .7157363544252189, (4, 1) = -1.2299999999798032, (5, 0) = -1.2299999999798032, (5, 1) = .8410187366856818, (6, 0) = .8410187366856818, (6, 1) = -1.2299999999798035}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = a3(z)]`; if -.1826283878*Y[1]^2 < -1 then YP[1] := undefined; return 0 end if; if -.6609822196*Y[1]^2 < -1 then YP[1] := undefined; return 0 end if; YP[1] := -12.174786*evalf((1-.6609822196*Y[1]^2)^(1/2))*evalf((1-.1826283878*Y[1]^2)^(1/2)); 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (Array(1..305, 0..1, {(1, 1) = .0, (2, 0) = .0, (2, 1) = .0, (3, 0) = .0, (3, 1) = 0.10364984558746136e-3, (4, 0) = 0.10364984558746136e-3, (4, 1) = -0.12619144064209943e-2, (5, 0) = -0.12619144064209943e-2, (5, 1) = 0.20729969117492272e-3, (6, 0) = 0.20729969117492272e-3, (6, 1) = -0.25238271175993233e-2, (7, 0) = -0.25238271175993233e-2, (7, 1) = 0.3109495367623841e-3, (8, 0) = 0.3109495367623841e-3, (8, 1) = -0.37857364383006736e-2, (9, 0) = -0.37857364383006736e-2, (9, 1) = 0.41459938234984545e-3, (10, 0) = 0.41459938234984545e-3, (10, 1) = -0.5047640673306076e-2, (11, 0) = -0.5047640673306076e-2, (11, 1) = 0.34833604325880127e-2, (12, 0) = 0.34833604325880127e-2, (12, 1) = -0.4239847049454727e-1, (13, 0) = -0.4239847049454727e-1, (13, 1) = 0.6552121482826179e-2, (14, 0) = 0.6552121482826179e-2, (14, 1) = -0.7969938711162143e-1, (15, 0) = -0.7969938711162143e-1, (15, 1) = 0.9620882533064347e-2, (16, 0) = 0.9620882533064347e-2, (16, 1) = -.11690662945170899, (17, 0) = -.11690662945170899, (17, 1) = 0.12689643583302514e-1, (18, 0) = 0.12689643583302514e-1, (18, 1) = -.15397680385648438, (19, 0) = -.15397680385648438, (19, 1) = 0.16122517299749697e-1, (20, 0) = 0.16122517299749697e-1, (20, 1) = -.19523010696600793, (21, 0) = -.19523010696600793, (21, 1) = 0.1955539101619688e-1, (22, 0) = 0.1955539101619688e-1, (22, 1) = -.2361989204547763, (23, 0) = -.2361989204547763, (23, 1) = 0.22988264732644058e-1, (24, 0) = 0.22988264732644058e-1, (24, 1) = -.2768253183563231, (25, 0) = -.2768253183563231, (25, 1) = 0.2642113844909124e-1, (26, 0) = 0.2642113844909124e-1, (26, 1) = -.3170529026160392, (27, 0) = -.3170529026160392, (27, 1) = 0.30427773082452543e-1, (28, 0) = 0.30427773082452543e-1, (28, 1) = -.36342665103366517, (29, 0) = -.36342665103366517, (29, 1) = 0.3443440771581385e-1, (30, 0) = 0.3443440771581385e-1, (30, 1) = -.4090987177000091, (31, 0) = -.4090987177000091, (31, 1) = 0.38441042349175156e-1, (32, 0) = 0.38441042349175156e-1, (32, 1) = -.4539892102217507, (33, 0) = -.4539892102217507, (33, 1) = 0.4244767698253646e-1, (34, 0) = 0.4244767698253646e-1, (34, 1) = -.4980225521011199, (35, 0) = -.4980225521011199, (35, 1) = 0.4707601510214095e-1, (36, 0) = 0.4707601510214095e-1, (36, 1) = -.5477285925071234, (37, 0) = -.5477285925071234, (37, 1) = 0.51704353221745436e-1, (38, 0) = 0.51704353221745436e-1, (38, 1) = -.5960941342765066, (39, 0) = -.5960941342765066, (39, 1) = 0.5633269134134993e-1, (40, 0) = 0.5633269134134993e-1, (40, 1) = -.6430260154216566, (41, 0) = -.6430260154216566, (41, 1) = 0.6096102946095442e-1, (42, 0) = 0.6096102946095442e-1, (42, 1) = -.6884400249508265, (43, 0) = -.6884400249508265, (43, 1) = 0.6591253046747747e-1, (44, 0) = 0.6591253046747747e-1, (44, 1) = -.7352591969666599, (45, 0) = -.7352591969666599, (45, 1) = 0.708640314740005e-1, (46, 0) = 0.708640314740005e-1, (46, 1) = -.7801740469048234, (47, 0) = -.7801740469048234, (47, 1) = 0.7581553248052356e-1, (48, 0) = 0.7581553248052356e-1, (48, 1) = -.823114685686209, (49, 0) = -.823114685686209, (49, 1) = 0.8076703348704659e-1, (50, 0) = 0.8076703348704659e-1, (50, 1) = -.8640218768230735, (51, 0) = -.8640218768230735, (51, 1) = 0.8484972859805444e-1, (52, 0) = 0.8484972859805444e-1, (52, 1) = -.8961869115432886, (53, 0) = -.8961869115432886, (53, 1) = 0.8893242370906229e-1, (54, 0) = 0.8893242370906229e-1, (54, 1) = -.9269138185229375, (55, 0) = -.9269138185229375, (55, 1) = 0.9301511882007014e-1, (56, 0) = 0.9301511882007014e-1, (56, 1) = -.9561842211393343, (57, 0) = -.9561842211393343, (57, 1) = 0.97097813931078e-1, (58, 0) = 0.97097813931078e-1, (58, 1) = -.9839833318110524, (59, 0) = -.9839833318110524, (59, 1) = .10062211222744795, (60, 0) = .10062211222744795, (60, 1) = -1.00678840597731, (61, 0) = -1.00678840597731, (61, 1) = .1041464105238179, (62, 0) = .1041464105238179, (62, 1) = -1.0284836142185316, (63, 0) = -1.0284836142185316, (63, 1) = .10767070882018787, (64, 0) = .10767070882018787, (64, 1) = -1.0490652330689563, (65, 0) = -1.0490652330689563, (65, 1) = .11119500711655782, (66, 0) = .11119500711655782, (66, 1) = -1.0685307214599358, (67, 0) = -1.0685307214599358, (67, 1) = .11431310243096295, (68, 0) = .11431310243096295, (68, 1) = -1.0848210444115174, (69, 0) = -1.0848210444115174, (69, 1) = .11743119774536806, (70, 0) = .11743119774536806, (70, 1) = -1.1002368207021282, (71, 0) = -1.1002368207021282, (71, 1) = .12054929305977319, (72, 0) = .12054929305977319, (72, 1) = -1.1147785546137297, (73, 0) = -1.1147785546137297, (73, 1) = .12366738837417832, (74, 0) = .12366738837417832, (74, 1) = -1.1284469811889712, (75, 0) = -1.1284469811889712, (75, 1) = .12640348363182397, (76, 0) = .12640348363182397, (76, 1) = -1.139722414348455, (77, 0) = -1.139722414348455, (77, 1) = .12913957888946964, (78, 0) = .12913957888946964, (78, 1) = -1.1503276739740989, (79, 0) = -1.1503276739740989, (79, 1) = .1318756741471153, (80, 0) = .1318756741471153, (80, 1) = -1.1602640698603375, (81, 0) = -1.1602640698603375, (81, 1) = .13461176940476094, (82, 0) = .13461176940476094, (82, 1) = -1.1695327191271392, (83, 0) = -1.1695327191271392, (83, 1) = .13696983833440174, (84, 0) = .13696983833440174, (84, 1) = -1.176985980079806, (85, 0) = -1.176985980079806, (85, 1) = .13932790726404257, (86, 0) = .13932790726404257, (86, 1) = -1.1839454443423434, (87, 0) = -1.1839454443423434, (87, 1) = .1416859761936834, (88, 0) = .1416859761936834, (88, 1) = -1.1904121474004363, (89, 0) = -1.1904121474004363, (89, 1) = .1440440451233242, (90, 0) = .1440440451233242, (90, 1) = -1.1963867324416086, (91, 0) = -1.1963867324416086, (91, 1) = .1460273668718157, (92, 0) = .1460273668718157, (92, 1) = -1.2010312720960516, (93, 0) = -1.2010312720960516, (93, 1) = .1480106886203072, (94, 0) = .1480106886203072, (94, 1) = -1.205329097234056, (95, 0) = -1.205329097234056, (95, 1) = .14999401036879872, (96, 0) = .14999401036879872, (96, 1) = -1.2092808344125783, (97, 0) = -1.2092808344125783, (97, 1) = .15197733211729023, (98, 0) = .15197733211729023, (98, 1) = -1.21288657506271, (99, 0) = -1.21288657506271, (99, 1) = .1532419292906772, (100, 0) = .1532419292906772, (100, 1) = -1.2150050645441022, (101, 0) = -1.2150050645441022, (101, 1) = .15450652646406415, (102, 0) = .15450652646406415, (102, 1) = -1.216983271248191, (103, 0) = -1.216983271248191, (103, 1) = .1557711236374511, (104, 0) = .1557711236374511, (104, 1) = -1.2188213276236481, (105, 0) = -1.2188213276236481, (105, 1) = .15703572081083808, (106, 0) = .15703572081083808, (106, 1) = -1.2205192473023838, (107, 0) = -1.2205192473023838, (107, 1) = .15830031798422506, (108, 0) = .15830031798422506, (108, 1) = -1.2220769857946303, (109, 0) = -1.2220769857946303, (109, 1) = .159564915157612, (110, 0) = .159564915157612, (110, 1) = -1.223494984245328, (111, 0) = -1.223494984245328, (111, 1) = .16082951233099896, (112, 0) = .16082951233099896, (112, 1) = -1.2247734508793204, (113, 0) = -1.2247734508793204, (113, 1) = .16209410950438594, (114, 0) = .16209410950438594, (114, 1) = -1.2259120009260447, (115, 0) = -1.2259120009260447, (115, 1) = .1628087356247067, (116, 0) = .1628087356247067, (116, 1) = -1.2264933895115786, (117, 0) = -1.2264933895115786, (117, 1) = .16352336174502746, (118, 0) = .16352336174502746, (118, 1) = -1.227030198914469, (119, 0) = -1.227030198914469, (119, 1) = .16423798786534824, (120, 0) = .16423798786534824, (120, 1) = -1.227522456457831, (121, 0) = -1.227522456457831, (121, 1) = .16495261398566902, (122, 0) = .16495261398566902, (122, 1) = -1.2279701111889383, (123, 0) = -1.2279701111889383, (123, 1) = .1656672401059898, (124, 0) = .1656672401059898, (124, 1) = -1.228373025983339, (125, 0) = -1.228373025983339, (125, 1) = .16638186622631057, (126, 0) = .16638186622631057, (126, 1) = -1.228731595465733, (127, 0) = -1.228731595465733, (127, 1) = .16709649234663132, (128, 0) = .16709649234663132, (128, 1) = -1.2290459781678713, (129, 0) = -1.2290459781678713, (129, 1) = .1678111184669521, (130, 0) = .1678111184669521, (130, 1) = -1.2293156267116803, (131, 0) = -1.2293156267116803, (131, 1) = .1679534238328131, (132, 0) = .1679534238328131, (132, 1) = -1.2293639137091485, (133, 0) = -1.2293639137091485, (133, 1) = .16809572919867405, (134, 0) = .16809572919867405, (134, 1) = -1.2294104341378878, (135, 0) = -1.2294104341378878, (135, 1) = .168238034564535, (136, 0) = .168238034564535, (136, 1) = -1.2294551880605513, (137, 0) = -1.2294551880605513, (137, 1) = .168380339930396, (138, 0) = .168380339930396, (138, 1) = -1.2294981754390624, (139, 0) = -1.2294981754390624, (139, 1) = .168522645296257, (140, 0) = .168522645296257, (140, 1) = -1.2295393962288792, (141, 0) = -1.2295393962288792, (141, 1) = .16866495066211795, (142, 0) = .16866495066211795, (142, 1) = -1.2295788506121428, (143, 0) = -1.2295788506121428, (143, 1) = .1688072560279789, (144, 0) = .1688072560279789, (144, 1) = -1.2296165386693825, (145, 0) = -1.2296165386693825, (145, 1) = .1689495613938399, (146, 0) = .1689495613938399, (146, 1) = -1.229652460286641, (147, 0) = -1.229652460286641, (147, 1) = .16907526225452904, (148, 0) = .16907526225452904, (148, 1) = -1.2296827211078099, (149, 0) = -1.2296827211078099, (149, 1) = .16920096311521818, (150, 0) = .16920096311521818, (150, 1) = -1.2297116037888527, (151, 0) = -1.2297116037888527, (151, 1) = .16932666397590732, (152, 0) = .16932666397590732, (152, 1) = -1.2297391083972988, (153, 0) = -1.2297391083972988, (153, 1) = .16945236483659645, (154, 0) = .16945236483659645, (154, 1) = -1.2297652348036772, (155, 0) = -1.2297652348036772, (155, 1) = .1695780656972856, (156, 0) = .1695780656972856, (156, 1) = -1.2297899828602528, (157, 0) = -1.2297899828602528, (157, 1) = .16970376655797473, (158, 0) = .16970376655797473, (158, 1) = -1.2298133529350264, (159, 0) = -1.2298133529350264, (159, 1) = .16982946741866387, (160, 0) = .16982946741866387, (160, 1) = -1.2298353451592163, (161, 0) = -1.2298353451592163, (161, 1) = .169955168279353, (162, 0) = .169955168279353, (162, 1) = -1.2298559591771572, (163, 0) = -1.2298559591771572, (163, 1) = .17005074170152765, (164, 0) = .17005074170152765, (164, 1) = -1.2298707100162487, (165, 0) = -1.2298707100162487, (165, 1) = .1701463151237023, (166, 0) = .1701463151237023, (166, 1) = -1.2298846642688712, (167, 0) = -1.2298846642688712, (167, 1) = .17024188854587693, (168, 0) = .17024188854587693, (168, 1) = -1.2298978219992245, (169, 0) = -1.2298978219992245, (169, 1) = .17033746196805158, (170, 0) = .17033746196805158, (170, 1) = -1.2299101830274517, (171, 0) = -1.2299101830274517, (171, 1) = .17043303539022622, (172, 0) = .17043303539022622, (172, 1) = -1.2299217471301178, (173, 0) = -1.2299217471301178, (173, 1) = .17052860881240087, (174, 0) = .17052860881240087, (174, 1) = -1.2299325148446703, (175, 0) = -1.2299325148446703, (175, 1) = .1706241822345755, (176, 0) = .1706241822345755, (176, 1) = -1.2299424863556812, (177, 0) = -1.2299424863556812, (177, 1) = .17071975565675016, (178, 0) = .17071975565675016, (178, 1) = -1.22995166105592, (179, 0) = -1.22995166105592, (179, 1) = .1707878199736112, (180, 0) = .1707878199736112, (180, 1) = -1.229957708889341, (181, 0) = -1.229957708889341, (181, 1) = .17085588429047222, (182, 0) = .17085588429047222, (182, 1) = -1.2299633527954243, (183, 0) = -1.2299633527954243, (183, 1) = .17092394860733326, (184, 0) = .17092394860733326, (184, 1) = -1.229968592854273, (185, 0) = -1.229968592854273, (185, 1) = .1709920129241943, (186, 0) = .1709920129241943, (186, 1) = -1.2299734288021393, (187, 0) = -1.2299734288021393, (187, 1) = .17106007724105535, (188, 0) = .17106007724105535, (188, 1) = -1.229977860204638, (189, 0) = -1.229977860204638, (189, 1) = .1711281415579164, (190, 0) = .1711281415579164, (190, 1) = -1.2299818881769908, (191, 0) = -1.2299818881769908, (191, 1) = .1711962058747774, (192, 0) = .1711962058747774, (192, 1) = -1.2299855131258797, (193, 0) = -1.2299855131258797, (193, 1) = .17126427019163845, (194, 0) = .17126427019163845, (194, 1) = -1.2299887335965314, (195, 0) = -1.2299887335965314, (195, 1) = .17130712356026062, (196, 0) = .17130712356026062, (196, 1) = -1.2299905531051283, (197, 0) = -1.2299905531051283, (197, 1) = .1713499769288828, (198, 0) = .1713499769288828, (198, 1) = -1.2299922126515708, (199, 0) = -1.2299922126515708, (199, 1) = .17139283029750496, (200, 0) = .17139283029750496, (200, 1) = -1.2299937123674143, (201, 0) = -1.2299937123674143, (201, 1) = .17143568366612713, (202, 0) = .17143568366612713, (202, 1) = -1.229995051783982, (203, 0) = -1.229995051783982, (203, 1) = .1714785370347493, (204, 0) = .1714785370347493, (204, 1) = -1.2299962284419177, (205, 0) = -1.2299962284419177, (205, 1) = .17152139040337147, (206, 0) = .17152139040337147, (206, 1) = -1.229997248910503, (207, 0) = -1.229997248910503, (207, 1) = .17156424377199364, (208, 0) = .17156424377199364, (208, 1) = -1.2299981160027433, (209, 0) = -1.2299981160027433, (209, 1) = .1716070971406158, (210, 0) = .1716070971406158, (210, 1) = -1.2299988194501785, (211, 0) = -1.2299988194501785, (211, 1) = .1716268701874011, (212, 0) = .1716268701874011, (212, 1) = -1.2299990858657959, (213, 0) = -1.2299990858657959, (213, 1) = .17164664323418638, (214, 0) = .17164664323418638, (214, 1) = -1.2299993186848845, (215, 0) = -1.2299993186848845, (215, 1) = .17166641628097168, (216, 0) = .17166641628097168, (216, 1) = -1.2299995182568093, (217, 0) = -1.2299995182568093, (217, 1) = .17168618932775698, (218, 0) = .17168618932775698, (218, 1) = -1.2299996833052242, (219, 0) = -1.2299996833052242, (219, 1) = .1716909618897975, (220, 0) = .1716909618897975, (220, 1) = -1.2299997177818978, (221, 0) = -1.2299997177818978, (221, 1) = .17169573445183803, (222, 0) = .17169573445183803, (222, 1) = -1.2299997502723743, (223, 0) = -1.2299997502723743, (223, 1) = .17170050701387857, (224, 0) = .17170050701387857, (224, 1) = -1.2299997807768501, (225, 0) = -1.2299997807768501, (225, 1) = .17170527957591908, (226, 0) = .17170527957591908, (226, 1) = -1.2299998092947007, (227, 0) = -1.2299998092947007, (227, 1) = .1717100521379596, (228, 0) = .1717100521379596, (228, 1) = -1.2299998358251032, (229, 0) = -1.2299998358251032, (229, 1) = .17171482470000013, (230, 0) = .17171482470000013, (230, 1) = -1.2299998603700404, (231, 0) = -1.2299998603700404, (231, 1) = .17171959726204067, (232, 0) = .17171959726204067, (232, 1) = -1.2299998829301964, (233, 0) = -1.2299998829301964, (233, 1) = .17172436982408118, (234, 0) = .17172436982408118, (234, 1) = -1.2299999035032116, (235, 0) = -1.2299999035032116, (235, 1) = .17172897758102992, (236, 0) = .17172897758102992, (236, 1) = -1.2299999214762944, (237, 0) = -1.2299999214762944, (237, 1) = .17173358533797867, (238, 0) = .17173358533797867, (238, 1) = -1.2299999376035862, (239, 0) = -1.2299999376035862, (239, 1) = .1717381930949274, (240, 0) = .1717381930949274, (240, 1) = -1.2299999518890175, (241, 0) = -1.2299999518890175, (241, 1) = .17174280085187615, (242, 0) = .17174280085187615, (242, 1) = -1.2299999643183166, (243, 0) = -1.2299999643183166, (243, 1) = .17174391301259395, (244, 0) = .17174391301259395, (244, 1) = -1.229999967038151, (245, 0) = -1.229999967038151, (245, 1) = .17174502517331172, (246, 0) = .17174502517331172, (246, 1) = -1.2299999696501147, (247, 0) = -1.2299999696501147, (247, 1) = .1717461373340295, (248, 0) = .1717461373340295, (248, 1) = -1.2299999721542105, (249, 0) = -1.2299999721542105, (249, 1) = .1717472494947473, (250, 0) = .1717472494947473, (250, 1) = -1.2299999745504322, (251, 0) = -1.2299999745504322, (251, 1) = .1717483616554651, (252, 0) = .1717483616554651, (252, 1) = -1.229999976838773, (253, 0) = -1.229999976838773, (253, 1) = .17174947381618286, (254, 0) = .17174947381618286, (254, 1) = -1.2299999790192464, (255, 0) = -1.2299999790192464, (255, 1) = .17175058597690063, (256, 0) = .17175058597690063, (256, 1) = -1.2299999810918574, (257, 0) = -1.2299999810918574, (257, 1) = .17175169813761842, (258, 0) = .17175169813761842, (258, 1) = -1.229999983056592, (259, 0) = -1.229999983056592, (259, 1) = .17175277189344013, (260, 0) = .17175277189344013, (260, 1) = -1.2299999848511147, (261, 0) = -1.2299999848511147, (261, 1) = .17175384564926183, (262, 0) = .17175384564926183, (262, 1) = -1.2299999865450983, (263, 0) = -1.2299999865450983, (263, 1) = .17175491940508356, (264, 0) = .17175491940508356, (264, 1) = -1.2299999881385508, (265, 0) = -1.2299999881385508, (265, 1) = .17175599316090526, (266, 0) = .17175599316090526, (266, 1) = -1.2299999896314446, (267, 0) = -1.2299999896314446, (267, 1) = .17175706691672696, (268, 0) = .17175706691672696, (268, 1) = -1.2299999910237447, (269, 0) = -1.2299999910237447, (269, 1) = .1717581406725487, (270, 0) = .1717581406725487, (270, 1) = -1.229999992315535, (271, 0) = -1.229999992315535, (271, 1) = .1717592144283704, (272, 0) = .1717592144283704, (272, 1) = -1.229999993506844, (273, 0) = -1.229999993506844, (273, 1) = .1717602881841921, (274, 0) = .1717602881841921, (274, 1) = -1.2299999945975741, (275, 0) = -1.2299999945975741, (275, 1) = .17176132486130738, (276, 0) = .17176132486130738, (276, 1) = -1.2299999955550536, (277, 0) = -1.2299999955550536, (277, 1) = .17176236153842267, (278, 0) = .17176236153842267, (278, 1) = -1.2299999964190211, (279, 0) = -1.2299999964190211, (279, 1) = .17176339821553796, (280, 0) = .17176339821553796, (280, 1) = -1.229999997189621, (281, 0) = -1.229999997189621, (281, 1) = .17176443489265325, (282, 0) = .17176443489265325, (282, 1) = -1.229999997866332, (283, 0) = -1.229999997866332, (283, 1) = .17176468511235396, (284, 0) = .17176468511235396, (284, 1) = -1.2299999980155247, (285, 0) = -1.2299999980155247, (285, 1) = .17176493533205467, (286, 0) = .17176493533205467, (286, 1) = -1.229999998159257, (287, 0) = -1.229999998159257, (287, 1) = .17176518555175538, (288, 0) = .17176518555175538, (288, 1) = -1.2299999982975294, (289, 0) = -1.2299999982975294, (289, 1) = .1717654357714561, (290, 0) = .1717654357714561, (290, 1) = -1.2299999984303411, (291, 0) = -1.2299999984303411, (291, 1) = .1717656859911568, (292, 0) = .1717656859911568, (292, 1) = -1.229999998557692, (293, 0) = -1.229999998557692, (293, 1) = .1717659362108575, (294, 0) = .1717659362108575, (294, 1) = -1.2299999986795835, (295, 0) = -1.2299999986795835, (295, 1) = .17176618643055822, (296, 0) = .17176618643055822, (296, 1) = -1.2299999987960144, (297, 0) = -1.2299999987960144, (297, 1) = .17176643665025892, (298, 0) = .17176643665025892, (298, 1) = -1.229999998906985, (299, 0) = -1.229999998906985, (299, 1) = .17176667822942643, (300, 0) = .17176667822942643, (300, 1) = -1.229999999008942, (301, 0) = -1.229999999008942, (301, 1) = .17176691980859393, (302, 0) = .17176691980859393, (302, 1) = -1.22999999910581, (303, 0) = -1.22999999910581, (303, 1) = .17176716138776146, (304, 0) = .17176716138776146, (304, 1) = -1.229999999197589, (305, 0) = -1.229999999197589, (305, 1) = .17176740296692897}, datatype = float[8], order = C_order)), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = a3(z)]`; if -.1826283878*Y[1]^2 < -1 then YP[1] := undefined; return 0 end if; if -.6609822196*Y[1]^2 < -1 then YP[1] := undefined; return 0 end if; YP[1] := -12.174786*evalf((1-.6609822196*Y[1]^2)^(1/2))*evalf((1-.1826283878*Y[1]^2)^(1/2)); 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 26 ) = (Array(1..0, {})), ( 25 ) = (Array(1..0, {})), ( 24 ) = (0) ] )), ( 4 ) = (3) ] ); _y0 := Array(0..1, {(1) = 0.}); _vmap := array( 1 .. 1, [( 1 ) = (1) ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); _i := false; if _par <> [] then _i := `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then _i := `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) or _i end if; if _i then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch: end try else try procname(procname("right")) catch: end try end if end if end if; return elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch: end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch: end try end if end if end if; return elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error end try; if _dat[17] <> _dtbl[1][17] then _dtbl[1][17] := _dat[17]; _dtbl[1][10] := _dat[10] end if; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; if type(_EnvDSNumericSaveDigits, 'posint') then _dat[4][26] := _EnvDSNumericSaveDigits else _dat[4][26] := Digits end if; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [z, a3(z)], (4) = []}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol := 1; _ndsol := _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error end try end proc

> 

p1 := plots:-odeplot(dsol, gridlines=false);

 

 

Maple's symbolic solution:

> 

dsolve({de, ic}) assuming positive;
eval(eval(a3(z), %), {params});
p2 := plot(%, z=0..1, gridlines=false);

a3(z) = -JacobiSN(z*K*a10, a20/a10)*a20

-2.34*JacobiSN(5.2029*z, 1.902439024)

 

 

 

 

Let's plot the numeric and symbolic solutions together.

The numeric and symbolic solutions agree up to z=0.2.  Beyond that

the solution of the ODE is non-unique -- the numeric solution remains

constant, while the symbolic solution oscillates.  Both are true solutios of the ODE.

 

EDIT: Oops!  I take that back.  Both solutions are correct up to z=0.2.  Beyond that the numeric

solution is correct but the symbolic solution is not, since according to the ODE,

the derivative of the solution should be less than or equal to zero!

> 

plots:-display([p1,p2], color=[red,blue]);

 

 

 

 

Mathematica's solution(?) does not agree with the numeric solution.

> 

MM := a3(z) = -JacobiSN(z*a10*K, a20^2/a10^2)*a20;
eval(eval(a3(z), %), {params});
plot(%, z=0..1, gridlines=false):
plots:-display([p1, %], color=[red,blue]);

a3(z) = -JacobiSN(z*K*a10, a20^2/a10^2)*a20

-2.34*JacobiSN(5.2029*z, 3.619274242)

 

 
 

Download mw.mw

 

Use implictplot...

Rouben Rostamian 7711
April 30 2020
0 1

This worksheet shows how to plot those curves.  These do not agree with the picture that you have shown.  Perhaps the picture corresponds to a different choice of parameters.

> 

restart;
> 

with(plots):
> 

eq1 := M^2*(1-sqrt(1-2*x/M**2))+2*f/(3*q-1)*(1-(1+(q-1)*x)^((3*q-1)/(2*q-2)))+2*(1-f)/(b*(3*q-1))*(1-(1+b*(q-1)*x)^((3*q-1)/(2*q-2)))=0;

M^2*(1-(1-2*x/M^2)^(1/2))+2*f*(1-(1+(q-1)*x)^((3*q-1)/(2*q-2)))/(3*q-1)+2*(1-f)*(1-(1+b*(q-1)*x)^((3*q-1)/(2*q-2)))/(b*(3*q-1)) = 0

> 

eq2 := 1/sqrt(1-2*x/M^2)-f*(1+(q-1)*x)^((q+1)/(2*q-2))/(q-1)-(1-f)*(1+b*(q-1)*x)^((q+1)/(2*q-2))/(b*(q-1))=0;

1/(1-2*x/M^2)^(1/2)-f*(1+(q-1)*x)^((q+1)/(2*q-2))/(q-1)-(1-f)*(1+b*(q-1)*x)^((q+1)/(2*q-2))/(b*(q-1)) = 0

> 

sys := eval([eq1,eq2], [b=0.2,q=0.3,M=1.73583]);

[3.013105789-3.013105789*(1-.6637669368*x)^(1/2)-0.2e2*f*(1-(1-.7*x)^0.7142857143e-1)-100.0000000*(1-f)*(1-(1-.14*x)^0.7142857143e-1) = 0, 1/(1-.6637669368*x)^(1/2)+1.428571429*f/(1-.7*x)^.9285714286+7.142857143*(1-f)/(1-.14*x)^.9285714286 = 0]

> 

fsolve(sys);
p1 := pointplot(eval([x,f], %), symbol=circle, symbolsize=25, color=black):

{f = 1.424999999, x = -0.3290821902e-8}

> 

p2 := implicitplot(sys, x=-4..4, f=-4..4, color=[red,blue],
  axes=boxed, numpoints=10000, legend=['eq1','eq2']):
> 

display([p1,p2]);

 

The circle marks the solution of the system {eq1, eq2}.

 

Download mw.mw

 

Try 2D to 2D...

Rouben Rostamian 7711
April 27 2020
0 0

I think the proper way to doing this is to work out the 2D to 2D transformation that deforms your equirectangular projection (I will refer to is as "the picture") onto the final result, without going through the intermediate 3D representation.

The following worksheet can serve as a starting point.  Let 0 < s < 2*Pi and -Pi/2 < t < Pi/.2 be the horizontal and vertical coordinates of the picture.  The mapping
(s, t) -> [cos(t)*cos(s), cos(t)*sin(s), sin(t)]
maps the picture into the sphere in 3D.  This is shown in plot p1.

The mapping
(s, t) -> [sqrt(2*(1+sin(t)))*cos(s), sqrt(2*(1+sin(t)))*sin(s)]
(which is 2D to 2D) maps the picture onto the desired Lambert projection.  This is shown in plot p2. (Well, actually in plot p2 I have embedded the result into 3D in order to later form a composite plot, but the embedded plot is inherently 2D.)

I think it is possible, although I haven't done this myself since it is late at night, to go from the 2D picture to the 3D sphere, rotate the sphere as needed, and then apply the Lambert transform, and figure out the overall map from 2D to 2D.

> 

restart;
> 

p1 := plot3d([cos(t)*cos(s), cos(t)*sin(s), sin(t)],
    s=0..2*Pi, t=-Pi/2..Pi/2,
    image="earth.jpg",
    labels=[x,y,z], size=[300,300]):
> 

p2 := plot3d([sqrt(2*(1+sin(t)))*cos(s), sqrt(2*(1+sin(t)))*sin(s), -2],
    s=0..2*Pi, t=-Pi/2..Pi/2,
    image="earth.jpg", size=[300,300]):
> 

plots:-display([p1,p2],
    scaling=constrained, orientation=[145,72,0],
    axes=none);

 

Download lambert-projection.mw

 

Incorrect solution...

Rouben Rostamian 7711
April 18 2020
1 3

You say "after solving..." but you haven't asked Maple to solve the equation.  The solution u:=... that you have entered manually is incorrect.  If you insist on entering your own solution, then apply the quadratic formula with a little bit more care.  Alternatively, ask Maple to produce the solution for you.

 

See if this is useful...

Rouben Rostamian 7711
April 17 2020
3 2

See if this is useful.

 Download mw2.mw 

External tools...

Rouben Rostamian 7711
April 15 2020
1 0

What you are asking is a job for your operating system, not Maple.

Install fswatch which is available for Linux, Mac, and Windows. Run it in a terminal and watch its notifications.  Then do what needs to be done in Maple.

 

 

First 16 17 18 19 20 21 22 Last Page 18 of 53