dharr

Dr. David Harrington

8340 Reputation

22 Badges

21 years, 5 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

@nm If a window has the focus, then holding down the windows key while you hit the right-arrow key puts the window on the right half of the screen; with left arrow to get on left half of the screen

@Alfred_F I immediately thought of adding "assuming m<>n". But that does not remove the terms for which m=n. So assuming adds information. Your idea of a constraint in this example is really a modification of the sum itself, and the solution is to modify the summation into two pieces that omit the offending terms. In that sense @mmcdara and I gave the correct answer. Maple does not have a way of automating this, and sum itself does not have this capability.

In terms of more general commands, in some cases adding an assuming clause will work, where an answer may be obtainable for some conditions on the parameters. For example, there are cases in my worksheet where adding assuming n::posint worked. This is probably as close as it gets to a real answer to your question.

In general, Maple gives generic solutions, which apply for most cases but not all. Several solvers have the option parametric, which breaks down the various cases (Summation, solve etc). For differential equations there is casesplit.

In some cases the conditions/constraints are passed as arguments to the commands, such as solve({eqn, a>0}) or commands in the Optimization package.

Bottom line is it depends on what you are doing.

@salim-barzani If you have exact instructions that can reproduce this, it might be possible to diagnose it.

@Alfred_F The last line of Carl's code does this, but only for a given value of k. So

combine(sin(3*x)*cos(x)^3)

gives

1/8*sin(6*x)+3/8*sin(4*x)+3/8*sin(2*x)

As usual you can go the other way with expand. But few Maple commands will deal with the general case; diff is an exception. It is this issue that makes your example hard. For example

int(sin(3*x)*cos(x)^3/x, x=0..infinity);

immediately gives 7*Pi/16.

