acer

30914 Reputation

29 Badges

19 years, 35 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

Your instincts are right. It can be done using seq(), and a loop which looks like,

S := S, ...

is to be avoided.

But I ought to ask, do you really need to create all 10 operators, which act by evaluating a piecewise?

Or do you just need that f[i](j) returns the value in M[i,j] , or possibly that f[i](y) returns M[i,trunc[y]] ?

There are two parts of the above question. The first is this.  Do you really need 10 operators, f[1]..f[10] (where f was presumably just a Maple table) or would it suffice to have one procedure f which was able be called in that indexed manner?

The second is this. Since the operators (10 individual ones, or 1 indexed proc, no matter) are apparently only going to be used to evaluate a piecewise (in x) at x=y then do you really need to utilize piecewise constructs for that? Maybe these operators don't need to use piecewise constructs at all (since those are, relatively speaking, expensive to evaluate at a specific point).

Some ideas follow, which you might find useful.

M := LinearAlgebra:-RandomMatrix(10,18):

f := proc(y)
#global M; # or by lexical scoping
local i;
  if type(procname,indexed) and type(op(1,procname),posint) then
    i:=op(1,procname);
  end if;
  M[i,trunc(y)];
end proc:

f[2](1);
f[2](1.3);
M[2,1];

Perhaps you really did need what you described, for some other reasons other than just the returned value M[i,j]. You noticed that this invocation below doesn't work, since seq() sees all the arguments flattened out (and hence invalid).

seq( j < x, MM[j,j]*`if`(j=1,0,M[i,j]), j=1..3 );

But maybe you could do something like this,

B := j < x, MM[j,j]*`if`(j=1,0,MM[i,j]);
seq( B, j=1..3 );

acer

Very nice, Axel.

But it seems to me that you have helped Maple out more than you should have had to, by "giving" it the FTOC result.

How about this: instead of,

D(F)(x): int(%,x=0..t) = F(t) - F(0);

you could have,

D(F)(x): Int(%,x=0..t) = int(%,x=0..t,continuous);

so that Maple does more of the work.

 acer

You forgot to define beta[1] as  value*deg .

By the way, I really prefer Units[Standard] to Units[Natural]. The Units[Natural] environment robs the global namespace of too many common symbols which one might want otherwise to use as variables.

> restart:

> with(Units[Standard]):

> alpha[1]:=4.25 * Unit(deg):
> theta[1]:=11.31 * Unit(deg):
> theta[2]:=78.69 * Unit(deg):
> alpha[2]:=-9.25 * Unit(deg):
> a:=1400.00:
> b:=509.90:

> beta[1]:=33.33 * Unit(deg): # I made that up.

> -b*sin(alpha[2]+alpha[1]-beta[1]-theta[2])+b*sin(beta[1]-theta[2]);

                              91.4313519

acer

Type cs\_e instead of cs_e to get it in 2D Math input.

acer

Check this for correctness. Since (3) is missing in the scan, some deduction was necessary. And sometimes I make mistakes.

The request is to make ETsolver accept a general function f of two arguments (which returns dy/dx). That's fine. But why embed into ETsolver the initial condition y[1]=a when that is only true for the particular supplied initial condition y(0)=0 ? And why not also let ETsolver accept a general endpoint b, instead of hard-coding it to only compute an approximation for y(1)?

> ETsolver := proc( f, a, y_a, h, b)
> local N,y,n;
> N := trunc( (b-a)/h );
> y[1]:=f(a,y_a):
> for n from 1 to N do
> y[n+1]:=y[n]+h/2 * ( f(a+(n-1)*h,y[n]) + f(a+n*h,y[n]+h*f(a+(n-1)*h,y[n])) );
> end do;
> return y[N+1];
> end proc:

> ETsolver( (x,y)->x+y, 0, 0, 0.1, 1);
0.7140808465

> ETsolver( (x,y)->x+y, 0, 0, 0.01, 1);
0.7182368623

> ETsolver( (x,y)->x+y, 0, 0, 0.001, 1);
0.7182813769

Like I mentioned before, the key to understanding this algorithm is to draw the exact curve, and the trapezoid with base from x[n] to x[n+1]=x[n]+h , and to see how it relates (as an approximation) to the fundamental theorem of Calculus.

acer

What have you got, so far, with this?

Do you expect complete solutions to your assignments exercises, without making the major contribution yourself? If not, then please just post your instructor's email address so that we can mail in our solutions directly.

It's asking you to write a procedure that implements Euler's method for numerical integration.

That procedure will accept input function f, initial point a, and stepsize h. You could also make it accept final target point b, so as to more generally be able to approximate y(b) instead of y(1).

It refers to equation (3), which is not shown. It's pretty clear, though, that equation (3) is equivalent to (Maple syntax),

