pagan

5127 Reputation

23 Badges

16 years, 194 days

 

 

"A map that tried to pin down a sheep trail was just credible,

 but it was an optimistic map that tried to fix a the path made by the wind,

 or a path made across the grass by the shadow of flying birds."

                                                                 - _A Walk through H_, Peter Greenaway

 

MaplePrimes Activity


These are answers submitted by pagan

Have you tried raising the iteration limit, by passing the optional argument iterationlimit=N to the `QPSolve` command, for some positive integer N?

If you apply the command  infolevel[Optimization]:=2 before calling `QPSolve` then you should be able to see the value for `iterationlimit` that it is using by default.

If that doesn't help, then could you upload a worksheet containing the minimal set up to reproduce the problem? For example, it could contain the QPSolve command being uses, and the two 4x4 Matrices and anything else needed to run it.

Since the `sum` command does not have special evaluation rules the argument MatrixPower(x,k) is being computed for nonnumeric symbol `k` before `sum` gets its hands on it.

That is Maple's usual mode of evaluating arguments to a procedures calls, up front.

This is a particular instance of a common case of misuse of the `sum` command. More generally, it is a frequently asked question.

Two popular workarounds are to either delay evaluation of the argument, using single-right uneval quotes, or to instead use the `add` command which does have special evalution rules.

add(MatrixPower(x,k), k = 0..3);  # likely better here than using uneval quotes

sum('MatrixPower'(x,k), k = 0..3);

Both of the above are likely better in general than using the `fnormal` command to remove small imaginary components after the computation (since their size, and thus the cut-off, will depend upon the scale of the entries of Matrix x).

It's possible that you were really trying to compute and use the generic result of MatrixPower(x,k), to try and speed up the repeated powering of Matrix x while computing the summation. That's not such a bad idea -- in itself -- for a float Matrix x. But in that scenario it might be better to diagonalize Matrix x in a more usual fashion (using the Eigenvectors command, say).

 

The name :-n is the global n. Using that protects against picking up the local n at higher scope.

But :-n might be assigned. So safest is to use ':-n' which should also protect against global :-n being assigned.

In other words, inside procedure g1, call g2 like so

  g2(':-n'=n);

Put a colon after the `od`.

@rembremading If you start Maple using one of the interfaces other than OpenMaple (commandline, Standard GUI, Classic) then you should be able to examine LD_LIBRARY_PATH and see what Maple itself thinks is necessary (or best, from what is available shipped).

For example (modify as you need it) too pick off the first directory on Linux:

   StringTools:-Split(getenv(LD_LIBRARY_PATH),":")[1];

You can use whatever Vector (RHS of A.X=B) that you want:

v := Vector( [54*a+3*b-c+d, 32*a-c+d, 96*b-69*c+85*d, 6*a+9*b+3*c+9*d] ):

RHS := Vector( [1.0, 2.0, 7.0, -1.0] ):

Equate( v, RHS );

        [54 a + 3 b - c + d = 1.0, 32 a - c + d = 2.0, 96 b - 69 c + 85 d = 7.0, 

          6 a + 9 b + 3 c + 9 d = -1.0]

{seq(eq, eq in Equate( v, RHS ) )};

        {32 a - c + d = 2.0, 96 b - 69 c + 85 d = 7.0, 6 a + 9 b + 3 c + 9 d = -1.0, 

          54 a + 3 b - c + d = 1.0}

solve( {seq(eq, eq in Equate( v, RHS ) )} );

       {a = 0.03704954955, b = -0.6050300300, c = -0.2588963964, d = 0.5555180180}

A, B := LinearAlgebra:-GenerateMatrix( Equate(v,RHS), [a,b,c,d] );

                              [54   3   -1   1]  [ 1.0]
                              [               ]  [    ]
                              [32   0   -1   1]  [ 2.0]
                      A, B := [               ], [    ]
                              [ 0  96  -69  85]  [ 7.0]
                              [               ]  [    ]
                              [ 6   9    3   9]  [-1.0]

Equate( Vector([a,b,c,d]), LinearAlgebra:-LinearSolve(A, B) );

       [a = 0.0370495495495495, b = -0.605030030030030, c = -0.258896396396396, 

         d = 0.555518018018018]

RHS := Vector( [0.0, 0.0, 0.0, 0.0] ):

Equate( v, RHS );

        [54 a + 3 b - c + d = 0., 32 a - c + d = 0., 96 b - 69 c + 85 d = 0., 

          6 a + 9 b + 3 c + 9 d = 0.]

{seq(eq, eq in % )};

       {32 a - c + d = 0., 96 b - 69 c + 85 d = 0., 6 a + 9 b + 3 c + 9 d = 0., 

         54 a + 3 b - c + d = 0.}

solve( % );

                           {a = 0., b = 0., c = 0., d = 0.}

