CODEWITHSUNDEEP
pythonmathsympysymbolic-mathintegral
Ryxsen
Ryxsen

Reputation: 15

Harder equations (with Derivatives and Integrals) and ConditionSet in SymPy

I am planning to calculate b (Its also x on Xo axis) for which curve (function) length from 0 to x is equal 1.

By knowing: https://www.mathsisfun.com/calculus/arc-length.html

(integral from 0 to b) ∫ (1 + ((f’(x))^2)^(1/2) dx = 1

and that:

(integral from a to b) ∫ f(x)dx = F(b) - F(a)

We can calculate it by

1 - F(0) + F(b) = 0 , where this is now an Equation in terms of x, beacuse b as I said is x on Xo axis.

So now I tried it for f(x) = x**3 (full code will be below)

F(b) is equal to this monster: https://www.wolframalpha.com/input/?i=integral&assumption=%7B%22C%22%2C+%22integral%22%7D+-%3E+%7B%22Calculator%22%7D&assumption=%7B%22F%22%2C+%22Integral%22%2C+%22integrand%22%7D+-%3E%22%281+%2B+9x%5E4%29%5E%281%2F2%29%22

All I get from SymPy is ConditionSet but not number. Of course ConditionSet cannot be evauluated by evalf()

So here are my questions:

Code:

from __future__ import division
import matplotlib.pyplot as plt
from sympy import *

x, y, z = symbols('x y z', real=True)

function1 = x**3

Antiderivative1 = integrate((1+(diff(function1))**2)**(1/2), x)

b = solveset(Eq(1 + Antiderivative1.subs(x, 0).evalf() - Antiderivative1, 0), x)

print(b)

Thats the output:

ConditionSet(x, Eq(x*hyper((-0.5, 1/4), (5/4,), 9*x**4*exp_polar(I*pi)) - 4.0*gamma(5/4)/gamma(1/4), 0), Complexes)

Thanks in advance and sorry for mistakes in grammar.

Upvotes: 0

Views: 240

Answers (1)

Oscar Benjamin
Oscar Benjamin

Reputation: 14500

Note that you should use S(1)/2 or Rational(1, 2) (or sqrt) rather than 1/2 which will give you a float in Python. With that we have

In [16]: integrand = sqrt(1 + ((x**3).diff(x))**2)                                                                                

In [17]: integrand                                                                                                                
Out[17]: 
   __________
  ╱    4     
╲╱  9⋅x  + 1 

In [18]: antiderivative = integrand.integrate(x)                                                                                  

In [19]: antiderivative                                                                                                           
Out[19]: 
          ┌─  ⎛-1/2, 1/4 │    4  ⅈ⋅π⎞
x⋅Γ(1/4)⋅ ├─  ⎜          │ 9⋅x ⋅ℯ   ⎟
         2╵ 1 ⎝   5/4    │          ⎠
─────────────────────────────────────
               4⋅Γ(5/4)

While that isn't the same form as the result from Wolfram Alpha it could easily be the same function (up to an additive constant). From this result or the one on Wolfram Alpha I very much doubt that you will find an analytic solution (using SymPy or anything else).

You can however find a numerical solution. Unfortunately there is a bug in SymPy's lambdify function that means nsolve doesn't work with this function:

In [22]: nsolve(equation, x, 1)                                                                                                   
...
NameError: name 'exp_polar' is not defined

We can do it ourselves with Newton steps though:

In [76]: f = equation.lhs                                                                                                         

In [77]: fd = f.diff(x)                                                                                                           

In [78]: newton = lambda xi: (xi - f.subs(x, xi)/fd.subs(x, xi)).evalf()                                                          

In [79]: xj = 1.0                                                                                                                 

In [80]: xj = newton(xj); print(xj)                                                                                               
0.826749667942050

In [81]: xj = newton(xj); print(xj)                                                                                               
0.791950624620750

In [82]: xj = newton(xj); print(xj)                                                                                               
0.790708415511451

In [83]: xj = newton(xj); print(xj)                                                                                               
0.790706893629886

In [84]: xj = newton(xj); print(xj)                                                                                               
0.790706893627605

In [85]: xj = newton(xj); print(xj)                                                                                               
0.790706893627605

Upvotes: 1

Related Questions