@salim-barzani For indefinite integration there is always an arbitrary constant (which Maple doesn't show). If the antiderivative is F(x), then indefinite untegration gives F(x) + c, or just F(x) in Maple. Definite integration from -infinity to x gives F(x) - F(-infinity) and the assumption is that F(-infinity) = 0 since we haven't accumulated any area yet at the start. Maple has Intat to describe integration at a point - you could use Intat(f(a),a=x) instead of indefinite integration.

The main Maplesoft site has said "site under maintenance" with possibly reduced functionality for about the same period (still true today), so looks like the revamp took Mapleprimes offline. The application center was also down for at least some of it.

@nm @Alfred_F Since odetest gives 0 for y(0)=a, it is probably doing something simple like below, though I didn't check its code.

dsolve.mw (won't display right now)

@Alfred_F Maple agrees with you.

expand(PDEtools:-dchange(y(t)=sqrt(u(t)),-ode/t,[u(t)]));
DEtools:-odeadvisor(%);

@acer Thanks for the vote and tips, especially for the limit. I knew there would be another answer for the labelling issue - I tried a few things with typeset and uneval quotes, before moving on.

I've become sold on the advantages of equations that are expressions implicitly equal to zero (I know @nm holds the opposite view). Then there aren't things to cancel, just things to divide by. Many things don't map across the `=` and I find normal and numer etc convenient in this context. It has the disadvantage that it is less readable, but I can live with that.

I solve this problem in Maple in a post here.

Incidently, evalf(limit(...)) finds the wrong value of the limit - seems to find the maximum point. I've submitted an SCR.

restart

S := simplify(n/(2^n)*sum(2^k/k, k = 1 .. n));
L := limit(S, n = infinity);
evalf(L);

-LerchPhi(2, 1, n)*n-I*2^(-n)*Pi*n+1

limit(-LerchPhi(2, 1, n)*n-I*2^(-n)*Pi*n+1, n = infinity)

2.666666667+0.4183075556e-21*I

plot(S, n = 0 .. 200)

NULL

Download limit.mw

@MichalKvasnicka Here is a routine to find the matrix exponential by solving the ode system dX/dt = M.X. I haven't optimized it since it seems you might be wanting to port to another language. If you want it further optimized, I can do that, especially through the parameter mechanism, where you only find the procedure to solve the odes once, and then evaluate it for multiple different parameter sets.
[Edit: now calls dsolve only once.]

Evaluating a matrix exponential by solving odes.

Mapleprimes

restart

Make the display precision 20, so we only see this many digits even if we calculate with more

with(LinearAlgebra); interface(displayprecision = 20)

Conditions

c1 := `&omega;__1`+`&omega;__2` = 1; c2 := f__p+f__d1+f__d2 = 1

omega__1+omega__2 = 1

f__p+f__d1+f__d2 = 1

Choose some exact parameters that satisfy the above

params_exact := {T = 7/10, f__d1 = 2/11, f__d2 = 6/11, f__p = 3/11, `&lambda;__1` = sqrt(2), `&lambda;__2` = sqrt(3), `&omega;__1` = 2/9, `&omega;__2` = 7/9}; eval([c1, c2], params_exact)

{T = 7/10, f__d1 = 2/11, f__d2 = 6/11, f__p = 3/11, lambda__1 = 2^(1/2), lambda__2 = 3^(1/2), omega__1 = 2/9, omega__2 = 7/9}

[1 = 1, 1 = 1]

The matrix we need the exponential of

FT := T*Matrix([[-(f__d1/f__p+1)*`&lambda;__1`, -`&omega;__1`*f__d2*`&lambda;__2`/f__p, `&omega;__1`, 0], [-(`&omega;__2`*f__d1/(`&omega;__1`*f__p)-1)*`&lambda;__1`, -(`&omega;__2`*f__d2/f__p+1)*`&lambda;__2`, `&omega;__2`, 0], [0, 0, 0, 1/T], [0, 0, 0, 0]])

Matrix(%id = 36893490500266617660)

FT_exact := eval(FT, params_exact)

Matrix(%id = 36893490500266591404)

Symbolic exp from Maple

expFT := MatrixExponential(FT)

Substitute exact values in, and evaluate at 100 digit accuracy to give the "real" numerical answer.

expFTexact := eval(expFT, params_exact); expFTnum := evalf[100](expFTexact)

Matrix(%id = 36893490500418890916)

Substitute in numerical values at default Digits setting

Digits; params_num := evalf(params_exact); expFTnum1 := eval(expFT, params_num)

10

params_num := {T = .7000000000, f__d1 = .1818181818, f__d2 = .5454545455, f__p = .2727272727, `&lambda;__1` = 1.414213562, `&lambda;__2` = 1.732050808, `&omega;__1` = .2222222222, `&omega;__2` = .7777777778}

Matrix(%id = 36893490500412336476)

Give Maple's MatrixExponential a numerical matrix - it may use a different algorithm

FTnum := eval(FT, params_num); MatrixExponential(FTnum)

Matrix(%id = 36893490500412326844)

Routine to find it by solving a system of ODEs (works at hardware precision)

expM := proc(M::Matrix(square, numeric))
  local n, X, X0, eM, tvar, des, ics, i, xvars, ans, x;
  n := upperbound(M)[1];
  eM := Matrix(n, n);
  xvars := [seq(x[i](tvar), i = 1..n)];
  X := Vector(n, xvars);
  X0 := Vector(n, 0);
  des := seq(diff(X, tvar) - M.X);
  ics := seq(eval(X, tvar = 0) =~ X0);
  ans := dsolve({des, ics}, numeric, 'stiff' = true);
  for i to n do
    X0[i] := 1;
    ans('initial' = [0, seq(X0)]);
    eM[.., i] := eval(X, ans(1));
    X0[i] := 0;
  end do;
  eM
end proc:

expFTnum2 := expM(FTnum)NULL

Matrix(%id = 36893490500412323108)

Maximum error of any entry

Norm(expFTnum2-expFTnum)

0.225022661212714858e-6

NULL

Download ExpByODEs.mw

@MichalKvasnicka So at the end you want to supply numerical values of the parameters and output a numerical matrix, is that correct?

In general the methods Maple chooses for symbolic solutions are different from those for numerical calculations, and putting numbers in at the end may give numerical errors that would not arise from methods optimized for numerical issues. This is especially the case for matrix exponentials. An instructive read is "Nineteen Dubious Ways to Compute the Exponential of a Matrix"doi: 10.1137/S00361445024180 or here. The summary is unless you have a normal matrix there are no general reliable methods, and special difficulties can arise if your matrix is non-diagonalizable (defective; missing eigenvectors) as is true for your case.

If you continue with your strategy and you use Maple at the end of the day, experiment with increasing Digits to make sure you have convergence.

Conceptually, the ode solution method is probably the easiest to implement (with a stiff method - in Maple use stiff=true).

But since the need for the matrix exponential come about from a way to solve some linear equations, can you solve them directly without needing the matrix exponential at all?

For more help, it would be more useful to know more specifically what you intend to do.

@acer  - thanks for your question that prompted the clarification - I had assumed it was for compact display, such as would be required in a description of a method or proof.

@MichalKvasnicka That's a good idea, but you have many variables, so I'm not sure if this meets your definition of compact.

In answer to your questions, you can just pass the whole Matrix:

optimize(expFT, 'tryhard')

 

@MichalKvasnicka In my experience, it is always better to remove dependent variables early. (Non-dimensionalization is one way to do this.) For example, compare 

CharacteristicPolynomial(F, x);
CharacteristicPolynomial(F2, x);

and you will see that 

has been simplified to 1 in the x^2 coefficient. I chose to leave omega__2 in and not omega__1 because some other expressions were simpler. Of course you may lose some symmetry this way, but I would argue that extra variables are more likely to be a problem - also computationally this speeds things up.

You can try other substitutions, and I would be guided by nondimensionalization in selecting these.

In some commands you can use implicit=true, and get your answer in terms of RootOfs, which can then be given a symbolic label, but this is not available for MatrixExponential.

You can try algsubs, e.g. algsubs(P=above expression, F) but the type of expressions it can handle is limited, e.g., it won't do things with sqrts in.

subs will sometimes work, e.g. subs(f__d2*lambda__2*omega__2 + f__d1*lambda__1 + f__p*lambda__1 + f__p*lambda__2 = P, g) does work where g is the output of CharacteristicPolynomial(F, x); But this depends on how the expression is represented natively, which is not always the way it appears to be, so one needs to have a good knowledge of Maple to do this.

A major problem is picking out a subexpression amongst many. subsindets can help with this, but is more advanced Maple; I'm not sure of your expertise here. Some people like pattern matching, which I don't have experience in.

So the real answer to the question is to understand what Maple is doing in terms of the algorithm. If your matrix were diagonalizable one can formulate it in terms of exps of the eigenvalues, from which the awkward square roots arise. But we are short one eigenvector, so the exp is done through the Jordon form and the primary matrix function, involving derivatives evaluated at the eigenvalues. That seemed too complcated, so I didn't pursue it.

I'd say this sort of thing is in general very difficult, as you supposed.

First 9 10 11 12 13 14 15 Last Page 11 of 87