CODEWITHSUNDEEP
pythonsympy
user672139
user672139

Reputation: 11

Writing a function for x * sin(3/x) in python

I have to write a function, s(x) = x * sin(3/x) in python that is capable of taking single values or vectors/arrays, but I'm having a little trouble handling the cases when x is zero (or has an element that's zero). This is what I have so far:

def s(x):
    result = zeros(size(x))
    for a in range(0,size(x)):
        if (x[a] == 0):
            result[a] = 0
        else:
            result[a] = float(x[a] * sin(3.0/x[a]))
    return result

Which...doesn't work for x = 0. And it's kinda messy. Even worse, I'm unable to use sympy's integrate function on it, or use it in my own simpson/trapezoidal rule code. Any ideas?

When I use integrate() on this function, I get the following error message: "Symbol" object does not support indexing.

Upvotes: 0

Views: 908

Answers (2)

asmeurer
asmeurer

Reputation: 91620

If you want to use SymPy's integrate, you need a symbolic function. A wrong value at a point doesn't really matter for integration (at least mathematically), so you shouldn't worry about it.

It seems there is a bug in SymPy that gives an answer in terms of zoo at 0, because it isn't using limit correctly. You'll need to compute the limits manually. For example, the integral from 0 to 1:

In [14]: res = integrate(x*sin(3/x), x)

In [15]: ans = limit(res, x, 1) - limit(res, x, 0)

In [16]: ans
Out[16]:
  9⋅π   3⋅cos(3)   sin(3)   9⋅Si(3)
- ─── + ──────── + ────── + ───────
   4       2         2         2

In [17]: ans.evalf()
Out[17]: -0.164075835450162

Upvotes: 0

ErikR
ErikR

Reputation: 52049

This takes about 30 seconds per integrate call:

import sympy as sp

x = sp.Symbol('x')

int2 = sp.integrate(x*sp.sin(3./x),(x,0.000001,2)).evalf(8)
print int2

int1 = sp.integrate(x*sp.sin(3./x),(x,0,2)).evalf(8)
print int1

The results are:

1.0996940
-4.5*Si(zoo) + 8.1682775

Clearly you want to start the integration from a small positive number to avoid the problem at x = 0.

You can also assign x*sin(3./x) to a variable, e.g.:

s = x*sin(3./x)
int1 = sp.integrate(s, (x, 0.00001, 2))

My original answer using scipy to compute the integral:

import scipy.integrate
import math

def s(x):
  if abs(x) < 0.00001:
    return 0
  else:
    return x*math.sin(3.0/x)

s_exact = scipy.integrate.quad(s,  0, 2)

print s_exact

See the scipy docs for more integration options.

Upvotes: 1

Related Questions