home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #18 / NN_1992_18.iso / spool / sci / math / 10348 < prev    next >
Encoding:
Internet Message Format  |  1992-08-17  |  22.6 KB

  1. Xref: sparky sci.math:10348 sci.math.symbolic:2226 sci.math.num-analysis:2479
  2. Path: sparky!uunet!gatech!emory!riddle
  3. From: riddle@mathcs.emory.edu (Larry Riddle)
  4. Newsgroups: sci.math,sci.math.symbolic,sci.math.num-analysis
  5. Subject: Re: Some comparison of Mma, Maple and SymbMath (REBUTTAL)
  6. Keywords: Mathematica, Maple, SymbMath
  7. Message-ID: <9315@emory.mathcs.emory.edu>
  8. Date: 18 Aug 92 03:37:38 GMT
  9. Followup-To: sci.math.symbolic
  10. Organization: Emory University, Dept of Math and CS
  11. Lines: 859
  12.  
  13. In a previous article, Weiguang Huang writes:
  14.  
  15. >    Here are some problems that Maple and Mathematica cannot 
  16. >solve, but SymbMath can do. 
  17. >    The following examples came from news on the sci.math.symbolic 
  18. >newsgroup in 1991, and were run in Maple V, Mathematica 2.0, or  
  19. >SymbMath 2.1. 
  20.  
  21. Unfortunately, some aspects of this article were inaccurate and
  22. misleading. Since I do not use Maple, I will only comment on the
  23. examples with respect to Mathematica.
  24.  
  25. Summary of comments:
  26. (1) Of the 10 examples that do not involve improper integrals or
  27. differential equations, Mathematica can do every example once we define
  28. some simple rules or load some standard packages.
  29.  
  30. (2) Version 2.1 of Mathematica can solve the 2 differential equations
  31. examples using the new package Calculus'DSolve'.
  32.  
  33. (3) SymbMath appears to evaluate improper integrals by evaluating the
  34. indefinite integral at the endpoints. This correctly gives the Cauchy
  35. Principal Value in some cases, but is *not* the way that the CPV is
  36. defined. Simple translations of these "correct" examples produce the
  37. same incorrect answers as given by Matheamtica. In several examples,
  38. Huang breaks the integral into two parts at the singularity to get
  39. SymbMath to give the correct answer, but *does not* do the same thing
  40. with Mathematica, then claims that SymbMath can do the example but
  41. Mathematica cannot. When SymbMath is told to compute the integral over
  42. *one* interval, it gives the same incorrect answer as Matheamtica.
  43.  
  44. Conclusion: With regard to these 20 examples, there is no difference
  45. between Mathematica and SymbMath.
  46.  
  47. Disclaimer: I did not write any part of Mathematica. In fact, I'm a
  48. rather novice Mathematica user. I think Huang has done an impressive
  49. job with SymbMath. I just think some misconceptions need to be corrected.
  50.  
  51. >************************ Example 1 ******************************
  52. >> int(exp(-a * x^2), x=0..infinity);
  53.  
  54. >    SymbMath :
  55. >    Input:
  56. >inte(exp(-a*x^2), x from 0 to inf)
  57. >assume(sqrt(a) > 0)
  58. >inte(exp(-a*x^2), x from 0 to inf)
  59. >    Output:
  60. >1/2*a^(-0.5)*sqrt(pi)*erf(inf*sgn(sqrt(a)))
  61. >assumed
  62. >1/2*a^(-0.5)*sqrt(pi)
  63. >
  64.  
  65. Mathematica 2.0 for SPARC
  66. Copyright 1988-91 Wolfram Research, Inc.
  67.  -- NeWS graphics initialized --
  68.  
  69. In[1]:= Integrate[Exp[-a x^2],{x,0,Infinity}]
  70.  
  71.         Sqrt[Pi]
  72. Out[1]= ---------
  73.         2 Sqrt[a]
  74.  
  75. Well, this is really only true if a is positive. But let's take a
  76. look at SymbMath's answer:
  77.  
  78. 1/2*a^(-0.5)*sqrt(pi)*erf(inf*sgn(sqrt(a)))
  79.  
  80. What does sgn(sqrt(a)) mean if a is not positive? Is this answer really
  81. any better than Mathematica's answer? In fact, if we integrate from 0
  82. to b, rather than from 0 to Infinity, Mathematica gives
  83.  
  84. In[2]:= Integrate[Exp[-a x^2],{x,0,b}]
  85.  
  86.         Sqrt[Pi] Erf[Sqrt[a] b]
  87. Out[2]= -----------------------
  88.                2 Sqrt[a]
  89.  
  90. which is exactly SymbMath's answer with b = inf. Thus I see no
  91. difference between Mathematica and SymbMath here, with slight
  92. preference to Maple's requirement of demanding to know if a is positive
  93. before returning an answer.
  94.  
  95. >************************** Example 2 *********************************
  96. >> int(x^n, x=0..1);
  97. >
  98. >    SymbMath :
  99. >    Input:
  100. >assume(n > -1)
  101. >inte(x^n, x from 0 to 1)
  102. >    Output:
  103. >assumed
  104. >1/(1 + n)
  105. >
  106.  
  107. Mathematica:
  108. In[3]:= Integrate[x^n,{x,0,1}]
  109.  
  110.                  1 + n
  111.           1     0
  112. Out[3]= ----- - ------
  113.         1 + n   1 + n
  114.  
  115. Ok, this looks a bit strange because Mathematica has left the power of
  116. 0 unevaluated. That is perhaps not surprising because we have made no
  117. claims about what type of number n is. 
  118.  
  119. What would SymbMath give if there were no assumptions about n? (Notice
  120. that Huang made the assumption n > -1 *before* he evaluated the
  121. integral.)
  122.  
  123. SymbMath:
  124. Input:
  125. inte(x^n, x from 0 to 1)
  126. Output:
  127. -0^sgn(1+n)/(1+n) + 1/(1+n)
  128.  
  129. Why, it looks just like Mathematica's answer!
  130.  
  131. This is easy to take care of in Matheamtica.
  132. First we load the package Declare.m written by Pekka
  133. Jahunnen that will allow us to declare n to be positive.
  134.  
  135. In[4]:= <<Declare.m
  136. {Declare, NewDeclare, NonPositive, RealQ}
  137.  
  138. In[5]:= Declare[n,Positive]
  139.  
  140. Now we must tell Mathematica how to evaluate positive powers of 0.
  141. Again, easy to do.
  142.  
  143. In[6]:= Unprotect[Power]
  144.  
  145. Out[6]= {Power}
  146.  
  147. In[7]:= 0^x_?Positive = 0
  148.  
  149. Out[7]= 0
  150.  
  151. In[8]:= Integrate[x^n,{x,0,1}]
  152.  
  153.           1
  154. Out[8]= -----
  155.         1 + n 
  156.  
  157. The same answer as SymbMath which also had to declare n > -1. The only
  158. extra we really had to add in Mathematica was that 0^n = 0 when n is
  159. positive, which was trivial to do.
  160.  
  161. >************************** Example 3 *****************************
  162. >> int(x^n, x=eps..1);
  163. >
  164. >                                          (n + 1)
  165. >                                 1     eps
  166. >                               ----- - ----------
  167. >                               n + 1      n + 1
  168. >
  169. >> limit(", eps=0, right);
  170. >
  171. >    SymbMath:
  172. >    Input:
  173. >assume(n > -1)
  174. >inte(x^n, x from eps to 1)
  175. >subs(eps=0 to last)
  176. >    Output:
  177. >assumed
  178. >1/(1 + n) - eps^(1 + n)/(1 + n)
  179. >1/(1 + n)
  180.  
  181. Mathematica:
  182. In[9]:=Declare[n,Positive]
  183.  
  184. In[10]:= Integrate[x^n, {x,eps,1}]
  185.  
  186.                     1 + n
  187.            1     eps
  188. Out[10]= ----- - --------
  189.          1 + n    1 + n
  190.  
  191. In[11]:= Limit[%,eps->0]
  192.  
  193.            1
  194. Out[11]= -----
  195.          1 + n
  196.  
  197. Again, the exact same answers as given by SymbMath. We are still using
  198. the rules defined for powers of 0 and that n is positive, but didn't
  199. have to do anything special for this example once those rules were made.
  200.  
  201. >*************************** Example 4 ***************************
  202. >> 0^n;
  203. >            0
  204. >
  205. ># Maple flags the error only when 'n' is replaced by the constant 0:
  206. >
  207. >> 0^0;
  208. >Error, 0^0 is undefined
  209. >
  210. >    SymbMath:
  211. >    Input:
  212. >assume(n > 0)
  213. >0^n
  214. >0^-n
  215. >0^0
  216. >    Output:
  217. >assumed
  218. >0
  219. >discont
  220. >undefined
  221. >
  222.  
  223. Huang complains about Maple's answer for 0^n when no assumption is made
  224. about n, but then assumes n > 0 before showing what SymbMath gives for
  225. 0^n. In fact, with no assumption about n, SymbMath will answer
  226. 0^sgn(n) for 0^n. 
  227.  
  228. Mathematica gives
  229. In[12]:= Declare[n, Positive]
  230.  
  231. In[13]:= 0^n  (* Remember we defined a rule to handle this *)
  232.  
  233. Out[13]= 0
  234.  
  235. In[14]:= 0^0
  236.  
  237.                                         0
  238. Power::indet: Indeterminate expression 0  encountered.
  239.  
  240. Out[14]= Indeterminate 
  241.  
  242. No difference between Mathematica and SymbMath here. I could easily
  243. have defined the rule 0^x_?Negative = Infinity and then Mathematica
  244. would have also returned 0^(-n) = Infinity.
  245.  
  246. >**************************** Example 5 ******************************
  247. >    Maple:
  248. >> int(x^k, x);
  249. >
  250. >                                     (k + 1)
  251. >                                    x
  252. >                                    --------
  253. >                                      k + 1
  254. >
  255. >    SymbMath:
  256. >    Input:
  257. >inte(x^k*d(x))
  258. >subs(k=-1 to last)
  259. >    Output:
  260. >constant + x^(1 + k)/(1 + k)
  261. >discont
  262.  
  263. I'm not really sure what the point of this example was. What did Maple
  264. do "wrong" that SymbMath did better? Anyway, here is Mathematica:
  265.  
  266. In[15]:= Integrate[x^k,x]
  267.  
  268.           1 + k
  269.          x
  270. Out[15]= ------
  271.          1 + k
  272.  
  273. In[15]:= % /. k -> -1
  274.  
  275.                                  1
  276. Power::infy: Infinite expression - encountered.
  277.                                  0
  278.  
  279. Out[16]= ComplexInfinity 
  280.  
  281. Seems Mathematica can do this example also.
  282.  
  283. >**************************** Example 6 *****************************
  284. >> limit(x^k/exp(x), x=infinity);
  285. >
  286. >    SymbMath:
  287. >    Input:
  288. >lim(x=inf, x^k/exp(x))
  289. >lim(x=inf, x^(10^10)/exp(x))
  290. >lim(x=inf, x^(10^10000)/exp(x))
  291. >    Output:
  292. >0
  293. >0
  294. >0
  295.  
  296. Mathematica does have trouble doing these limits for k > 3 with the
  297. built-in Limit command. However, if we load the package Calculus'Limit'
  298. then Mathematica can do this example also.
  299.  
  300. In[17]:= <<Limit.m
  301.  
  302. In[18]:= Limit[x^k/Exp[x],x->Infinity]
  303.  
  304. Out[18]= 0
  305.  
  306. In[19]:= Limit[x^(10^10)/Exp[x],x->Infinity]
  307.  
  308. Out[19]= 0
  309.  
  310. In[20]:= Limit[x^(10^1000)/Exp[x],x->Infinity]
  311.  
  312. Out[20]= 0
  313.  
  314. >****************************** Example 7 ****************************
  315. >    Maple:
  316. >> int(x^m * exp (-b * x), x=0..infinity);
  317. >
  318. >                            infinity
  319. >                               /
  320. >                              |       m
  321. >                              |      x  exp(- b x) dx
  322. >                              |
  323. >                             /
  324. >                             0
  325. >
  326. >
  327. ># As expected, the integral remains unevaluated.  Now declare the signs 
  328. ># of the symbolic parameters:
  329. >
  330. >> signum(b) := 1;
  331. >
  332. >                                 signum(b) := 1
  333. >
  334. >> signum(m) := 1;
  335. >
  336. >                                 signum(m) := 1
  337. >
  338. >
  339. ># Upon attempting to compute the integral a second time...
  340. >
  341. >> int(x^m * exp (-b * x), x=0..infinity);
  342. >
  343. >                            infinity                        
  344. >                               /                            
  345. >                              |       m                     
  346. >                              |      x  exp(- b x) dx       
  347. >                              |                             
  348. >                             /
  349. >                             0                              
  350. >
  351. ># ...THE INTEGRAL CONTINUES TO REMAIN UNEVALUATED, despite the fact that
  352. ># the signs of the parameters 'b' and 'm' were declared PRIOR to the second
  353. ># attempt at computation.
  354. >
  355. >    SymbMath:
  356. >    Input:
  357. >inte(x^n*exp(-a*x), x from 0 to inf)
  358. >    Output:
  359. >inte(x^n*exp(-a*x), x, 0, inf)
  360.  
  361. Huh? It seems SymbMath cannot do this integral either.
  362. Mathematica:
  363.  
  364. In[21]:= Integrate[x^m Exp[-b x],{x,0,Infinity}]
  365.  
  366.           -1 - m
  367. Out[21]= b       Gamma[1 + m]
  368.  
  369. >************************ Example 8 *********************************
  370. >    Mathematica:
  371. >In[1]:= Integrate[1/x,{x,-1,1}]
  372. >Out[1]= -Log[-1]
  373. >
  374. >    Maple:
  375. >has the same problem.
  376. >
  377. >    SymbMath:
  378. >    Input:
  379. >inte(1/x, x from -1 to 1)
  380. >inte(1/x, x from -1 to 2)
  381. >    Output:
  382. >0
  383. >ln(2)
  384. >
  385. The problems of Mathematica with improper integrals with an
  386. internal singularity are well known. But is SymbMath's answer any
  387. better? The standard definition of the inproper integral of f(x) from a
  388. to b when f has a singularity at c between a and b is
  389.  
  390. limit int(f[x],a to c-e1) + limit int(f[x], c+e2 to b)   (1)
  391. e1->0+                       e2->0+
  392.  
  393. and with this definition int(1/x, x from -1 to 1) does not exist, so
  394. SymbMath's answer is wrong. What SymbMath has computed is actually the
  395. Cauchy Principal Value, defined by
  396.  
  397. limit [ int(f[x],a to c-e) + int(f[x],c+e to b} ]        (2)
  398. e->0+
  399.  
  400. If the limit in (1) exists, then so does the limit in (2), but the
  401. converse is not true. Which definition to use is, of course, up to the
  402. user and depends on the reason for computing the integral. However, the
  403. first definition is certainly the one most often encountered,
  404. particularly in undergraduate calculus courses. In fact, Watson Fulks
  405. writes in "Advanced Calculus: An Introduction to Analysis" that with
  406. the Cauchy Principal Value "we are not dealing with a proper integral
  407. nor even with an ordinary or garden-variety improper one" (p419).
  408.  
  409. Mathematica *can* compute Cauchy Principal Values also, but the user
  410. must supply the exact location of the singularity. Here are the two
  411. examples that Huang gives:
  412.  
  413. In[4]:= <<CauchyPrincipalValue.m
  414.  
  415. In[5]:= CauchyPrincipalValue[1/x,{x,-1,{0},1}]
  416.  
  417. NIntegrate::ploss:
  418.    Numerical integration stopping due to loss of precision. Achieved neither
  419.      the requested PrecisionGoal nor AccuracyGoal; suspect one of the
  420.      following: highly oscillatory integrand or the true value of the integral
  421.      is 0.
  422.  
  423. Out[5]= 0.
  424.  
  425. In[6]:= CauchyPrincipalValue[1/x,{x,-1,{0},2}]
  426.  
  427. NIntegrate::ploss:
  428.    Numerical integration stopping due to loss of precision. Achieved neither
  429.      the requested PrecisionGoal nor AccuracyGoal; suspect one of the
  430.      following: highly oscillatory integrand or the true value of the integral
  431.      is 0.
  432.  
  433. Out[6]= 0.693147
  434.  
  435. >*************************** Example 9 *******************************
  436. >    Maple:
  437. >has a problem for int(tan(x), x=0..pi).
  438. >
  439. >    SymbMath:
  440. >    Input:
  441. >inte(tan(x), x from 0 to pi)
  442. >    Output:
  443. >0
  444.  
  445. Again, SymbMath has computed the Principal Value. Under the first
  446. definition of an improper integral given above, this integral does not
  447. exist.
  448.  
  449. In Mathematica:
  450. In[2]:= CauchyPrincipalValue[Tan[x],{x,0,{Pi/2},Pi}]
  451.  
  452. NIntegrate::ploss:
  453.    Numerical integration stopping due to loss of precision. Achieved neither
  454.      the requested PrecisionGoal nor AccuracyGoal; suspect one of the
  455.      following: highly oscillatory integrand or the true value of the integral
  456.      is 0.
  457.  
  458.                    -16
  459. Out[2]= -1.11022 10   
  460.  
  461. Well, this looks close enough to 0 for me.
  462.  
  463. >**************************** Example 10 ****************************
  464. >    Mathematica:
  465. >cannot evaluate integral of sgn(x).
  466. >
  467. This is incorrect
  468.  
  469. In[3]:= Integrate[Sign[x],{x,-1,1}]
  470.  
  471. Out[3]= 0
  472.  
  473. In[4]:= Integrate[Sign[x],{x,-1,2}]
  474.  
  475. Out[4]= 1 
  476.  
  477. >    SymbMath:
  478. >    Input:
  479. >inte(sgn(x), x from -1 to 1)
  480. >inte(sgn(x), x from -1 to 2)
  481. >    Output:
  482. >0
  483. >1
  484. >
  485. >************************* Example 11 ********************************
  486. >Implicit diff. gives 1+y'[x](1+1/y[x])==0; y'[x]==-y[x]/(y[x]+1).
  487. >
  488. >    Mathematica:
  489. > given this eq. as input to DSolve says that built-in 
  490. > procedure can't solve it.
  491. >
  492. DSolve in version 2.0 cannot solve this equation. However, version 2.1
  493. comes with an additional DSolve package that *can* solve this equation.
  494. From the Macintosh 2.1 version we get
  495.  
  496. In[1]:=DSolve[y'[x]==-y[x]/(y[x]+1),y[x],x]
  497.  
  498. Solve::tdep:
  499.  
  500.   The equations appear to involve trancendental functions
  501.     of the variables in an essentially non-algebraic way.
  502.  
  503. Out[1]= Solve[-Log[y[x]] - y[x] == x + C[1] , y[x]]
  504. (Note, what this means is that Mathematica was unable to find an 
  505. *explicit* solution for y[x], but neither could SymbMath. They just
  506. indicated this in different ways)
  507.  
  508. >
  509. >    SymbMath:
  510. >solve the differential equation by integration inte() or by dsolve().
  511. >    Input:
  512. >d(y)/d(x)*(1+1/y) === -1
  513. >inte(last*d(x))
  514. >Expand=On
  515. >dsolve(d(y)/d(x) === -y/(y+1), y)
  516. >    Output:
  517. >(1 + 1/y)*d(y)/d(x) === -1
  518. >y + ln(y*sgn(y)) === constant - x
  519. >Expand = On
  520. >-y - ln(y*sgn(y)) === constant + x
  521. >
  522. >
  523. >************************* Example 12 ********************************
  524. >    Mathematica:
  525. >                 y'[x] = y[x]^(1/2)
  526. >                  y[0] = 0
  527. >
  528. >DSolve could not handle it. (It rarely solves anything!!), and the 
  529. >RungeKutta package only gave me the solution
  530. >             
  531. >                    y[x]=0
  532. >
  533. >Obviously, there is another solution viz.
  534. >
  535. >                   y[x] = (x/2)^2   
  536.  
  537. Again, version 2.1 can solve this equation:
  538.  
  539. In[2]:= DSolve[{y'[x]==y[x]^(1/2), y[0]==0},y[x],x]
  540.  
  541.                  2
  542.                 x
  543. Out[2]= {{y[x]->--}}
  544.                 4
  545.  
  546. >
  547. >    SymbMath:
  548. >    Input:
  549. >dsolve(d(y)/d(x) === sqrt(y), y) [Can SymbMath handle initial conditions?]
  550. >(last/2)^2
  551. >    Output:
  552. >2*sqrt(y) === constant + x
  553. >y === 1/4*(constant + x)^2
  554. >
  555. >************************* Example 13 ********************************
  556. >    Mathematica:
  557. > Sqrt[a^2] evaluates to Sqrt[a^2].
  558. >
  559. >    SymbMath:
  560. >    Input:
  561. >sqrt(x^2)
  562. >assume(a > 0)
  563. >sqrt(a^2)
  564. >assume(b <0 )
  565. >sqrt(b^2)
  566. >    Output:
  567. >x*sgn(x)
  568. >assumed
  569. >a
  570. >assumed
  571. >-b
  572. >
  573.  
  574. If we tell Mathematica what kind of variables we are working with, it
  575. has no problem with this example:
  576.  
  577. In[4]:= <<Declare.m
  578. {Declare, NewDeclare, NonPositive, RealQ}
  579.  
  580. In[5]:= Declare[x,Real,a,Positive,b,Negative]
  581.  
  582. In[6]:= Sqrt[x^2]
  583.  
  584. Out[6]= Abs[x]
  585.  
  586. In[7]:= Sqrt[a^2]
  587.  
  588. Out[7]= a
  589.  
  590. In[8]:= Sqrt[b^2]
  591.  
  592. Out[8]= -b  
  593.  
  594. >********************** Example 14 **********************************
  595. >    Maple and Mathematica cannot find the integrals of abs(x).
  596. >
  597. >    SymbMath:
  598. >    Input:
  599. >inte(abs(x), x from -1 to 1)
  600. >inte(abs(x)^5*d(x))
  601. >    Output:
  602. >1
  603. >constant + 1/6*abs(x)^6*sgn(x)
  604. >
  605.  
  606. Huang is correct. But it's not hard to fix this. Let's define our
  607. own version of abs by abs[x] = x Sign[x] 
  608.  
  609. In[1]:= abs[x_] := x Sign[x]
  610.  
  611. General::spell1:
  612.    Possible spelling error: new symbol name "abs"
  613.      is similar to existing symbol "Abs".
  614.  
  615. In[2]:= Integrate[abs[x],{x,-1,1}]
  616.  
  617. Out[2]= 1   (* Same answer as SymbMath *)
  618.  
  619. In[3]:= Integrate[abs[x]^5,x]
  620.  
  621.                    5        5
  622. Out[3]= Integrate[x  Sign[x] , x]   (* Integration failed *)
  623.  
  624. In[4]:= Unprotect[Power,Integrate]
  625.  
  626. Out[4]= {Power, Integrate}
  627.  
  628. In[5]:= Sign[x_]^n_?EvenQ = 1    (* Give rules for powers of Sign[x] *)
  629.  
  630. Out[5]= 1
  631.  
  632. In[6]:= Sign[x_]^n_?OddQ = Sign[x]
  633.  
  634. Out[6]= Sign[x]
  635.  
  636. In[7]:= Integrate[x_^n_?OddQ Sign[x_], x_] = x^(n+1)/(n+1) Sign[x]
  637.  
  638.          1 + n
  639.         x      Sign[x]
  640. Out[7]= --------------   (* Give rule for integration of odd power *)
  641.             1 + n
  642.  
  643. In[8]:= Integrate[abs[x]^5,x]
  644.  
  645.          6
  646.         x  Sign[x]
  647. Out[8]= ----------  (* Same answer as SymbMath *)
  648.             6
  649.  
  650. In[9]:= Integrate[abs[x]^4,x]
  651.  
  652.          5
  653.         x
  654. Out[9]= --
  655.         5
  656.           
  657. We had to define our own rules, but Mathematica could do this example.
  658. Of course, I suspect that SymbMath has similar rules built-in.
  659.  
  660. >---------------------------------------------------------------------
  661. >
  662. >    The following problems are taken from Swokowski's Calculus 
  663. >book.  They cause Mathematica to fail because of a singularity in 
  664. >the interior of the interval of integration:
  665. >
  666. >        Section 10.4 Problems 3, 12, 15, 16, 23, 29.
  667. >
  668. >The comments "INTEGRAL IS DIVERGENT" and "Principal Value" come from
  669. >Macsyma.  Mma gives no indication that anything is amiss.
  670. >
  671. >************************ Problem3 *********************************
  672. >    Mathematica:
  673. >In[4]:= Integrate[1/x^2,{x,-3,1}]
  674. >
  675. >          4
  676. >Out[4]= -(-)                            (* INTEGRAL IS DIVERGENT *)
  677. >          3
  678. >
  679. >    SymbMath:
  680. >    Input:
  681. >inte(1/x^2, x from -3 to 1)
  682. >    Output:
  683. >inf
  684.  
  685. This is impressive. Was SymbMath able to detect the singularity and
  686. compute the appropriate improper integral? Let's translate the
  687. function and interval one unit to the right. This should not affect the
  688. value of the integral
  689.  
  690.   SymbMath:
  691.   Input:
  692. inte(1/(x-1)^2, x from -2 to 2)
  693.   Output:
  694. -4/3
  695.  
  696. Oops. Looks like SymbMath got the same answer as Mathematica.
  697.  
  698. >
  699. >************************** Problem12 ***************************
  700. >    Mathemtica:
  701. >In[6]:= Integrate[x^(-4/3),{x,-1,1}]
  702. >
  703. >Out[6]= -6                              (* INTEGRAL IS DIVERGENT *)
  704. >
  705. >    SymbMath:
  706. >    Input:
  707. >inte(x^(-4/3), x from -1 to 1)
  708. >        Output:
  709. >inf
  710. >
  711. Let's try translating the function and interval again. 
  712.  
  713.   SymbMath:
  714.   Input:
  715. inte((x-1)^(-4/3), x from 0 to 2)
  716.   Output:
  717. -6
  718.  
  719. Well, again SymbMath got the same incorrect answer as Mathematica.
  720. My suspicion is that SymbMath has built-in rules to tell it how to
  721. integrate x^(-n) from a to b when n > 1 and a < 0 < b. However, when
  722. SymbMath cannot use one of these rules, as when we translated the
  723. function, then it simply integrates and evaluates at the endpoints,
  724. *just as Mathematica does*.
  725.  
  726. >**************************** Problem15  ************************
  727. >    Mathematica:
  728. >In[7]:= Integrate[1/x,{x,-1,2}]
  729. >
  730. >Out[7]= -Log[-1] + Log[2]
  731. >
  732. >    Maple:
  733. >has the same problem.
  734. >    
  735. >    Macsyma:
  736. >    (c7) integrate(1/x,x,-1,2);
  737. >    Principal Value
  738. >    (d7)                                log(2)
  739. >
  740. >    SymbMath:
  741. >    Input:
  742. >inte(1/x, x from -1 to 2)
  743. >    Output:
  744. >ln(2)
  745. >
  746. Here SymbMath has computed the Cauchy Principal Value. The improper
  747. integral as Swokowski would define it is divergent.
  748.  
  749. >
  750. >********************** Problem16 ************************************
  751. >    Mathematica:
  752. >In[8]:= Integrate[1/(x^2-x-2),{x,0,4}]
  753. >
  754. >        -Log[-2]   Log[2]   Log[5]
  755. >Out[8]= -------- + ------ - ------
  756. >           3         3        3
  757. >
  758. >    Macsyma:
  759. >    (c8) integrate(1/(x^2-x-2),x,0,4);
  760. >    Principal Value
  761. >                                         log(5)
  762. >    (d8)                               - ------
  763. >                                           3
  764. >
  765. >
  766. >    SymbMath:
  767. >    Input:
  768. >inte(1/(x^2-x-2), x from 0 to 4)
  769. >    Output:
  770. >-1/3*ln(2) + 1/3*ln(2/5)
  771. >
  772. Here is an interesting answer. Why would it be of this form? It might
  773. be interesting to see how SymbMath would do the indefinite integral.
  774.  
  775.   SymbMath:
  776.   Input:
  777. inte(1/(x^2-x-2)*d(x))
  778. subs(x=4 to last) - subs(x=0 to last)
  779.   Output:
  780. constant + 1/3*ln(abs(-4+2*x)/abs(2+2*x))
  781. -1/3*ln(2) + 1/3*ln(2/5)
  782.  
  783. It appears that SymbMath evaluated this improper integral by simply
  784. integrating and evaluating at the endpoints, which is not,
  785. in general, a valid procedure. It works here because if
  786.  
  787. int(f(x)*d(x)) = ln(abs(g(x))  where g(x) is continuous 
  788.  
  789. then the Cauchy Principal Value between a and b will equal
  790. ln|g(b)| - ln|g(a)|. So the one advantage SymbMath has over Mathematica
  791. in this context is that it has been programmed to use
  792.  
  793. inte(1/x*d(x)) = ln (x * sgn(x))
  794.  
  795.  
  796. >***************************** Problem23 **********************
  797. >    Mathematica:
  798. >In[10]:= Integrate[(1/x^2)Cos[1/x],{x,-1,2}]
  799. >
  800. >                       1
  801. >Out[10]= Sin[-1] - Sin[-]              (* INTEGRAL IS DIVERGENT *)
  802. >               2
  803. >
  804. >    SymbMath:
  805. >    Input:
  806. >y=1/x^2*cos(1/x)
  807. >inte(y, x from -1 to 0-zero) + inte(y, x from 0+zero to 2)
  808. >    Output:
  809. >y = x^(-2)*cos(1/x)
  810. >sin(-1) - sin(1/2) + 2*sin(inf)
  811. >
  812. Now wait a minute. Huang makes Mathematica integrate over an interval
  813. that includes the singularity, but has SymbMath break the integral into
  814. its two part as the definition of this improper integral requires. That
  815. doesn't seem quite fair. Let's have Mathematica do the integral the
  816. same way that SymbMath does:
  817.  
  818. In[1]:= Integrate[1/x^2 Cos[1/x],{x,-1,0}] + Integrate[1/x^2 Cos[1/x],{x,0,Infinity}]
  819.  
  820.                                  1
  821. Power::infy: Infinite expression - encountered.
  822.                                  0
  823.  
  824. Out[1]= Infinity
  825.                
  826. Now Mathematica tells us that the integral is divergent. It would have
  827. been interesting to see what answer SymbMath gave if we asked it to
  828. integrate from -1 to 2, but the shareware version of SymbMath, which I
  829. was using, does not include the trig functions.
  830.  
  831. >************************** Problem29 ********************************
  832. >    Mathematica:
  833. >In[12]:= Integrate[1/(x-4)^2,{x,0,Infinity}]
  834. >
  835. >           1
  836. >Out[12]= -(-)                          (* INTEGRAL IS DIVERGENT *)
  837. >           4
  838. >
  839. >    SymbMath:
  840. >    Input:
  841. >y=1/(x-4)^2
  842. >inte(y, x from 0 to 4-zero) + inte(y, x from 4+zero to inf)
  843. >    Output:
  844. >y = (-4 + x)^(-2)
  845. >inf
  846. >
  847. Here we go again, asking Mathematica to integrate over the singularity
  848. while letting SymbMath integrate over two intervals. Let's have
  849. SymbMath do it the same way as Mathematica:
  850.  
  851.   SymbMath:
  852.   Input:
  853. inte(1/(x-4)^2, x from 0 to inf)
  854.   Output:
  855. -1/4
  856.  
  857. The same *wrong* answer as Mathematica. Now let's have Mathematica do
  858. the integral over two intervals:
  859.  
  860. In[2]:= Integrate[1/(x-4)^2,{x,0,4}] + Integrate[1/(x-4)^2,{x,4,Infinity}]
  861.  
  862. Infinity::indet: Indeterminate expression -Infinity + Infinity encountered.
  863.  
  864. Out[2]= Indeterminate
  865.  
  866. Mathematica now correctly determines that the integral is divergent.
  867. -- 
  868. Larry Riddle        | riddle@mathcs.emory.edu         PREFERRED
  869. Agnes Scott College | {rutgers,gatech}!emory!riddle   UUCP 
  870. Dept of Math        | riddle@emory.bitnet             NON-DOMAIN BITNET
  871. Decatur, GA 30030   | Phone: Voice 404-371-6222, FAX 404-371-6177
  872.