y(1) = Int( f(x), x=a..1 );

or

y(b) = Int( f(x), x=a..b );

The procedure could figure out how many steps of length h it will take to get from point a to point 1 (or b). It needs to run a loop, for each such step. Each time through the loop, it will compute y[n+1] using y[n]. It starts off with y[1]:=a .

Here's a hint. To gets from initial point y[1]=a to final point y[N]=b you need to add N increments of stepsize h = (b-a)/N . So you can deduce that N = (b-a)/h . So, you can figure out how many times to go through the loop.

Here's another hint. Try it without a loop, for practice. Try to get it to compute y[2] . From your equation (4) that should equal y[1] + h/2 * (F[1,n]+F[2,n]) .

Try to draw it by hand, as a picture on paper. What do the F[1,n] and F[2,n] mean? Where is the trapezoid? Is the area of the trapezoid almost equal to the area under the curve f between x=a and x=a+h ? Draw it all.

acer

What else can you say about the nonlinear equations? For example, are they (or their subsets that you created) multivariable polynomials?

Are your data sets exact or floating-point approximations, or formulas involving mixed data? Would you consider conversion of float data to exact rational prior to systematic separation of equations?

It sounds like a dimensional effect, as you describe it. For systems of equations in more than one variable, fsolve uses Newton's method. With many variables and hence a much "larger" space, it often becomes more and more unlikely that a starting point will converge. I don't think that fsolve adjusts its number of attempts (distinct initial points) according to the size of the system, and there's no option to specify the number of attempts. So, yes, as the size of the system gets large it might become less and less likely that it will find a n approximate root.

Your idea of separating and reducing the system sounds good, then. But, do you do it "by hand"? Maple's `solve` routine might help you automate it. Or you might be able to use Groebner, depending on the type of equations.

It's possible that someone might be able to do something more with it, if you can upload the equations and some data values to this site.

acer

