dharr

Dr. David Harrington

8477 Reputation

22 Badges

21 years, 33 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
I am a retired professor of chemistry at the University of Victoria, BC, Canada. My research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

MaplePrimes Activity


These are replies submitted by dharr

@Math-dashti Perhaps this is what you want. It finds the equilibrium points and the Jacobian.

bifur.mw

This version returns the quadratic form in a sum-of-squares form:

restart;

SumSquares takes a quadratic form with real coefficients and returns an explicitly positive sum of squares form if possible, i.e., if the matrix is positive semidefinite.

Set infolevel[SumSquares]:=2 to output the Matrix and its eigenvalues.

SumSquares:=proc(p::polynom)  # polynomial with names as variables and constants as coefficients
  local i,j,term,p1,terms,nvars,vars,tvars,A,a,tf,L,v2,P;
  uses LinearAlgebra;
  p1:=expand(p);
  vars:=[indets(p1,assignable(name))[]];
  nvars:=nops(vars);
  A:=Matrix(nvars,nvars,shape=symmetric);
  if type(p1,`+`) then terms:=[op(p1)] else terms:=[p1] end if; # make list of terms
  for term in terms do
    if degree(term)<>2 then error "not a quadratic form" end if;
    tvars:=indets(term,assignable(name));  # one indet for a*x^2; two for a*x*y
    if nops(tvars)=1 then
      member(tvars[1],vars,'i');
      a:=eval(term,tvars[1]=1);
      if is(Im(a)=0) then
        A[i,i]:=a;
      else
        error "coefficient %1 may have imaginary part",a
      end if;
    else
      member(tvars[1],vars,'i');
      member(tvars[2],vars,'j');
      a:=eval(term,{tvars[1]=1,tvars[2]=1});
      if is(Im(a)=0) then
        A[i,j]:=a/2;
      else
        error "coefficient %1 may have imaginary part",a
      end if;    
    end if;    
  end do;
  if ormap(type,[entries(A,'nolist')],float) then
    WARNING("result may be unreliable for floating point coefficients");
    if Rank(A)<nvars then WARNING("floating point singular matrix increases probability of unreliable result") end if;
  end if;
  userinfo(2,'procname',"Matrix and Eigenvalues:",print(A),Eigenvalues(A,output='list')[]);
  tf:=IsDefinite(A, query = 'positive_semidefinite'); # eigenvalues non-negative
print(tf);
  if not type(tf,truefalse) then error "Matrix is only positive semidefinite if %1",tf end if;
  if not tf then error "Matrix is not positive semidefinite" end if;
  L:=LUDecomposition(A, method='Cholesky');
  v2:=L^+.Vector(vars);
  P:=v2^+.v2;
  if type(P,`+`) then map(simplify,P) else P end if;
end proc:

P:=2*x^2-2*x*y-2*x*z+2*y^2-2*y*z+2*z^2;

2*x^2-2*x*y-2*x*z+2*y^2-2*y*z+2*z^2

SumSquares(P);

true

Warning, Matrix is not positive-definite

(1/2)*(2*x-y-z)^2+(3/2)*(y-z)^2

NULL

Download SumSquares.mw

@Alfred_F The routine I gave here writes a quadratic form in the sum of squares form. I'm not aware of a specific Maple command that does that. But your polynomial is not a quadratic form, so that wouldn't help here.

If by "same line" you mean "in the same execution group (same > prompt) but on the next line" then try shift-enter. If you really mean same line, then don't press enter, just keep typing, separating your commands with semicolons.

Otherwise, please provide a more detailed description of the problem and what you want.

@Suryakanth Did you make those changes?

@salim-barzani I think you want to get rid of z:=0; and later put eval(eqt1, z = 0).

@salim-barzani Sorry, I thought it was wrong to just put z=0 in eq17 by any means, assignment or not. I don't understand what you want to do and have no further comment.

@Suryakanth This time OdeSys and Cond are sets, so to combine them in the dsolve call you need dsolve(OdeSys union Cond,...). The D@@2 need parentheses: (D@@2).

Then you call ContoursWithLabels with argument y=yL..yR but the error message says that it received 0=-0.8..1,5, so you must have set y=0 somewhere.

Look at eq17. A few lines later you assign z:=0. So now eq17 has functions with (x,y,0,t) which presumably don't make sense. This carries through and leads to the D[1](f)(x, y, 0, t) - if you still had the z it would display as an indexed derivative.

This sort of assigning is bad practice in Maple for this reason. Better to use eval(..,z=0) to set z zero only in the things you want and not disrupt all earlier expressions. I'm guessing that you only want z=0 for plotting purposes, so perhaps it should be later? But I'm not sure what is intended here.

@Suryakanth I don't know anything about streamlines. In particular you seem to need psi(x,y) when the odes only gave psi(y).

If you are asking about how to make the contourplot work for the psi_fun you give then note that eval(psi(y), Ans[1]) is a procedure so you need to invoke it at y with eval(psi(y), Ans[1])(y). Also y ranged from -1..1 in that function, not 0..400.

So the following makes the plot work.

pulatile_flow_error.mw

@acer I like it!

@salim-barzani You changed the ode so now there is no w. For the other one {A[0],A[1],A[2],B[0],B[1],B[2],w} takes a very long time and I didn't wait. With 7 equations you have several choices for the 7 variables, so you could try some other options. I don't know how to get a solution.

For the long wave length limit I looked at this before and gave part of the solution, but didn't understand how to get the expression with the theta's. If you have some instructions on how to do that I can take a look.

Do you only have coordinate data or do you also know the identity of the atoms? Are you trying to find the identities from the bond lengths? What class of molecules are you interested in? The variety of bond lengths and hence the tolerance question depends significantly on the class of molecules - what does "heavy-atom molecule" mean?

@Alfred_F Add frames = 20. To run the animation click on the plot and then on the play icon on the toolbar.

@sand15 Yes, as I said the wolf only captures the goat in an asymptotic sense. From a more practical point of view the finite-sized wolf catches the finite-size goat after the separation reaches a small number, which I took to be 0.05. Or the 0.05 can just be a number after which you are fairly certain (without an analytical solution) that the wolf will catch the goat at infinite time. The real question about good/bad starting locations has an uninteresting answer that says it doesn't matter - the wolf either never or always catches the goat, depending on your point of view.

I originally programmed the wolf chasing the pre-programmed goat moving in a circle, just to see I had the right physics, thinking I would need to make the goat more clever. So I then made the odes for the goat, with motion on the unit circle, which leads to just the same thing. Once @Alfred_F confirmed the two speeds were each equal and constant, there was not anything else to do.

If one made the two speed maxima, and allowed some acceleration, the wolf and goat to move in arbitrary directions, and allowed the goat to move within the unit circle, the wolf will still catch the goat. Assuming the wolf doesn't start at (1,0) but say from (1+d,0), then the goat and wolf just chase each other in the -x direction until the wolf catches the goat at (-1,0). If the goat deviates, then I think the wolf can anticipate and move to force the goat back toward the x-axis, and catch it after a few oscillations.

1 2 3 4 5 6 7 Last Page 1 of 89