Robert Israel

6577 Reputation

21 Badges

18 years, 215 days
University of British Columbia
Associate Professor Emeritus
North York, Ontario, Canada

MaplePrimes Activity


These are answers submitted by Robert Israel

This seems to be a bug in Maple 11 Standard, not present in Classic.  Don't bother submitting a bug report, though, because it has been fixed in Maple 12.  The Classic version of loglogplot is still available in Standard, as `plots/loglogplot`.  So, e.g.:

> `plots/loglogplot`(exp(ln(x)^2),x=1..100);

 

In any differential equations model you can run time backwards.  That is, instead of specifying an initial condition and seeing how the system evolves into the future, you specify a final condition and see how the system must have evolved in the past to attain that final state.  If you're examining an unstable equilibrium, for example, it may be practically impossible to arrange an initial condition far away from that equilibrium that will approach the equilibrium very closely in the future.  It's much easier to take a final condition near the equilibrium and run time backwards.

If you change each diff(x[j](t), t) to -diff(x[j](t),t) in an autonomous system of DE's, you're in effect making time run backwards instead of forwards.  An unstable node or spiral (where all eigenvalues have positive real part) then becomes a stable one (where all eigenvalues have negative real part).  On the other hand, a saddle (where some eigenvalues are positive and others are negative) remains a saddle.

Changing two equations out of three in a system may or may not make sense: maybe there's a symmetry that changes the direction of time and also changes the sign of
some of the variables.

The equation for y in terms of x is

eq:= diff(y(x),x) = e*( ( a + b* y(x) + c* x+ d* x* y(x)+ f *y(x)^2 )  / ( a + b*x + c*y(x) + d* x*y(x) + f* x^2 ) );

dsolve(eq);

seems to be difficult for Maple 12: it was well into virtual memory and slowing my computer to a crawl when I aborted it.

I tried the special case with e=2 and all other parameters = 1, and no closed-form solution was found.

DEtools[odeadvisor] classifies the equation as Abel type 2B.  I'm not sure what transformation it has in mind.  In general, these don't have known closed-form solutions.  There are methods that work for some of them, but apparently not this one.

The special case with all parameters = 1 does have closed-form (implicit) solutions:
 

ln(x^2+2*x+2*x*y(x)+2+y(x)^2+2*y(x))-2*arctan(x+y(x)+1)-2*ln(-x+y(x))-_C1 = 0

Presumably, your second order differential equations express those second derivatives in terms of the variables and first derivatives.  So those expressions are what you give to odeplot.  For example, something like this:

> sys := { diff(x(t),t$2) + x(t)*diff(y(t),t$2) - 
      x(t)*diff(y(t),t) = 1,   
         diff(x(t),t$2) - diff(y(t),t$2)+diff(x(t),t)+y(t)=0};
Rsys:= solve(sys,{diff(x(t),t$2),diff(y(t),t$2)});
S:= dsolve(sys union {x(0)=0,D(x)(0)=1,y(0)=0,D(y)(0)=1},numeric);

plots[odeplot](S,[y(t),subs(Rsys,diff(x(t),t$2))],t=0..5);

The "saddle paths" are, I think, the trajectories entering and leaving the equilibrium point.  These go in the direction of the eigenvectors of the linearized system: in this case the system is actually linear so they are just straight lines in the direction of the eigenvectors of the matrix.  For your system:

>  A := <<1|1>,<4|1>>;
  E,V:= LinearAlgebra[Eigenvectors](A);
 

Here E[1] is negative and E[2] is positive, so the incoming trajectories correspond to the first column of V, namely <-1/2,1>, and the outgoing trajectories correspond to the second column of V, namely <1/2, 1>.
 

If you're using the t interval -10..10, you might use initial conditions for the incoming trajectories near the beginning of this interval and for the outgoing trajectories near the end of this interval, e.g.

> incoming := {[x(-8)=-1/2, y(-8) = 1], [x(-8) = 1/2, y(-8) = -1]};
   outgoing := {[x(8) = 1/2, y(8) = 1], [x(8) = -1/2, y(8)=-1]};

And then try your DEplot using incoming union outgoing instead of your initial.

 

 


 

Here are your equations, with I replaced by i as Alec suggests.

> eqs:= {Mt=M+K1*M*i+K2*M*i^2+K3*M*Z,
  It=i+K1*M*i+K2*M*i^2,
  Zt=Z+K3*M*Z,
  650*i+36800*K1*M*i+51000*K2*M*i^2=Aobs};

These four equations can be solved for four unknowns, say K3, M, Z, i:

> solve(eqs, {K3,M,Z,i});

The result is rather complicated.  Since it contains RootOf a quadratic polynomial, there will generally be two solutions, which may or may not be real, depending on the values of the parameters.  The discriminant of the quadratic is

-129100*K2*Aobs*K1+3724220000*K2*K1*It+Aobs^2*K2^2-102000*Aobs*K2^2*It+2601000000*K2^2*It^2+1306822500*K1^2

Oops: the <maple> tag is acting up again... that should be

-129100*K2*Aobs*K1+3724220000*K2*K1*It+Aobs^2*K2^2-102000*Aobs*K2^2*It
+2601000000*K2^2*It^2+1306822500*K1^2

Assuming the parameters are real, there are two real solutions if this is positive, one if it is 0, two if it is negative.

There certainly is a Stop button, but it does not always work: see recent discussions here and here.  It does sometimes happen that a Maple kernel is still running even though the Maple window has closed.  If you're in Windows XP, Ctrl-Alt-Del should bring up a Task Manager, where you can end the mserver.exe process.

As for Matrix in the sidebar, here's a step-by-step example.

1) Place the cursor in a 2D input region.