Please post the example, the actual polynomial. (If you'd like, you can lprint it and put that up here in a reply.)

The fsolve routine is at least sometime capable of doing what you'd like. In Maple 8, 10, and 11 the example below returns  a result with multiple root 3.0 shown twice.

p := expand( (x-3)^2 * (x+4)^3 );
fsolve(p,x=2.5..3.5);

Maybe there's somthing trickier going on with your example. It might be hard for some here to pin down what that might be, if it is specific to Maple 8. Or it could just be something simple, such as all your roots getting put into a maple set at some point (which would "uniquify" them).

Here's a routine that should return the roots of a univariate polynomial which lie in a complex box. The box is specified by lower left corner `a` and upper right corner `b`. In your case, a and b would be the endpoints of a range on the real line. This is similar to how Matlab's `roots` function works, I believe. It should show the multiplicity.

R := (p, a, b) -> select(
    x -> evalb(
        Re(x) <= Re(b) and Im(x) <= Im(b)
        and Re(a) <= Re(x) and Im(a) <= Im(x)),
    convert(LinearAlgebra:-Eigenvalues(
              evalf[trunc(evalhf(Digits))](
                  LinearAlgebra:-CompanionMatrix(p))),
            list));

That routine R won't round or polish the roots, remove small imaginary components, or remove duplicates, etc. That's on purpose, so that it can be you who judges the roots.

y:=randpoly(x,degree=100):
fsolve(y,x=-10..10);
R(y,-10,10);

 You might also be interested in searching the web for Bezout's theorem, roots, and Companion matrix. Here's one link that might interest you.

acer

Please post your example, and maybe it could be improved. acer
First, consider creating an operator which simplifies under your given equation. use_assump := x -> simplify(simplify(x,{a^2+b^2=1})): Now let's look at this Matrix, m := Matrix([[a^2+b^2-1,1,0],[1-a^2-b^2,1,0],[1,0,1]]); with(LinearAlgebra): evals,evecs := Eigenvectors(m); If we now apply the simplification on the eigenvectors Matrix (stored as columns of evecs) then it goes wrong. Presumably this is because it accidentally used as a pivot in the NullSpace computation something that was really zero (eg. an element like 1-a^2-b^2). It fails to simplify the eigenvectors, because they have hidden zeros in the denominators. map(use_assump, evecs); One possibility is to apply the simplification to Matrix m prior to computing the eigenvectors. That works here. But just by luck I think. For another example it too could go wrong as above, I think. But I give it next, anyway. evals,evecs:= Eigenvectors(map(use_assump,m)); map(use_assump, m.evecs - evecs.DiagonalMatrix(evals)); Now, what if the NullSpace computation itself could use the assumption, and so avoid selection of a "hidden" zero as any pivot? One might hope that NullSpace uses Testzero & Normalizer just as LUDecomposition does, and that this benefit also goes for Eigenvectors too. Normalizer:=use_assump: evals,evecs := Eigenvectors(m); map(use_assump, m.evecs - evecs.DiagonalMatrix(evals)); Normalizer := normal: So, maybe that is an Ok method. acer
Jacques reply got quite to the centre of it, I thought. The important part for me is the statement which begins, "It is only various contexts which give an operational semantics to = ". Poor old = , it has too many uses in Maple. For example,
if i = 9 then  ...  end if;

seq( x, x = [1,3,5] );

subs( x = z/2 , cos(x) );

eqn := T = 13/Q - P;

LinearSolve( A, B, method = LU );

G:=table([x = 4]);
As pointed out, the operational semantics of those first three are determined by the context. And even in the fourth example, it is the context which gives a semantic meaning as an inert form to = (if that makes sense to say). I'd guess that there are other interesting examples of use of = , which I've missed. Now for the fun part. What happens if you create a Maple package which exports, or overload, = ? Which of the above examples should or should not be affected by rebinding = ? In which examples above is = in fact so rebound by loading such a package? Are there always workarounds, in the above examples, to force use of the (global, toplevel) postfix :-`=`() operator after such a rebinding of the = symbol? acer
It sounds to me as if you would like to compute the residuals and the sum of squares of the residuals, given a forced least-squares function (model). You could of course compute the sums of squares of the residuals easily in Maple. Simply evaluate the model at each X value, and square the difference with the Y value, and add those up. That's probably one or two lines of Maple code. But you asked whether the Statistics routines could be used for the purpose, in a canned way. You could specify a parameter that doesn't actually appear in your candidate model equation. And then use the solutionmodule output to get the extra details. Here's an example, based on one from the ?Statistics,Fit help-page. Notice that the variable `dummy` doesn't appear in the model, which is the basis of the trick. The model equation is thus fully supplied by the user, with no parameters to actually be estimated (except some incidental estimation of the dummy, which is of no consequence).
> X := Vector([1, 2, 3, 4, 5, 6], datatype=float):
> Y := Vector([2, 3, 4.8, 10.2, 15.6, 30.9], datatype=float):
> Digits:=trunc(evalhf(Digits)):

> # Here is the forced model.

> eq:=0.887576142919275224+0.606352318207692531*exp(0.649251558313310717*t):

> # Here are the simple commands to get the results,
> # without using Statistics.

> add( (eval(eq,t=X[i])-Y[i])^2, i=1..op(1,X) );

                               2.29446465108687
 
> seq( (eval(eq,t=X[i])-Y[i]), i=1..op(1,X) );

0.04819978094862, 0.10913477921538, 0.33987862305974, -1.17305895930964,
 
    0.8671971243695, -0.1913514547231

> # And here is Statistics:-Fit computing those results.

> sol := Statistics:-Fit( eq, X, Y, t, parameternames=[dummy], output=solutionmodule):

> sol:-Results(["residualsumofsquares","residuals"]);

2.29446465108658, [0.0481997809486180984, 0.109134779215386946,
 
    0.339878623059754581, -1.17305895930959636, 0.867197124369345929,
 
    -0.191351454723299952]
acer
Compare these two forms,
> 10 = 0 mod 5;

                                    10 = 0

> (10 = 0) mod 5;

                                     0 = 0
In the first of those, the `mod` operator works only on the 0 and not on the 10. See the help-page ?operators,precedence which may help clarify the relative precedence of various operators. That page shows that mod has a higher binding strength than does =. The round-bracket delimiters in the second example above get around the higher precedence of mod over =. Your code could be written as follows.
printlevel:=2:
for n from 1 to 10 do
if ( (n=0) mod 5 ) then 0;
end if;
end do;
Note that the brackets around the conditional are not strictly necessary, so this should also work here.
printlevel:=2:
for n from 1 to 10 do
if (n=0) mod 5 then 0;
end if;
end do;
acer
Remove the ( round bracket between "if" and "isprime". acer
See ?dsolve[classical] for details. Hopefully I got it right. sys:={diff(y(t),t) = ((2-t)*(2-y(t)))/y(t) , y(0)=1}: sol:=dsolve(sys,numeric,method=classical[foreuler],stepsize=0.2,output=listprocedure): fy[0.2] := eval( y(t), sol): sol:=dsolve(sys,numeric,method=classical[foreuler],stepsize=0.1,output=listprocedure): fy[0.1] := eval( y(t), sol): sol:=dsolve(sys,numeric,method=classical[foreuler],stepsize=0.05,output=listprocedure): fy[0.05] := eval( y(t), sol): fy[0.05](1.0); fy[0.1](1.0); fy[0.2](1.0); plot([fy[0.2],fy[0.1],fy[0.05]],0..8,legend=[fy[0.2],fy[0.1],fy[0.05]]); acer
First 307 308 309 310 311 312 313 Last Page 309 of 322