Carl Love

Carl Love

28095 Reputation

25 Badges

13 years, 100 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

No, your approach is not correct. Your limits of integration are wrong. [1/2 <= x*y <= 2, 1 <= x <= 3] means that the limits of integration should be [y= 1/(2*x)..2/x, x= 1..3]. The integral should be entered as

int(exp(1/(x*y))/(y^2*(x+1)^2), [y= 1/(2*x)..2/x, x= 1..3]);

Note that you should enter this using a single int to avoid an irrelevant warning message.

Note that 1/x*y is not the same as 1/(x*y). Indeed, the former is the same as y/x and the latter is the same as 1/x/y or 1/y/x. So I made it exp(1/(x*y)) because the other forms are not elementary integrals. This parentheses placement is just a guess on my part; you'll need to check it.

Maple's differential operator, D, can find the derivative of quite complicated procedurally defined functions.

f:=x->max(x^2, sqrt(abs(x))):
D(f);

eqn:= diff(y(x),x$4)-20*diff(y(x),x$3)+158*diff(y(x),x$2)-580*diff(y(x),x)+841*y(x);
eval(eqn, y(x)= x*exp(5*x)*cos(2*x));

The above returns 0, which verifies the solution.

The implicitdiff command that you gave works for me (Maple 16.02), so I wonder what version of Maple you are using. Try changing y(x) to y; that also works for me.

g:= (x,y)-> x^3+y^3=1:
implicitdiff(g(x,y), y, x);

Get second derivative:

implicitdiff(%, y, x);

 

 

string1:= cat(seq(sprintf("(%a,%a,)==", pt3[k,1], pt3[k,2]), k= 1..op([1,1], pt3)))[1..-3];

A few characters shorter and a bit more cryptic (but requiring no local variable) is

string1:= sprintf(cat("(%a,%a,)==" $ op([1,1], pt3)), convert(pt3^%T, list)[])[1..-3];

This is not available in Maple. I think that implementing it would require fundamental changes to the GUI.

The problem with you original command is that your didn't use op. Your command was

ptslist := convert([1, 1], points, listlist);

It should be

ptslist:= convert(op([1, 1], points), listlist);

A somewhat more robust command is plottools:-getdata, because it doesn't rely on the "magic numbers" [1,1]:

convert(plottools:-getdata(points, element= "curve")[-1], listlist);

You generally can't rely on a 2d plot command to give you equally spaced data, although some will with some options (such as plot(..., adaptive= false)). 3d plots are generally evenly spaced.

If you want a 3d plot, then you need three Vectors of the same length, one for the x-coordinate data, one for the y-coordinate data, and one for the z-coordinate data. Make sure that you spell Vector with a capital V, unlike in your Question, or that you use the angle-bracket Vector constructor: < 1, 2, 3 >. If you want a 2d plot, then you just need two Vectors. Since you list two independent variables [p, sec], you are specifying a 3d plot.

To combine multiple plots, I usually prefer plots:-display over plots:-multiple, although the latter is sometimes better. Regardless of which you use, you'll need to pass it valid plotting commands. (So, note that display and multiple are both "meta" plotting commands: their primary arguments need to be plotting commands, not algebraic expressions.) The usual command for plotting equations in two variables is plots:-implicitplot. However, when it works, I prefer the command algcurves[plot_real_curve]. This command only works for some polynomial equations, but it usually gives a better plot than implicitplot when both plotting commands are passed minimal options. For example, plot_real_curves doesn't require you to specify ranges for the variables.

eqns:= [x^3-4*x=y, y^3-3*y=x]:

plots:-display(map(algcurves[plot_real_curve], (lhs-rhs)~(eqns), x, y));

From that plot we can see that all the "action" is in the ranges x= -3..3, y= -4..4. This information can be used to refine the plot using implicitplot, as done by Kitonum.

