dharr

Dr. David Harrington

5973 Reputation

21 Badges

19 years, 300 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 professor of chemistry at the University of Victoria, BC, Canada, where 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

restart;

IsQformNonNeg takes a quadratic form with real coefficients and returns true if it is non-negative for any real values of the variables.

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

IsQformNonNeg:=proc(p::polynom)  # polynomial with names as variables and constants as coefficients
  local i,j,term,p1,terms,nvars,vars,tvars,A,a,tf;
  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
  if not type(tf,truefalse) then error "only true if %1",tf end if;
  tf
end proc:

with(LinearAlgebra):

Two cases that are sums of squares and so explicitly non-negative

p1:=expand((x - y/2)^2 + 3/4*y^2);
p2:=expand((x - y)^2 + (y - z)^2 + (z - x)^2);

x^2-x*y+y^2

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

IsQformNonNeg(p1);
IsQformNonNeg(p2);

true

true

is(p1>=0) assuming real;
is(p2>=0) assuming real;

FAIL

FAIL

But if the coefficients are floats, then the (semi)positive definiteness cannot be reliably determined - here the Matrix is singular and the zero eigenvalue has a small negative floating value.

infolevel[IsQformNonNeg]:=2:
IsQformNonNeg(evalf(p2));
infolevel[IsQformNonNeg]:=0:

Warning, result may be unreliable for floating point coefficients

Warning, floating point singular matrix increases probability of unreliable result

IsQformNonNeg: Matrix and Eigenvalues:

Matrix([[2., -1.000000000, -1.000000000], [-1.000000000, 2., -1.000000000], [-1.000000000, -1.000000000, 2.]])

HFloat(-1.1102230246251565e-16) HFloat(3.0) HFloat(3.0)

false

Error testing

IsQformNonNeg(2);

Error, (in IsQformNonNeg) not a quadratic form

IsQformNonNeg(RootOf(z^2+1)*x^2);

Error, (in IsQformNonNeg) coefficient RootOf(_Z^2+1) may have imaginary part

IsQformNonNeg(Pi*x^2);

true

IsQformNonNeg(arcsin(2)*x^2);

Error, (in IsQformNonNeg) coefficient arcsin(2) may have imaginary part

evalf(arcsin(2));

1.570796327-1.316957897*I

IsQformNonNeg(RootOf(z^2-1)*x^2);

Error, (in IsQformNonNeg) only true if 0 <= RootOf(_Z^2-1)

 

 

Download QformProcedure.mw

@Maple_lover1 From a statistics point of view, the maximum likelihood for normally-distrbuted errors leads to least squares minimization. From a calculation point of view, minimizing least squares is easy because the equations you solve (the normal equations) are linear. I seem to recall that a truncated Fourier series has a minimum least squares error. But without these rationales, and perhaps for your application, you can choose whatever you like.

@acer Agreed, but I have no idea what the OP really wants to do. The point of my post was that you could do something without learning MathML, and that leads to the above procedure. As usual, one can get a better result with some more work, but as an educator, I expect the OP to do some optimization, based on their specific requiremenents. My view of MathML is as for postscript, it is something that you could program, but really it is something for a computer to produce. 

@Ali Hassani 

subsup:=(base,sub,sup)->cat(`#msubsup(mi("`,base,`"),mi("`,sub,`"),mi("`,sup,`"))`);
subsup(x,5,7);

@Zeineb LinearSolve appears to be working now, but your new algorithm is diverging. I'm not able to see anything obviously wrong.

@Maple_lover1 You are setting all the coefficients zero to make the whole thing zero. So if the whole thing is zero, the numerator zero will be zero, and you can just work with that (provided the denominator doesn't go to zero).

As for the P values set to 1, I have no idea - I just noticed that that worked. Your problem doesn't say what to do about the sines, but somehow they disappear.

The d_1 is in the third equation in your post - perhaps it is a typo? or something else about the problem?

In terms of solve, there are 3 equations, so I always choose which 3 variables to solve for. In your set 1 answers, you had a[1]=a[1], meaning a[1] can be anything and the others can be in terms of a[1]. So for set 1, I left out a[1].

In your set 2 answers, there is a[2]=a[2], so I left out a[2], but it didn't work. 

In general solving 3 equations for 4 unknowns won't give anything, but if it does (as here) you don't get to choose what depends on what. Here there is some dependence between the variables that I didn't think about, but which has something to do with why set 2 didn't come out as I expected.

@robertocooper The initial or boundary conditions are determined by the physical problem, and the solution depends on them, so you need to know them for your case. Here I plot for U(0)=0, D(U)(0)=1.

PDE2.mw

 

@robertocooper The solution still has two unknown constants _C1 and _C2, whose values can only be found if you specify two boundary/initial conditons, e.g. U(0)=0, D(U)(0)=1. Then plots can be made, but since the PDE solution is complex what exactly do you want to plot?

This is hard to diagnose unless you upload a .mw file that illustrates the two results. Use the fat green up-arrow, choose the file and upload with "insert link".

@jat741 unapply turns a result you already have into a function, which you can then use for plots etc. In general you should not cut and paste from output to input or retype results.

PDEPlot.mw

@robertocooper I missed the point that the eqn for v was to be derived from the imaginary part, and have added an edited version transform4.mw to my answer above. In your second example, you divided by exp(I*mu), but that left you with a mix of old and new variables. I tidied it up a bit. (I usually leave off the =0 in eqns since dsolve, solve etc understand this as default, and then I do not have to worry whether various simplifications apply automatically to both sides of an eqn or not.)

For the subs, I'm not sure if you are asking a mathematical or Maple question. For the mathematical question, abs(z) is the magnitude of the complex quantity z=r*exp(I*theta), i.e., abs(z)=r. From the Maple point of view, using subs does a substitution without doing any simplification, i.e.,  just syntactic substitution. I didn't want any attempt to try to apply diff to abs. More usually I would used eval(expr,a=b) rather than subs(a=b,expr) since it is smarter about math things. But here eval messed with the PDE too much. 

newcode2.mw

Use I for sqrt(-1) in Maple (or use interface(imaginaryunit=i) to change the default). To force evaluation to real and imaginary parts use evalc. I assumed everything was positive, but you can use some other assumptions if that is not right, You can get simpler forms using simplify, but the abs() is problematic here, 

transform.mw

Providing an approximate solution with the approxsoln option often works here - it can be anything with approximately the right shape, but you have to know something about what you expect.

@dharr I realised that the Q value is intended to select the right phase; the Maple help and coolprop website are not very clear on this.

So this is a better way to do this (worksheets not displaying these days):

Thermo.mw
 

 

@Ronan I edited my original answer with an improved code. Your code agreed with the Matrix. This may be useful in my work, so thanks for introducing me to the Boole-Mobius transform.

First 39 40 41 42 43 44 45 Last Page 41 of 62