A, B := LinearAlgebra:-GenerateMatrix( Equate(v,RHS), [a,b,c,d] );

                               [54   3   -1   1]  [0]
                               [               ]  [ ]
                               [32   0   -1   1]  [0]
                       A, B := [               ], [ ]
                               [ 0  96  -69  85]  [0]
                               [               ]  [ ]
                               [ 6   9    3   9]  [0]

Equate( Vector([a,b,c,d]), LinearAlgebra:-LinearSolve(A, B) );

                             [a = 0, b = 0, c = 0, d = 0]

Your approach is quite the tangle of escaped locals and assumptions, I'm afraid.

Your subP3 procedure assigns various procedures to global names, and these procs reference a local C1 upon which your subP3 first places assumptions. When you call those procs after running subP3 then your results (quite naturally) still contain references the (now, escaped) local C1.

Ideally, it could all be rewritten to do things quite differently. But there are too many tangles to do that quickly. Other issues involve use of `int` rather than `Int`, the looped list concatenation, etc. But the central issues revolve around the references to C1 in the bodies of the procedures which get assigned globally (for subsequent use).

One quick workaround would be to instead declare C1 as global to subP3, then query the results from the later calls for local names, and substitute into the result using a mapping of global=local equations. (This might entail conversion of the assumed local to a string, then chopping off the tilde, etc...)

An alternative, also having C1 declared as global within subP3, is to place assumptions on some other name CC1 instead of C1. The assumpton appears to be needed only to get the double integral in the argument to fsolve to compute successfully. Attached is something in this manner, which seems to succeed. Other members might see an easier way.

for_MP.mw

You could interpolate the data, to find the intersection(s). That means that the fitted curves will pass directly through the data points.

One alternative approach is to fit the data through approximation, but without necessarily passing directly through each data point. You could see the Statistics:-Fit command, as one way to try doing that.

 

restart:

with(plots):

ax := [13300, 17400, 12500, 20300, 12100, 22100, 11800, 23100,
       11600, 23800, 11400, 24300, 11300, 24700, 11200, 25100]:
ay := [.2, .2, .3, .3, .4, .4, .5, .5,
       .6, .6, .7, .7, .8, .8, .9, .9]:
cx := [17700, 17900, 18200, 18500, 18900,
       19300, 19900, 20700, 22700, 29600]:
cy := [.1, .2, .3, .4, .5, .6, .7, .8, .9, .9]:

AXY:=Matrix(sort([seq([ax[i],ay[i]],i=1..nops(ax))],(s,t)->s[1] CXY:=Matrix(sort([seq([cx[i],cy[i]],i=1..nops(cx))],(s,t)->s[1]
display({plot(AXY, style = line, color = blue),
         plot(CXY, style = line, color = red)});

# Use method=linear, if you prefer.
#methoption:=method=linear:
methoption:=method=spline:

axycurve:=x->CurveFitting:-ArrayInterpolation(AXY,x,methoption,extrapolate):
cxycurve:=x->CurveFitting:-ArrayInterpolation(CXY,x,methoption,extrapolate):

display({plot(axycurve, min(ax)..max(ax),color=blue),
         plot(cxycurve, min(cx)..max(cx),color=red)});

methoption:=method=spline:

fsolve(axycurve-cxycurve,min(ax)..max(ax));

17944.44743

fsolve(axycurve-cxycurve,min(cx)..max(cx)); # possible extrapolation

25199.98161

methoption:=method=linear:

fsolve(axycurve-cxycurve,min(ax)..max(ax));

17957.69231

fsolve(axycurve-cxycurve,min(cx)..max(cx)); # possible extrapolation

HFloat(25100.0)

 

 

Download intersect.mw

If I understand your problem, then you shouldn't have to call `dsolve` repeatedly. Instead, you could make use of the ability of dsolve,numeric to allow parameters for IVPs. (It doesn't allow for parameters for BVPs, and if it did then maybe this problem could be done entirely by dsolve...!?)

You might even have a choice of rootfinders to use here, eg. fsolve, one from DirectSearch(?), or NextZero. You could probably use an optimization routine as well, instead of a dedicated rootfinder for the outer layer of solving (maybe better suited, as the `objective` used below does not cross the x-axis but merely gets close to it).

 

restart:

stlsl := {-(40.00000000*(0.4469633000e-3*(diff(y(t), t))^2
           +0.4469633000e-3*(diff(x(t), t))^2))*(diff(x(t), t))
           /sqrt((diff(y(t), t))^2+(diff(x(t), t))^2) = diff(x(t), t, t),
          (40.00000000*(0.4469633000e-3*(diff(y(t), t))^2
           +0.4469633000e-3*(diff(x(t), t))^2))*(diff(y(t), t))
           /sqrt((diff(y(t), t))^2+(diff(x(t), t))^2)-9.810000000
           = diff(y(t), t, t), x(0) = 0, y(0) = .80,
          (D(x))(0) = 8.5*cos(0.1744444444e-1*angle),
          (D(y))(0) = 8.5*sin(0.1744444444e-1*angle)}:

gensol:=dsolve(stlsl,numeric,parameters=[angle],output=listprocedure):
X:=eval(x(t),gensol):
Y:=eval(y(t),gensol):

err:=(X(t)-5.44)^2+(Y(t)-1.60)^2:

objective:=proc(ang)
  local res;
  if not type(ang,numeric) then return 'procname'(args); end if;
  X(parameters=[ang]);
  res:=Optimization:-Minimize(err, t=0..2)[1];
  if type(res,numeric) then
     res; #if`(res   else Float(infinity); end if;
end proc:

ans1:=RootFinding:-NextZero(objective,0.0);

62.24001424

X(parameters=[ans1]);
Optimization:-Minimize(err, t=0..2); # if we want to known at what t it attains
plot(err,t=0..2);

[angle = 62.24001424]

[HFloat(3.4608140151503127e-18), [t = HFloat(1.4774685603002573)]]

ans2:=fsolve(objective(A),A,avoid={A=ans1});

35.08802379

X(parameters=[ans2]);
Optimization:-Minimize(err, t=0..2); # if we want to known at what t it attains
plot(err,t=0..2);

[angle = 35.08802379]

[HFloat(1.2438283666257583e-16), [t = HFloat(0.8248851317654872)]]

 

 

Download dsolveroot.mw

 

If my understanding is ok, but the sheet is not clear enough, there are more ways to get at it:

ans:=Optimization:-Minimize( objective(A), method=nonlinearsimplex, initialpoint=[A=70] );
X(parameters=[rhs(ans[2][1])]);
Optimization:-Minimize(err, t=0..2); # if we want to known at what t it attains
ans:=Optimization:-Minimize( objective(A), method=nonlinearsimplex, initialpoint=[A=10] );
X(parameters=[rhs(ans[2][1])]);
Optimization:-Minimize(err, t=0..2); # if we want to known at what t it attains
#plot(objective, 0..90); # only two answers, but that is intuitively clear from the physics
X(parameters=[ans1]):
P1:=plot([X-5.44,Y-1.66],0..2,color=red):
X(parameters=[ans2]):
P2:=plot([X-5.44,Y-1.66],0..2,color=green):
plots:-display(P1,P2); # want both curves of a given color intercepting x-axis at same t

I'm not at all sure if this is what you want. But if not then it might still give you some usage ideas.

ode.mw

One of the problems in your original could be old fashioned premature evaluation, because this part of the arguments to the second dsolve(...) produces an immediate error (ie. simply by being passed). In essence, you had not set up your y[1] procedure so that it would return unevaluated if called with an argument that was not numeric.

y[1](t-1);
Error, (in y[1]) improper op or subscript selector

What is the purpose of the inner for-loop with index variable `i`? How is this inner loop supposed to affect the computation? Or is the presence of the inner loop a mistake?

You're already using `i` as the summation index. The clash is the cause and meaning of the error message you see.

You should consider using the `add` command, not the `sum` command, for adding up a finite number of terms like that.

if b>a and b>c then a; end if;

Is this what you wanted to get? (You have your own values in diagonal Matrix z, of course.)

restart:
with(LinearAlgebra):

z:=DiagonalMatrix([a,b,c,d,e,f]);

                                [a  0  0  0  0  0]
                                [                ]
                                [0  b  0  0  0  0]
                                [                ]
                                [0  0  c  0  0  0]
                           z := [                ]
                                [0  0  0  d  0  0]
                                [                ]
                                [0  0  0  0  e  0]
                                [                ]
                                [0  0  0  0  0  f]

Z:=Matrix(6):

node:=3:
for i from 1 to (node*2-1) do
for j from (i+1) to (node*2) do
Z[i,j]:=(z[i,i]-z[j,j])/2;
end do:
end do:

Z;

         [   1     1    1     1    1     1    1     1    1     1  ]
         [0  - a - - b  - a - - c  - a - - d  - a - - e  - a - - f]
         [   2     2    2     2    2     2    2     2    2     2  ]
         [                                                        ]
         [              1     1    1     1    1     1    1     1  ]
         [0      0      - b - - c  - b - - d  - b - - e  - b - - f]
         [              2     2    2     2    2     2    2     2  ]
         [                                                        ]
         [                         1     1    1     1    1     1  ]
         [0      0          0      - c - - d  - c - - e  - c - - f]
         [                         2     2    2     2    2     2  ]
         [                                                        ]
         [                                    1     1    1     1  ]
         [0      0          0          0      - d - - e  - d - - f]
         [                                    2     2    2     2  ]
         [                                                        ]
         [                                               1     1  ]
         [0      0          0          0          0      - e - - f]
         [                                               2     2  ]
         [                                                        ]
         [0      0          0          0          0          0    ]

For example,

f1:=sin(x): f2:=cos(x): f3:=x^2: f4:=sqrt(x):

plots:-display(
  plot([f1,f2], x=-0.5..0, color=[red,green]),
  plot([f3,f4], x=0..0.5, color=[blue,gold])
               );

Perhaps it is this one, which currently shows as 36443 views ("14698 unique viewers" which I think means distinct IPs).

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