If you substitute numeric values for the variables in an equation, the result is still an equation! Were you perhaps expecting a boolean value? To get such, you need to apply at least the most trivial of Maple's boolean evaluators, evalb. In this case, (assuming that your solutions are exact algebraic numbers such as provided by solve(..., explicit)) you'll also need something like evala to simplify the algebraic numbers so that their equality can be verified.

And what do you mean by "both equations get the same result"? Even if the solution is valid, why should that mean that each equation has the same result? And it is also possible that the solution is invalid and yet the substituted equations are the same. For example, if both equations evaluate to 1=0 that would be an invalid solution, yet both equations are the same. So I see no significance to "both equations get the same result".

Now if what you actually want is to verify the accuracy of the solutions, here's something that you can do:

eqns:= {x^3-4*x=y, y^3-4*y=x}:
solns1:= solve(eqns, [x,y], explicit):
for s in solns1 do evalb~(evala(subs(s, eqns))) end do;

You will get {true} for every correct solution. This is IMO an inelegant method.

A more elegant verification method begins by converting the equations to expressions presumed to be equated to 0. But I will assume that your solution technique (your alternative to my solve) requires them to be equations. So I will convert them to expression form after generating the solutions.

eqns:= {x^3-4*x=y, y^3-4*y=x}:
solns1:= solve(eqns, [x,y], explicit):
Eqns:= (lhs-rhs)~(eqns): #Convert to expression form.
remove(s-> evala(eval(Eqns, s)) = {0}, solns1);

The returned list will contain all the invalid solutions. If they are all valid, that will be the empty list. You can use simplify instead of evala if you wish.

 

 

 

You have numerous syntax errors and some other problems. The specific problem that caused your error message is irrelevant. I am reluctant to list all the errors because there are so many and because some of them interact with some others in ways that make them difficult to describe individually. So first I'll post my corrected version of the code, then I'll make some comments.

B:= 10;
A:= Array(1..B, 0);

#b is in your code, but you never gave it a value.
b:= 2;
#You'll get more meaningful results by using larger primes.
p:= 2^20;

for i from 1 to B do
     p:= nextprime(p); #Not nextprime(i)
     a:= numtheory:-primroot(p);
     #Garbage collecting leads to more accurate timings.
     gc();
     A[i]:= A[i]+
         CodeTools:-Usage(
               numtheory:-mlog(b,a,p, method= indcalc),
               output= cputime, quiet, iterations= 2^12
         )
end do;

A;

Additional comments:

1. Remember to use := for assignment statements rather than simply =
2. Commands in packages need the package prefix, i.e., numtheory:-primroot rather than simply primroot.
3. convert(..., decimal) is meaningless in this context: The values returned by Usage are already decimal.
4. When timing a small operation, you need to use the iterations option to get meaningful results. This is because the resolution of Maple's clock is much, much larger than the resolution of your CPU's clock.
5. I don't see any point in specifying output= [cputime, output] since you're discarding the output.

To get {-2,3}, use

eval({x,y}, result);

If you use sets like that (they're sets because they're in curly braces), then you cannot rely on the first element corresponding to x and the second corresponding to y. Indeed, if the values of x and y are the same, then your set will only have one element. For these reasons, you should use lists instead of sets, like this:

eval([x,y], result);

The $ (unlike seq) can work symbolically, i.e., for symbolic values of i. And diff can work with this information to produce derivatives of arbitrary order. (Of course, diff can only do this for somewhat simple functions.) Observe the output of your first command:

diff(x^7, x$i);

This expression gives a valid derivative for any i, including i = 0. In your first example, this is the expression that is passed to seq, which then simply plugs in the values of i.

The command seq has "special evaluation rules" similar to the first parameter being declared with uneval (see ?uneval) . (These rules also apply to add and mul.) You are correct that these rules are not described at ?seq. These rules mean that in your second example, the unevaluated expression diff(x^7, x$i) is passed to seq. Evaluating this expression at i=0 causes an error because the expression becomes diff(x^7, NULL)i.e., nothing is after the comma. This can be avoided by using diff(x^7, [x$i]), which is valid for any i, even if the [ ] are empty.

 