2) In the Matrix sidebar, set Rows to 3 and Columns to 3, click Insert Matrix.  A 3 x 3 Matrix appears in the input region, with magenta-coloured entries m[1,1] etc., and the top left entry highlighted.

3) Type the top left entry, followed by the Tab key to move to the next entry, and so on until you have all the entries.  Then press Enter.  The Matrix appears as output.

4) If you want to work with this the traditional way, assign it to a variable and use the usual LinearAlgebra functions.  Or if you prefer point-and-click, right-click and choose Standard Operations and then Determinant or Inverse.

 

 

 

1)

> Q:= subs(
solve({de1,de2,de3,de4},{seq(diff(theta[j](t),t$2),j=1..4)}),
     diff(theta[1](t),t$2)):
> plots[odeplot](res, [t, Q], t = 0 .. 2, refine = 1);
      

2)

> plots[odeplot](res, [t, V40], t = 0 .. 2, refine = 1);

 

Actually I think the best command in Optimization for this problem is LSSolve.

> with(Optimization):
   f:= x -> ((A*exp((2*Pi*x)/L))/(1+exp((2*Pi*x)/L)))+((B*exp((2*Pi*x)/L))/(1+exp((2*Pi*x)/L))^2);
  residuals := map(p -> (f(p[1])-p[2]), data):
  R:= LSSolve(residuals);

[.52235658391185e-4, [A = .944204083162349280e-2, B = .162589652654836492, L = 3.54572177530242394]]

Thus the fitted expression is

> ffit := eval(f(x), R[2]);


 

ffit := .944204083162349280e-2*exp(.5640600496*Pi*x)/(1+exp(.5640600496*Pi*x))+.162589652654836492*exp(.5640600496*Pi*x)/(1+exp(.5640600496*Pi*x))^2

 

To plot it and the points:

> with(plots):
   display([pointplot(data), plot(ffit, x=data[1,1]..data[-1,1])]);

Unfortunately it's not a great fit.

Specifically, the problem here is that in Maple / can't be followed directly by -.

Instead of   ... / -0.57... you want ... / (-0.57...).  Similarly, there are cases of

... - -0.5...  which should be ... - (-0.5...) or ... + 0.5... .

This works:

> (((((X1  * (((X3  / X1 ) + (X2  * ((((((X1  * (((X1  / X1 ) 
+ (X4  + X1 )) / (X4  / ((X2  / X3 ) / (-0.570992561856773))))) 
- (X1  - (-0.5457363063736358))) * X2 ) + X1 ) + X3 ) / X2 ))) 
/ (X4  / ((X4  + X2 ) / (-0.570992561856773))))) - (X1  - 
(-0.5457363063736358))) * X2 ) + X1 ) * (((0.2061949348454455 * 
((-0.6832734709973418 + X3 ) / (-0.570992561856773))) - X1 ) / 
(X4  - 0.2061949348454455)));

And to simplify it, you can try

> simplify(%, size);

but it doesn't really simplify much.

 

Maybe you didn't execute the commands.  Here it is in a worksheet.

View 4541_amm2bx.mw on MapleNet or Download 4541_amm2bx.mw
View file details

Suppose your text file consists of lines of the form

c, x, y, z

where c is 1 to 5 (for the colour code) and [x,y,z] is the point to plot, using commas as separators.

> M := ImportMatrix("c:/mypath/myfile.txt", source=csv);
   colours:= [red, green, blue, cyan, black];

The following will work in Standard, but not in Classic.

> plots[pointplot3d](M[1..-1,2..4],
      color=convert(map(t->colours[t], M[1..-1,1]),list));

Alternatively, you might separate out the points of each colour, plot them separately, and combine the plots with display.  This will work in both Standard and Classic.

>  with(plots):
   display([seq](pointplot3d(M[
   select(t -> (M[t,1]=j),[$1..LinearAlgebra[RowDimension](M)])
         , 2..4], 
       colour=colours[j]), j=1..5));

 

Why the TeX?

It might look nice this way, but if you want us to be able to test your result in Maple it's much better to enter Maple commands in ordinary text (1D input).  This way we have to enter the whole thing ourselves by hand rather than copy-and-paste. 

Anyway, the integral Maple ends up doing is

Int(-(exp(cos(t))*cos(sin(t))-1)/(exp(2*cos(t))-2*exp(cos(t))*cos(sin(t))+1)*sin(t)
+exp(cos(t))*sin(sin(t))/(exp(2*cos(t))-2*exp(cos(t))*cos(sin(t))+1)*cos(t),t = 0 .. 2*Pi)

Maple 11 finds an elementary antiderivative, namely

-1/2*exp(t*I)+1/2*ln(exp(2*I*sin(t))+(-exp(2*cos(t))-1)*exp(-cos(t))*exp(sin(t)*I)+1)

Unfortunately, this antiderivative has a jump discontinuity at t=Pi, where
exp(2*I*sin(t))+(-exp(2*cos(t))-1)*exp(-cos(t))*exp(sin(t)*I)+1 crosses the branch cut of ln on the negative real axis.  Failure to detect this discontinuity results in the incorrect value for the integral.

When I tried this with Maple 12, it found a different antiderivative:

-cos(t)+1/2*ln(exp(2*cos(t))-2*exp(cos(t))*cos(sin(t))+1)

which does not have this problem, and results in the correct answer of 0.

Do you mean like this?

> for Q in [P,P1,P2,P3] do
       fprintf("c:/mypath/myfile.txt","%10.10f\n", 
            Matrix(51,3,(i,j)->Q[i][j]))
  end do:
  close("c:/mypath/myfile.txt");
First 99 100 101 102 103 104 105 Last Page 101 of 138