7

How accurate do you want the answer ? How costly is evaluating your function ? If it is costly, then you dont want to use a rule with too many nodes. How many times do you want to do the quadrature ? If it has to be repeated many times, like in finite element methods, the cost of computing weights/nodes of gauss rule will not matter much since it has to be ...


5

There is no need for numerical computation here. First, $T(q)$ is a well-known function, the logarithmic integral. Repeated integration by parts gives an asymptotic expansion $$\mathrm{Li}(q) = \frac{q}{\log q}\sum_{k=0}^{K-1} \frac{k!}{\log^k q} + O\left(\frac{q}{\log^{K+1}q}\right).$$ There's also a fairly rapidly convergent representation due to ...


5

The calculation of $$ \int_{-1}^{1} L_k(x)\,\text{d} x $$ for the Lagrange polynomials $L_k$ defined on an arbitrary grid $x_k, k=0,\ldots,n$ can be performed by the following two steps: Calculate the Clenshaw-Curtis quadrature weights $w^{\text{cc}}_k$ on the Chebyshev extrema grid $y_k$ for $k=0,\ldots,n$: $$ y_k = \cos\left(\frac{k\pi}{n}\right)\\ w^{\...


4

This integral has a closed analytic solution. The trick is to write $$\frac{1}{x^3+ ax + 2a} = \frac{A}{x-x_1} + \frac{B}{x-x_2} + \frac{C}{x-x_3}$$ by a method called partial fraction decomposition. The values $x_i$ are the roots of the polynomial in the denominator. As your polynomial only has third degree, an analytic formula for the roots exists (aka "...


4

You have to write your second order equation as a system of two first order equations. Let $y' = v$, then your equation $$ y'' + \omega^2 y = 0 $$ becomes $$ \begin{pmatrix} y' \\ v' \end{pmatrix} = \begin{pmatrix} v \\ -\omega^2 y \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -\omega^2 & 0 \end{pmatrix} \begin{pmatrix} y \\ v \end{pmatrix} $$ If you ...


4

Your principal assumption is wrong: But on standing foot, nobody knows what wi and xi are. They need to be calculated. That takes time. It's just not true that "nobody knows" what the weights and points are: They are tabulated and don't need to be computed any more. We also don't use 4096 points -- we subdivide the interval of integration into smaller ...


4

You are almost there, just put $t$ under the square root, and it will become $ \displaystyle\int_0^{1/1095} {\frac{dt}{\sqrt{(1+644.153t)(4.17 \cdot 10^{-5} + 0.145 t)}}}, $ which eliminates the singularity at $t=0$


4

Is there a way to solve for multiple initial conditions without resorting to a slow python for loop? I think the solution to your problem might be parallelization, not vectorization, unfortunately. While Python is going to be slower than C or C++, the for loop in question isn't actually a part of the computation, it's a different problem numerically because ...


4

You write $S(q)$ and $T(q)$ as integrals, but it is easier to think of them as solutions of ODEs: $$ S'(q) = \sin^2\left(\frac{π\Gamma(q)}{2q}\right) $$ with initial conditions $$ S(2) = 0, $$ and similarly for $T(q)$. You can then use any of the common ODE integrators in matlab, mathematica, maple, ..., to solve and plot the solutions so that you can ...


3

The only chance you stand to deal with this problem from a numerical perspective is oscillatory integration methods. Filon/Levin-type methods can sometimes handle problems like this, particularly when they are of $\sin$ or $\cos$ type, though the $\Gamma$ function and its run-away growth may be prohibitive for the $q$ you are hoping for. In any case, I was ...


3

The identity you're missing from Gauss' divergence theorem is: $$ \int_\Omega \nabla \varphi \cdot\mathbf{v} = -\int_\Omega \varphi\nabla\cdot\mathbf{v} +\int_{\partial\Omega}\varphi\mathbf{v\cdot n} $$ where I've written $\varphi$ as an arbitrary scalar field. So, using the divergence of $\mathbf{u}$ as the scalar field you'd get $$ -\int_\Omega(\lambda+\mu)...


3

The general form of the equation is $$ \frac{\partial \sigma_{ij}}{\partial x_j} + F_i = \rho \frac{\partial^2 U_i}{\partial t^2} $$ where the stress is given by $$ \sigma_{ij} = \sigma_{ij}(U) = 2 \mu \varepsilon_{ij} + \lambda \varepsilon_{kk} \delta_{ij}, \qquad \varepsilon_{ij} = \varepsilon_{ij}(U) = \frac{1}{2}\left( \frac{\partial U_i}{\partial x_j} + ...


3

A typical use case is when you are fine with a fixed low number of nodes, for instance $n=4$ or $n=10$. In that case, you hardcode the weights in your code, and you get higher precision than the trapezoidal rule at basically no cost.


3

I would recommend 3 things to pursue: Use ode15s (a stiff solver). This will allow for better resolution near the jumps Rescale the problem. Your coefficients are very large and if they are changing by several orders of magnitude, this is likely causing undue conditioning problems with your solver so rescaling then undoing the scaling after the ode solve is ...


3

Yes, it is very well arguable that $\frac{e^x-1}{x}$ is ``as right as'' expm1. It appears frequently in applications; see for instance Sections 2.1 and 10.7.4 of Higham's Function of Matrices, which mention exponential integrators and matrix functions; this function is called $\psi_1(x)$ there. Your example, integrating exponentials, is another one, and it ...


2

I would either interpolate the samples of $u_0$ or do a curve-fit, so as to obtain an approximate $u_0(\eta)$ that is easy to integrate analytically. The choice of interpolation or fitting methods depends on what you know about $u_0(\eta)$. For example, if it is a band-limited signal, and the sample frequency is greater than the Nyquist rate, then ...


2

expm1 is a primitive because for small $x$, the Taylor series for $e^x=1+x+x^2+\cdots$ adds small terms to one -- which implies that if $x$ is small enough, you lose almost all digits of accuracy other than the leading one and so the expression $e^x-1$ can not be accurately computed by first evaluating the exponential. The fact that you happen to have an ...


2

The phenomenon you're dealing with here is called a "stiff system". Check here for a full explanation: http://www.scholarpedia.org/article/Stiff_systems Matlab’s ode45 is a non-stiff solver, thus ill-suited for your problem. You should try to use one stiff solver. Matlab has 4 stiff solvers but other tools like Scipy’s solve_ivp also possess non-stiff ...


2

Your function func(rho, gamma) takes two arguments while you pass only one (I think none). That is what your error tells you. Recheck where you use func(rho, gamma) and make sure you pass two arguments.


2

The given $$ \int_\Omega v \frac{\partial}{\partial x}\left( \frac{\partial u}{\partial x} + \frac{\partial u}{\partial y}\right) dx dy = 0 $$ becomes $$ \int_\Omega \frac{\partial}{\partial x}\left[v\left( \frac{\partial u}{\partial x} + \frac{\partial u}{\partial y}\right) \right] dx dy - \int_\Omega \frac{\partial v}{\partial x}\left( \frac{\partial u}{\...


2

There's a closed form solution. For $a,b,c,d,x_0 > 0$, $$ \int_{x_0}^\infty \frac{dx}{x\sqrt{(a+bx)(c+dx)}} = \\ \frac{\log(2 \sqrt{b d (a + b/x_0) (c + d/x_0)} + a d + b c + 2 b d/x_0) - \log(2 \sqrt{b d a c} + a d + b c)}{\sqrt{b d}} $$ a <- 1 b <- 644.153 c <- 4.17e-5 d <- 0.145 x <- 1/1095 (log(2*sqrt(b*d*(a+b*x)*(c+d*x)) + a*d + b*c + ...


1

I think that you don't need a change of variable (yourself) for this problem. Quadpack seems to work just fine for it. It uses a Gauss-Kronrod quadrature. I tried it (in Python) and it seems to work. import numpy as np from scipy.integrate import quad fun = lambda x: 1/(x*np.sqrt((x + 644.153)*(4.15e-5*x + 0.145))) inte, err = quad(fun, 1095, np.inf) ...


1

You are approximating a definite integral with cumtrapz it won't give you the same result as the integrated equation unless you add a constant and plot with the given x coordinates: import numpy as np import matplotlib.pyplot as plt import scipy.integrate as it x = np.arange(-10,10, 0.01) # start,stop,step f = x**2 f_int=it.cumtrapz(f,x, initial=0) plt....


1

First one should verify the numerical solution for $a(t)$ against the analytic solution. $ d_t a = 1/a + 1/a^2 = \frac{a+1}{a^2} $ Use $b=a+1$, then the ODE becomes $ \frac{(b-1)^2}{b} db = dt, $ from which we find the solution $ b^2/2 - 2 b - \ln(b) = t $ This relation (defined up to a constant to satisfy the initial conditions) is not invertible but ...


1

The reason you didn't find 2D quadrature is that we haven't implemented it yet. As to your compilation failure, this does the trick: #include <boost/math/quadrature/gauss_kronrod.hpp> #include <iostream> int main(int argc, char *argv[]) { using namespace boost::math::quadrature; auto f1 = [](double t, double s) { return std::exp(-(t*t+...


1

There is not enough information to determine if the criticism of the original paper is merited or not! In general, approximations are useless unless we can obtain a reliable and accurate error estimate or an error bound which is not too large. Our ability to do just that is compromised when the integrand is not sufficiently smooth. The demands on the ...


1

In essence, you are asking whether you have values for a function $h'(r)$ at points $r_i$, you can obtain an approximation of $h(r)$. The answer is of course yes: If you connect the points $(r_i,h'(r_i))$ by a piecewise linear curve, then you can integrate that to obtain a piecewise quadratic approximation of $h(r)$. You can be more accurate if you connect ...


1

I think the easiest way may simply be an outer loop over the outer integral, while interatively increasing the inner one. So you basically do a midpoint rule for the outer integral, and remember the value of your inner integral from the last step. pseudocode: double dx = 1e-12; double inner_int =0.0; double outer_int =0.0; for(double x = 0.0; x<=1.0;x+=...


1

Please have a look at the answers to this question on CrossValidated. There you can find a worked out solution (including R-code) based on Gauss-Hermite quadrature, which is especially designed for this problem (note that it requires a variable transformation) and -in a different answer- a comparsion of MonteCarlo (slow) with the R builtin function integrate ...


Only top voted, non community-wiki answers of a minimum length are eligible