You have a classic off-by-one error. In addition to the problem pointed out by sazroberson of not propagating in the other arrays, you are also not propagating to the next value of i.

Change the lines

S[i]:= S[i-1], etc. to

S[i+1]:= S[i], etc.

In the long if - elif block, change each line of the form

S[i]:= S[i]+1 to

S[i+1]:= S[i]+1.

Remove the if statement that prevents the execution of this code for i=0.

Jakub,

I have finished the project including extensive research of your course handouts from your professor's website (http://www.kosiorowski.edu.pl/wp-content/uploads/2014/10/Ia-NumMet-W-S-1-Error-theory.pdf) and a complete-from-scratch rewrite of your code including individually checking for zero denominators for each method (so that the procedure will perform all applicable methods even if some are inapplicable). I extended the code so that it works for functions of any number of variables (not just x, y, z) while at the same time vastly simplifying it.



Write the procedure of the name INVERROR, whose entries are:

• 

function of three variables,

• 

approximate values of its variables,

• 

upper bound for value of "f,"

• 

if some method cannot be applied, the information about it should be noted.

 

The procedure should calculate admissible errors for each variable, by three different methods in the inverse problem of error theory.

The exemplary output should be:

=====================================

equal upper bounds method:

e[x] = e[y] and e[y] = e[z] and e[z] = 0.2e-2

equal influence rule:

e[x] = .3, e[y] = 0.2e-2, e[z] = 0.1e-1

equal mesurement method:

e[x] = 0.1e-1, e[y] = 0.5e-2, e[z] = 0.2e-1

=====================================

It may be also in tabluar/vector form. The decision is yours.

 

 

The standard instructions for summation and derivatives do not require any special package.

 

 

errorf:= (f::procedure, X::list(float))->
     trunc(
          10.^Digits*
          evalf(add(abs(D[_](f)(X[])*(round(X[_])-X[_])), _= 1..nops([op(1,eval(f))])))
      )/10.^Digits
;     

proc (f::procedure, X::(list(float))) options operator, arrow; trunc(10.^Digits*evalf(add(abs((D[_](f))(X[])*(round(X[_])-X[_])), _ = 1 .. nops([op(1, eval(f))]))))/10.^Digits end proc

(1)

func := proc (x, y, z) options operator, arrow; 2/x+3/y+4/z end proc;

proc (x, y, z) options operator, arrow; 2/x+3/y+4/z end proc

(2)

Delta[f] := errorf(func, [2.01, 3.11, 2.1]);

.1297720743

(3)

INVERROR:= proc(f::procedure, x::list(float), E::float)
local
     n:= nops([op(1,eval(f))]), #arity of f
     d:= evalf(abs~([seq(D[_](f)(x[]), _= 1..n)])), #abs~grad(f(x))
     X:= abs~(x),
     S:= `+`(d[]) #1-norm of grad(f(x))
;
     if fnormal(S)=0 then
          error "Norm of gradient is 0. Cannot proceed with any method."
     end if;
 
     printf("\nequal upper bounds method:\n%g", E/S);

     if `*`(fnormal~(d)[]) = 0 then
          printf("\nGradient has a zero component. Cannot use equal influence rule.")
     else
          printf("\nequal influence rule:\n%g", <(E/n/~d)[]>)
     end if;

     S:= `+`((d*~X)[]);
     if fnormal(S)=0 then
          printf("\nCannot use equal measurement method due to zero denominator.\n")
     else
          printf("\nequal measurement method:\n%g\n", <(X*~E/~S)[]>)
     end if

end proc:

INVERROR(func, [2.01, 3.11, 2.1], Delta[f]);


equal upper bounds method:
0.0757909
equal influence rule:
0.087382 0.139463 0.0476912
equal measurement method:
0.0674984 0.104438 0.0705207

 

 



Download INVERROR.mw

 

First 277 278 279 280 281 282 283 Last Page 279 of 395