Recent Questions - Computational Science Stack Exchange most recent 30 from scicomp.stackexchange.com 2023-03-25T16:31:15Z https://scicomp.stackexchange.com/feeds https://creativecommons.org/licenses/by-sa/4.0/rdf https://scicomp.stackexchange.com/q/42639 0 How can I convert a C program into a CUDA program? user366312 https://scicomp.stackexchange.com/users/39670 2023-03-24T17:12:35Z 2023-03-24T17:24:51Z <p>Suppose I already have a C program that solves a specific computational problem. I want to convert that into a CUDA program.</p> <p>What steps should I follow to do that?</p> <p>For instance, <strong>can I think as follows?</strong></p> <p><strong>Step 1.</strong> select which part you want to parallelize.</p> <p><strong>Step 2.</strong> determine if a GPU-friendly algorithm could be adapted.</p> <p>E.g., in the case of matrix multiplication, we can use the sum of the outer product method.</p> <p><strong>Step 3.</strong> determine how much memory is needed in the GPU. Allocate that memory to the GPU.</p> <p><strong>Step 4.</strong> determine how many threads, blocks, and grids are needed.</p> <p>E.g., Assuming that the matrix multiplication is implemented using this approach, the number of threads, blocks, and grids required would be:</p> <ol> <li>The number of threads per block can be determined based on the hardware limits of the GPU and the size of the tile being computed. A common choice is to use a block size of 16x16 or 32x32, depending on the GPU architecture.</li> <li>The number of blocks per grid can be determined based on the size of the matrices being multiplied and the block size. Specifically, the number of blocks required in each dimension can be computed as ceil(n / block_size), where n is the size of the matrix in that dimension. For example, if the block size is 16 and the size of the first matrix is m*n, then the number of blocks required in the x dimension would be ceil(n / 16), and the number of blocks required in the y dimension would be ceil(m / 16).</li> <li>The number of grids required would depend on the size of the matrices being multiplied and the number of blocks per grid. Specifically, the number of grids required in each dimension can be computed as ceil(matrix_size / (block_size * num_blocks)), where matrix_size is the size of the matrix in that dimension and num_blocks is the number of blocks in that dimension.</li> </ol> <p><strong>Step 5.</strong> ... ... ...</p> <p>If yes, please give me a guideline.</p> https://scicomp.stackexchange.com/q/42638 0 Does exact diagonalization of a matrix allow for efficient computation of a Lanczos basis? miggle https://scicomp.stackexchange.com/users/45267 2023-03-24T13:08:24Z 2023-03-24T13:45:17Z <p>Suppose that we are given a large, real-symmetric matrix <span class="math-container">$L$</span>, which is simply too large to perform exact diagonalization on numerically. If we want to study its spectrum, one tool we can use is the Lanczos algorithm, which reduces <span class="math-container">$L$</span> to a tridiagonal form in some subspace of the full vector space in which <span class="math-container">$L$</span> is defined. Call the orthonormal basis which spans this subspace the &quot;Lanczos basis&quot;. Usually, this tridiagonal form and the Lanczos basis are just byproducts produced along the way towards computing eigenvalues of <span class="math-container">$L$</span>.</p> <p>I am in the strange position where my goal is to invert this procedure. For my problem, <span class="math-container">$L$</span> is a very large matrix that I happen to have good analytic control over, so I can compute its eigenvalues and eigenvectors by hand. My question is, is it possible to use this knowledge to efficiently construct the Lanczos basis and corresponding tridiagonalization of <span class="math-container">$L$</span>?</p> <p>It's hard to imagine that knowledge of the complete spectrum doesn't give a computational advantage, but so far I have not been able to figure out how to use it. It's an unfortunate position to be in, because the Lanczos algorithm suffers from severe numerical stability issues and I would rather avoid running it; there should be a basis transformation which relates the Lanczos basis of <span class="math-container">$L$</span> to its eigenspace, but I don't know how to construct it. Any ideas from linear algebra practitioners?</p> <p><strong>Edit:</strong> It was asked in the comments (1) if the initial Lanczos vector is fixed or arbitrary and (2) what exactly is needed from the output of my computations. The initial Lanczos vector (call it <span class="math-container">$q_0$</span>) I'll look at will depend on the precise application, but I think it's sufficiently generic to consider the case where <span class="math-container">$q_0$</span> is not an eigenvector of <span class="math-container">$L$</span>.</p> <p>As for the output that I need, what I am ultimately after are the off-diagonal elements of <span class="math-container">$L$</span> when it is put in tridiagonal form. That is, I don't need specific knowledge of the Lanczos basis per se, only the elements of <span class="math-container">$L$</span> in that basis. If that makes things simpler, that's great, I just don't see how we could get one without the other.</p> https://scicomp.stackexchange.com/q/42637 1 How to get the inverse FFt in this Fortran code? David https://scicomp.stackexchange.com/users/45263 2023-03-23T07:12:22Z 2023-03-23T07:36:55Z <p>I find this fft algorithm on the link</p> <p>The code looks simple and easy to implement. But it does not have inverse fast Fourier transformation. A brief search on the internet shows that to get the inverse of it, the data should be divided by the length of the array. Dividing this, however, does not produce the original input data. So which part is missing in getting the inverse of FFT?</p> <pre><code>module fft_mod implicit none integer, parameter :: dp=selected_real_kind(15,300) real(kind=dp), parameter :: pi=3.141592653589793238460_dp contains ! In place Cooley-Tukey FFT recursive subroutine fft(x) complex(kind=dp), dimension(:), intent(inout) :: x complex(kind=dp) :: t integer :: N integer :: i complex(kind=dp), dimension(:), allocatable :: even, odd N=size(x) if(N .le. 1) return allocate(odd((N+1)/2)) allocate(even(N/2)) ! divide odd =x(1:N:2) even=x(2:N:2) ! conquer call fft(odd) call fft(even) ! combine do i=1,N/2 t=exp(cmplx(0.0_dp,-2.0_dp*pi*real(i-1,dp)/real(N,dp),kind=dp))*even(i) x(i) = odd(i) + t x(i+N/2) = odd(i) - t end do deallocate(odd) deallocate(even) end subroutine fft end module fft_mod program test use fft_mod implicit none complex(kind=dp), dimension(8) :: data = (/1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0/) integer :: i call fft(data) do i=1,8 write(*,'(&quot;(&quot;, F20.15, &quot;,&quot;, F20.15, &quot;i )&quot;)') data(i) end do end program test </code></pre> <p>I did write</p> <pre><code>data = data/8 ! 8 is the dimension </code></pre> <p>In matlab, i could just use it like</p> <p><a href="https://i.stack.imgur.com/YdvNJ.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/YdvNJ.png" alt="enter image description here" /></a></p> <p>but working for it in Fortran is so complicated!</p> https://scicomp.stackexchange.com/q/42635 0 Is it normal for a genetic algorithm to not have mating phase? noob234 https://scicomp.stackexchange.com/users/45258 2023-03-22T15:44:07Z 2023-03-22T15:52:11Z <p>I know that is possible for a genetic algorithm to not have a mating phase because some genetic algorithms use a &quot;mutation-only&quot; approach, where new candidate solutions are generated by randomly mutating existing solutions without crossover.</p> <p>But what about creating a genetic algorithm that does not include crossover and mutation? Is that a good approach?</p> <p>I am thinking of creating a recommendation service that retrieves audio files from the database. This service uses an algorithm that calculates the fitness score of each individual based on some properties of that file and after that, I sort the individuals based on the fitness score and get the best candidate.</p> <p>But at this point isn't that just a &quot;scoring&quot; algorithm and not a genetic one? Or that's perfectly fine to call it like that?</p> https://scicomp.stackexchange.com/q/42634 2 Quadrature rules with the weight function $w(x) = |x|^\gamma$ Justin Dong https://scicomp.stackexchange.com/users/6643 2023-03-22T04:51:34Z 2023-03-22T06:29:47Z <p>I am interested in integrals of the form <span class="math-container">$$\int_{[0,1]^{d}} |x|^{\gamma}f(x)dx.$$</span></p> <p><span class="math-container">$\gamma&gt;0$</span> and <span class="math-container">$f$</span> has some singular behavior at <span class="math-container">$\vec{0}$</span>. The weight function <span class="math-container">$|x|^\gamma$</span> is commonly used in weighted Sobolev spaces which include reduced regularity functions.</p> <p>For <span class="math-container">$d=1$</span>, a Gauss-Jacobi quadrature rule with <span class="math-container">$(\alpha,\beta) = (0,\gamma)$</span> can be used to approximate the above integral. For higher dimensions, I don't think this is the case as <span class="math-container">$|x|^\gamma \neq \prod_{i=1}^{d} x_{i}^{\gamma}$</span>. The choice <span class="math-container">$(\alpha,\beta)=(0,0)$</span> (i.e. Gauss-Legendre) does not work well for me since <span class="math-container">$|x|^{\gamma}f(x)$</span> still contains a singularity (depending on the value of <span class="math-container">$\gamma$</span>).</p> <p>I'm wondering if anyone has literature suggestions for integrals of the above form, even just for 2D would be sufficient. I couldn't find anything in my search even though I thought this would be a weight function someone's looked at before.</p> https://scicomp.stackexchange.com/q/42631 0 How to design a sin and an arcsin function such that arcsin(sin(x))=x, where x is a finite precision floating point number user3677630 https://scicomp.stackexchange.com/users/37189 2023-03-21T12:16:56Z 2023-03-21T14:07:16Z <p>As commonly known for programming on computer, if x is a finite precision floating-point number such as double/float in C language, arcsin(sin(x)) is usually <strong>not</strong> equal to x due to the numerical issue. I'm wondering whether there is a way to design a sin function and an arcsin function such that arcsin(sin(x))=x where x is an fixed-size multi-precision floating-point number. Ideally, both sin and arcsin function would have error guarantee.</p> <p>On solution would be to approximate the sin using linear piecewise function. However, it would take a lot of memory to store all functions if we want to achieve a very accurate approximation.</p> https://scicomp.stackexchange.com/q/42629 0 Time and memory required to diagonalize a 18000 by 18000 matrix using numpy in python Snpr_Physics https://scicomp.stackexchange.com/users/44877 2023-03-20T18:28:33Z 2023-03-21T15:36:27Z <p>Can someone give an estimate of the Time and memory required to diagonalize a 20000 by 20000 complex hermitian matrix using numpy in python ?</p> https://scicomp.stackexchange.com/q/42628 2 Order of local error when integrating ODE with discontinous derivatives Tor https://scicomp.stackexchange.com/users/19544 2023-03-20T08:58:26Z 2023-03-20T08:58:26Z <p>I'm working with ODEs, <span class="math-container">$$\dot{x} = f(x, t),$$</span> where the (higher) derivatives of the right-hand side have discontinuities. In particular, <span class="math-container">$f(x, t)$</span> is obtained by interpolation of discrete samples, and the number of continious derivatives thus depend on the order of interpolation.</p> <p>I know from Theorem II.3.1 on page 157 of Hairer, Nørsett &amp; Wanner that if a Runge-Kutta method is of order <span class="math-container">$p$</span>, then the local error is bounded by <span class="math-container">$C h^{p+1}$</span> if all partial derivatives of <span class="math-container">$f(x,t)$</span> up to order <span class="math-container">$p$</span> exist and are continuous.</p> <p>Hence, I know that if I for example use cubic spline interpolation to obtain <span class="math-container">$f(x, t)$</span>, then only the first and second derivatives are continuous, and hence I cannot expect fourth-order accuracy from a fourth-order Runge-Kutta method. In practice, though, I <em>do</em> see fourth-order convergence when <span class="math-container">$f(x, t)$</span> is obtained by cubic spline interpolation, while if <span class="math-container">$f(x, t)$</span> is obtained with linear interpolation (where not even the first derivative is continuous) I get second-order convergence from fourth-order Runge-Kutta.</p> <p>Hairer, Nørsett &amp; Wanner only tells me that I cannot expect a method of order <span class="math-container">$p$</span> to work as advertised unless I have <span class="math-container">$p$</span> continous partial derivatives, but they don't tell me how large the error will be. I've done some numerical experiments, but I'm curious to see if there exists some more rigorous literature on this point.</p> <p>My question is therefore if there exists any (reasonably accessible) theory that I can use to reason about this.</p> https://scicomp.stackexchange.com/q/42626 2 Poisson equation solution in a semiconductor structure Photon https://scicomp.stackexchange.com/users/39675 2023-03-19T14:09:46Z 2023-03-19T20:17:41Z <p>I am trying to solve the <span class="math-container">$\textbf{1-D}$</span> Poisson equation for a semiconductor structure at equilibrium (There is no external bias applied).</p> <p><span class="math-container">$\textbf{Background}$</span></p> <p><span class="math-container">\begin{equation} \frac{d^2V}{dx^2} = -\frac{\rho(V)}{\epsilon}\\ \rho(V) = q(N_D-N_A+p(V)-n(V)) \end{equation}</span></p> <p>When there are two different semiconductors placed in contact, there will be a redistribution of charge, causing a dependence of charge on potential (the bands will bend).</p> <p><span class="math-container">$\textbf{Solving with the Newton-Rhapson method}$</span></p> <p>I am trying to solve the above non-linear equation using the newton-rhapson method. Taking the central difference, the second derivative becomes: <span class="math-container">\begin{equation} \frac{d^2V}{dx^2}|_i = \frac{(V_{i+1}-2V_i+V_{i-1})}{dx^2} \end{equation}</span></p> <p>The issue with solving this now is that this derivative cannot be computed for the two end points of the structure. What are the boundary conditions to be applied on the edges?</p> <p>Another issue I'm facing is when computing the charge densities by including a non-parabolic coefficient <span class="math-container">$\alpha$</span> (for non-parabolic bands), I have to compute these integrals in every iteration - <span class="math-container">\begin{equation} n = \int^{\infty}_{Ec}F_{1/2}(E)g_n(E)dE \hspace{10pt}p= \int^{E_v}_{-\infty}F_{1/2}(E)g_p(E)dE \end{equation}</span> where <span class="math-container">$F_{1/2}$</span> is the 1/2nd order Fermi-distribution and <span class="math-container">$g(E)$</span> is the density of states. I have to compute the integrals with <span class="math-container">$E_c\text{ and }E_v$</span>, but those are unknown. How do I initiate and update their values? How are each of these energies (<span class="math-container">$E_c,E_v,E_f$</span>) related to each other and the potential of the device?</p> https://scicomp.stackexchange.com/q/42625 0 why does the average temperature vary a lot from my setting temperature in nve? ciel https://scicomp.stackexchange.com/users/45241 2023-03-18T21:05:19Z 2023-03-18T21:05:19Z <p>I set my t_start and t_stop as 0.52 (lj unit), and I ran the nvt first.<a href="https://i.stack.imgur.com/ROv0I.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/ROv0I.png" alt="enter image description here" /></a> after it reached equilibration I ran the nve but found the temperature of the last 1000steps went up to 0.56 while my coworker still had the normal result around 0.52. <a href="https://i.stack.imgur.com/WecOd.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/WecOd.png" alt="enter image description here" /></a>we even used the same code but our results are quite different. why does the average temperature vary so much from the one I set before?</p> https://scicomp.stackexchange.com/q/42623 4 Rank-1 correction of matrix exponential Yaroslav Bulatov https://scicomp.stackexchange.com/users/18786 2023-03-18T06:03:03Z 2023-03-18T09:15:42Z <p>I need to approximate the following in <span class="math-container">$O(d)$</span> time for <span class="math-container">$d\times d$</span> diagonal <span class="math-container">$A$</span> and rank-1 <span class="math-container">$B$</span> <span class="math-container">$$u^T \exp(-A+B) v$$</span></p> <p>Here <span class="math-container">$u,v$</span> are vectors in <span class="math-container">$\mathbb{R^+}^d$</span>, <span class="math-container">$A,B$</span> are positive semi-definite and <span class="math-container">$B$</span> is relatively small</p> <p>The following approximations take <span class="math-container">$O(d)$</span> to compute and get me within factor of 2 of true value on sample data</p> <ul> <li>set <span class="math-container">$B$</span> to 0</li> <li>truncate Zassenhaus <a href="https://en.m.wikipedia.org/wiki/Baker%E2%80%93Campbell%E2%80%93Hausdorff_formula#Zassenhaus_formula" rel="nofollow noreferrer">formula</a> at first term</li> </ul> <p>However, adding more terms from the Zassenhaus expansion seems to make the approximation worse, any tips?</p> <p><sup><sub><a href="https://www.wolframcloud.com/obj/yaroslavvb/nn-linear/mathoverflow-DPR1-approx.nb" rel="nofollow noreferrer">Notebook</a></sub></sup></p> https://scicomp.stackexchange.com/q/42620 0 On solving a first order nonlinear differential equation Monster https://scicomp.stackexchange.com/users/45229 2023-03-17T10:04:40Z 2023-03-17T10:57:58Z <p>It all starts with this Cauchy problem:</p> <p><span class="math-container">$$\begin{cases} \sin(2x(t)) -\cos(3x'(t)) = x(t) + x'(t) \\ x(0) = 1 \\ \end{cases} \quad \quad \text{with} \; t \in [0,10]\,.$$</span></p> <p>Not knowing which way to turn, I started Mathematica:</p> <pre><code>Clear[t, x]; xsol = NDSolveValue[{Sin[2 x[t]] - Cos[3 x'[t]] == x[t] + x'[t], x == 1}, x, {t, 0, 10}, Method -&gt; {&quot;EquationSimplification&quot; -&gt; &quot;Residual&quot;}]; Plot[xsol[t], {t, 0, 10}, PlotRange -&gt; All] </code></pre> <p><a href="https://i.postimg.cc/6QFmxQbP/img.png" rel="nofollow noreferrer"><img src="https://i.postimg.cc/6QFmxQbP/img.png" alt="enter image description here" /></a></p> <p>So, seeing no hope of being able to apply Euler or Runge-Kutta-4, I thought of finite differences:</p> <pre><code>t = x = ConstantArray[0, 101]; {dt, t[], x[]} = {0.1, 0, 1}; Do[t[[k]] = t[[k - 1]] + dt; res = Sin[2 y] - Cos[3 (y - x[[k - 1]])/dt] - y - (y - x[[k - 1]])/dt; x[[k]] = NSolve[res == 0, y, Reals][[1, 1, 2]], {k, 2, 101}]; ListLinePlot[Transpose[{t, x}], PlotRange -&gt; All] </code></pre> <p><a href="https://i.postimg.cc/02XJcKB7/img.png" rel="nofollow noreferrer"><img src="https://i.postimg.cc/02XJcKB7/img.png" alt="enter image description here" /></a></p> <p>which works great but… can you do better? Is it really necessary to use <code>NSolve[]</code>?</p> <p>In environments without a similar solver how could this differential equation be solved?</p> <hr /> <p><strong>@Lutz Lehmann</strong>: I think I have to eat some pasta before I get such good ideas!</p> <p>Thanks a lot, I've already tried with Euler and everything works wonderfully.</p> https://scicomp.stackexchange.com/q/42618 1 Solving 2D Poisson Eq with mixed BC's in Python user82261 https://scicomp.stackexchange.com/users/45219 2023-03-16T00:45:40Z 2023-03-16T00:45:40Z <p>I am trying to numerically solve the Poisson's equation</p> <p><span class="math-container">$$u_{xx} + u_{yy} = - \cos(x) \quad \text{if} - \pi/2 \leq x \leq \pi/2 \quad \text{0 otherwise}$$</span></p> <p>The domain is the rectangle with vertices <span class="math-container">$(-π, 0), (-π,2), (π,0)$</span> and <span class="math-container">$(π,2)$</span>. The boundary conditions are mixed:</p> <p>Dirichlet: <span class="math-container">$u(x,0) = 0$</span> Neumann: <span class="math-container">$u_y(x,2) = 0$</span> Periodic: <span class="math-container">$u(-\pi,0) = u(\pi,0) = 0$</span>.</p> <p>I have the code written, but I can't figure out why my code isn't plotting the correct solution.</p> <p>My solution:</p> <p><a href="https://imgur.com/a/fnDEvvl" rel="nofollow noreferrer">https://imgur.com/a/fnDEvvl</a></p> <p>Correct solution:</p> <p><a href="https://capture.dropbox.com/4hzsLwatNsc05Fhl" rel="nofollow noreferrer">https://capture.dropbox.com/4hzsLwatNsc05Fhl</a></p> <p>Here is my code:</p> <pre><code>import numpy as np from scipy.sparse import diags, csc_matrix from scipy.sparse.linalg import spsolve import matplotlib.pyplot as plt import seaborn as sns import math x_min = -np.pi #Left endpoint of x interval x_max = np.pi #Right endpoint of x interval y_min = 0 #Left endpoint of y interval y_max = 2 #Right endpoint of y interval nx = 50 # Number of grid points in x ny = 50 # Number of grid points in y dx = (x_max - x_min)/nx # Spacing in x direction dy = (y_max - y_min)/ny # Spacing in y direction n = nx*ny # Dimension of system x = np.linspace(x_min, x_max-dx, nx) # Grid points in x-direction y = np.linspace(y_min+dy, y_max, ny) # Grid points in y-direction xg, yg = np.meshgrid(x, y) # Create the diagonals ones = np.full(n, 1) a=np.ones(ny-1,dtype='float') b=np.array([0.0],dtype='float') c=np.concatenate((a,b),axis=0) upper_diagonal = np.tile(c,nx) d=np.ones(ny-2,dtype='float') e=np.array([2.0,0.0],dtype='float') f=np.concatenate((d,e),axis=0) lower_diagonal = np.tile(f,nx) # Create the offsets offsets = [0, -1, 1, ny, -ny, (nx-1)*ny, -(nx-1)*ny] # Create the sparse matrix A = diags([-2*ones/dx**2 -2*ones/dy**2,lower_diagonal/dx**2,upper_diagonal/dx**2,ones/dy**2,ones/dy**2,ones/dy**2,ones/dy**2], offsets, shape=(n, n), dtype=float) # Construct the RHS b = np.zeros((ny,nx),dtype=float) for j in range(nx): #for i in range(ny): if abs(x[j]) &lt;= np.pi/2: b[:,j] = -np.cos( x[j] ) b1 = b.reshape(n,1) #Solve the linear system using a sparse matrix solver As = csc_matrix(A) bs = csc_matrix(b1) u = spsolve(As, bs).reshape(ny, nx) #Plot the solution fig = plt.figure(figsize=(8, 8)) ax = plt.axes(projection='3d') surf = ax.plot_surface(xg, yg, u) ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('u'); plt.show() </code></pre> <p>The only mistake I can think of is in the construction of the matrix <span class="math-container">$A$</span>. Because I have been able to write a similar code where I can solve a number of BVPs with all Dirichlet conditions. So I must not be taking into account the Neumann and periodic boundary conditions correctly into account. However, I think the matrices I generate are correct. For instance, here's the heatmap of the matrix when <span class="math-container">$n_y = 7, n_x = 4$</span>:</p> <p><a href="https://imgur.com/a/ep4Qun8" rel="nofollow noreferrer">https://imgur.com/a/ep4Qun8</a></p> <p>For simplicity, I have shown the matrix prior to having divided each entry by <span class="math-container">$\delta x, \delta y$</span>. The distribution of <span class="math-container">$-4, 1, 2$</span>'s seems right to me, so I don't know what's going wrong. Suggestions?</p> https://scicomp.stackexchange.com/q/42617 2 Schur complement formulation of linear system Leonardo https://scicomp.stackexchange.com/users/43096 2023-03-15T22:42:19Z 2023-03-16T18:11:58Z <p>Consider a system of the following form:</p> <p><span class="math-container">$$(A+K)x=b$$</span></p> <p>where <span class="math-container">$A$</span> is symmetric, positive definite and block diagonal (in fact, a block diagonal matrix made of stiffness matrices arising from FEM discretizations at multiple timesteps, grouped together). Multiplying our linear system by <span class="math-container">$P^t=K^tA^{-1}+I$</span>, we get <span class="math-container">$$(K+K^t + K^tA^{-1}K+A)x = P^tb$$</span></p> <p>We assume that <span class="math-container">$K$</span> is non-symmetric, has full rank and that <span class="math-container">$K+K^t$</span> is positive definite, so that this is a symmetric and positive definite system of equations. The only inconvenient here is the presence of <span class="math-container">$A^{-1}$</span>, and the fact that one needs to solve a (block diagonal) system to obtain <span class="math-container">$P^tb$</span>.</p> <p>Yet another way is to further rewrite everything recognizing a Schur complement:</p> <p><span class="math-container">$$\left\lbrack \begin{array}{cc} A&amp; -K^t\\ -K^t&amp; -K-K^t-A \end{array} \right\rbrack \left\lbrack \begin{array}{c} p\\ x\end{array} \right\rbrack= \left\lbrack \begin{array}{c}-b\\ -b\end{array}\right\rbrack$$</span></p> <p><em>Which one of the three systems would you rather work with, and with which solver would you work? I am interested in wall-clock time and parallelizability</em></p> <p>We assume our blocks to be large, e.g. of size <span class="math-container">$1e8-1e10$</span>.</p> <hr /> <p>To connect the first two formulations, an idea might be to solve formulation one with pre-conditioned GMRES, and use <span class="math-container">$P^t$</span> as preconditioner. Applying <span class="math-container">$P^t$</span> can be done with the (preconditioned) conjugate gradient, and this is highly parallelizable since <span class="math-container">$A$</span> has a diagonal block structure.</p> <p>But in this way, since we are solving essentially formulation two, we'd still need to come up with a good preconditioner for it.</p> <p>Edit. Also because formulation two seems much worse conditioned than formulation one, as numerical experiments show.</p> https://scicomp.stackexchange.com/q/42616 2 Solving 2D Poisson Eq with Dirichlet BC's in Python user45217 https://scicomp.stackexchange.com/users/45217 2023-03-15T19:25:39Z 2023-03-15T23:54:37Z <p>Question:</p> <p>I am trying to solve the following PDE:</p> <p><span class="math-container">\begin{align*} u_{xx} + u_{yy} = \begin{cases} - \cos(x) \quad -\pi/2 \leq x \leq \pi/2, \\ 0 \quad \text{otherwise} \end{cases} \end{align*}</span> The domain is the rectangle with vertices <span class="math-container">$(-π, 0), (-π,2), (π,0)$</span> and <span class="math-container">$(π,2)$</span>. The boundary conditions are Dirichlet: <span class="math-container">$u$</span> must be zero on the boundary.</p> <p>Code:</p> <pre><code>import numpy as np from scipy.sparse import diags, csr_matrix from scipy.sparse.linalg import spsolve import matplotlib.pyplot as plt import math x_min = -np.pi #Left endpoint of x interval x_max = np.pi #Right endpoint of x interval y_min = 0 #Left endpoint of y interval y_max = 2 #Right endpoint of y interval #mesh sizes and grid points dx = 0.05 dy = 0.01 nx = int((x_max - x_min)/dx)-1 ny = int((y_max - y_min)/dy)-1 n = nx*ny x = np.linspace(x_min+dx, x_max-dx, nx) # Grid points in x-direction y = np.linspace(y_min+dy, y_max-dy, ny) # Grid points in y-direction xg, yg = np.meshgrid(x, y) # Create the diagonals ones = np.full(n, 1) a=np.ones(ny,dtype='float') b=np.array([0.0],dtype='float') c=np.concatenate((a,b),axis=0) lower_diagonal = np.tile(c,nx) upper_diagonal = np.tile(c,nx) # Create the offsets offsets = [0, -1, 1,ny,-ny] # Create the sparse matrix A = diags([-2*ones/dx**2 + -2*ones/dy**2,lower_diagonal/dy**2,upper_diagonal/dy**2,ones/dx**2,ones/dx**2], offsets, shape=(n, n), dtype=float) # Construct the RHS b = np.zeros((ny,nx)) for i in range(nx): if abs(x[i]) &lt;= np.pi/2: b[:, i] = -np.cos( x[i] ) b1 = b.reshape(n,1) #Solve the linear system using a sparse matrix solver As = csr_matrix(A) bs = csr_matrix(b1) u = spsolve(As, bs).reshape(ny, nx) #Plot the solution fig = plt.figure(figsize=(8, 8)) ax = plt.axes(projection='3d') surf = ax.plot_surface(xg, yg, u) ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('u'); </code></pre> <p>Here is what my solution looks like:</p> <p><a href="https://imgur.com/a/kznhrZ5" rel="nofollow noreferrer">https://imgur.com/a/kznhrZ5</a></p> <p>Here's the correct solution:</p> <p><a href="https://imgur.com/a/oAWq38K" rel="nofollow noreferrer">https://imgur.com/a/oAWq38K</a></p> <p>However, I can't figure out where my code is wrong. The only non-trivial thing I have done is generated the matrix A using scipy. I think my matrix is correct, though. :/ With <code>dx = 1.2 (nx = 4), dy = 0.4 (ny = 4)</code>, my matrix is of the following form:</p> <p><a href="https://imgur.com/a/5elDvhG" rel="nofollow noreferrer">https://imgur.com/a/5elDvhG</a></p> <p>I think this is correct, no?</p> https://scicomp.stackexchange.com/q/42615 0 Odeint - Solve second order differential equation with a list as argument danial https://scicomp.stackexchange.com/users/43454 2023-03-15T16:23:40Z 2023-03-15T19:20:12Z <p>(<strong>please read the comments after the question</strong>) I'm trying to solve a second order differential equation with a list as argument. the equation is: <a href="https://i.stack.imgur.com/W1B9h.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/W1B9h.jpg" alt="diff-eq" /></a></p> <p>before this, I have obtained <code>q7</code> by solving a first order differential equation and <code>len(q7)=len(t)</code>. I need to put elements with same index from <code>t</code> and <code>q7</code> in the <code>odeint</code> to find the evolution of the second order diff-equation with respect to the evolution of <code>q7</code>; but that gives this error. I'm sure there is not such number <code> ValueError: 0.0001220704714498965 is not in list</code> in <code>q7</code> or <code>t</code> and <code>w</code>, so where this number came from? someone please be kind and help me; what is this error and what should I do?</p> <pre><code>import numpy as np from scipy.integrate import odeint import numpy as np # Define initial conditions and time points z0 = [4.0714*(10**(26)),0] t=10**31*np.logspace(0.247237,2.35443,10**6) t=t.tolist() w=t.copy() q7=np.load('q7.npy') q7=q7.tolist() def arg1_func(t): return q7[w.index(t)] def f(z,t,arg1_func): return [z,-2*(1.4441*1e-33)*(np.sqrt((0.31*((arg1_func(t))**(-4))/3400)+0.31*((arg1_func(t))**(-3))+0.69))*z-(((1e-29)*(arg1_func(t)))**2)*z] sol = odeint(f, z0, t, args= (arg1_func,)) </code></pre> <p><em><strong>and the error:</strong></em></p> <pre><code> --------------------------------------------------------------------------- ValueError Traceback (most recent call last) /tmp/ipykernel_61730/3491914321.py in &lt;module&gt; 20 # Solve the differential equation using odeint 21 ---&gt; 22 sol = odeint(f, z0, t, args= (arg1_func,)) /usr/lib/python3/dist-packages/scipy/integrate/_odepack_py.py in odeint(func, y0, t, args, Dfun, col_deriv, full_output, ml, mu, rtol, atol, tcrit, h0, hmax, hmin, ixpr, mxstep, mxhnil, mxordn, mxords, printmessg, tfirst) 239 t = copy(t) 240 y0 = copy(y0) --&gt; 241 output = _odepack.odeint(func, y0, t, args, Dfun, col_deriv, ml, mu, 242 full_output, rtol, atol, tcrit, h0, hmax, hmin, 243 ixpr, mxstep, mxhnil, mxordn, mxords, /tmp/ipykernel_61730/3491914321.py in f(z, t, arg1_func) 9 10 def f(z,t,arg1_func): ---&gt; 11 return [z,-2*(1.4441*1e-33)*(np.sqrt((0.31*((arg1_func(t))**(-4))/3400)+0.31*((arg1_func(t))**(-3))+0.69))*z-((m[j]*(arg1_func(t)))**2)*z] 12 13 /tmp/ipykernel_61730/3491914321.py in arg1_func(t) 4 5 def arg1_func(t): ----&gt; 6 return q7[w.index(t)] 7 8 ValueError: 0.0001220704714498965 is not in list </code></pre> <p><em><strong><strong>update</strong></strong></em></p> <p><strong>here is <code>q7</code></strong>: <a href="https://filebin.net/qedjynink4kr7395" rel="nofollow noreferrer">q7.npy</a></p> https://scicomp.stackexchange.com/q/42613 2 Need help to fully understand SciPy's odeint's reported step sizes, eval times, # of funct calls & total proc. time (re. question in Astronomy SE) uhoh https://scicomp.stackexchange.com/users/17869 2023-03-15T01:42:00Z 2023-03-21T13:15:36Z <p>A recent question in Astronomy SE <a href="https://astronomy.stackexchange.com/q/53137/7982">Numerical Programming using odeint takes more than 17 minutes</a> got me interested in looking closer at <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html" rel="nofollow noreferrer">SciPy's odeint</a>.</p> <p>The problem is a modified orbital mechanical problem in the solar system. I'd used odeint with great success in orbital mechanics without really understanding the details of how it works under the hood. I chose a reasonably small rtol, enough evaluation times for the output to make nice plots and if necessary rescaled from billions of kilometers and seconds to distances and times of order unity in AU and years.</p> <p>I had always assumed that the integrator chose its own step sizes internally, and when it was all finished <em>re-evaluated</em> the results at the requested times by interpolation.</p> <p>odeint kept the integration results and interpolation coefficients hidden from the user, and the newer <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp" rel="nofollow noreferrer">solve_ivp exposes these</a> if you set <code>dense_output=True</code>. This will allow one to choose a specific region of interest and do a fine-grid interpolation at a later time.</p> <p>It also allows for a bigger choice of integration routines.</p> <p>In the script below, I get the same basic answer (1E-08 variation in final positions and velocities) whether I request the output evaluate at 10 time steps or a million.</p> <p>In all cases there's about 500 function evaluations, which makes sense considering the final results are the same.</p> <ol> <li>But the number of nonzero values in <code>info['hu']</code> varies from 9 to 274, and I haven't a clue what that means.</li> <li>Total process time increases from 5 to 750 milliseconds. Is that due to the time spent interpolating being 100x longer than the integration for 1 million output values?</li> </ol> <p>I wrote a short script (included below) for a 1D strongly anharmonic oscillator, here are the results:</p> <p><a href="https://i.stack.imgur.com/u5IJF.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/u5IJF.png" alt="1D anharmonic oscillator with SciPy odeint and different number of output times" /></a></p> <p><a href="https://i.stack.imgur.com/YkAVI.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/YkAVI.png" alt="1D anharmonic oscillator with SciPy odeint and different number of output times" /></a></p> <pre><code>import numpy as np import matplotlib.pyplot as plt from scipy.integrate import solve_ivp, odeint import time def deriv(y, t, k1, k2): x, v = y a = -k1 * x - k2 * x**3 return np.hstack((v, a)) k1, k2 = 0.5, 2 y0 = np.array([1, 0], dtype=float) args = (k1, k2) n_steps = [10**i + 1 for i in range(1, 8)] rtol = 1E-10 atol = None t_evals, results, infos, process_times = [], [], [], [] for n_step in n_steps: t_eval = np.linspace(0, 10, n_step) time_start = time.process_time() result, info = odeint(deriv, y0, t_eval, args=args, full_output=True, rtol=rtol, atol=atol) results.append(result) infos.append(info) process_time = time.process_time() - time_start process_times.append(process_time) t_evals.append(t_eval) print(n_step, process_time * 1E+06) if True: fig, (row1, row2) = plt.subplots(2, 2) for t_eval, result in zip(t_evals, results): x, v = result.T row1.plot(t_eval, x) row2.plot(t_eval, v) row1.set_ylabel('position') row2.set_ylabel('velocity') x_finals, v_finals = list(zip(*[result[-1] for result in results])) row1.plot(x_finals - x_finals[-1]) row2.plot(v_finals - v_finals[-1]) plt.show() tot_func_evals = [info['nfe'][-1] for info in infos] eval_times = [np.array(info['tcur']) for info in infos] d_eval_times = [et[1:] - et[:-1] for et in eval_times] milli_seconds = [1000 * t for t in process_times] non_zero_d_eval_times = [[d for d in d_eval_time if d != 0] for d_eval_time in d_eval_times] if True: fig, ax = plt.subplots(1, 1) things = non_zero_d_eval_times, n_steps, tot_func_evals, milli_seconds for thing, n_step, n_func_eval, ms in zip(*things): non_zeros = len(thing) ax.plot(thing, label=str((n_step, non_zeros, n_func_eval, round(ms, 1)))) ax.set_title('non-zero changes in evaluation times') ax.legend() plt.show() </code></pre> https://scicomp.stackexchange.com/q/42609 0 Numerical integration library interfacing with eigen KeynesCoeFen https://scicomp.stackexchange.com/users/45205 2023-03-14T12:51:32Z 2023-03-17T12:32:01Z <p>I am looking for a numerical integration library like <a href="https://www.google.com/url?sa=t&amp;rct=j&amp;q=&amp;esrc=s&amp;source=web&amp;cd=&amp;cad=rja&amp;uact=8&amp;ved=2ahUKEwixw-TXttv9AhUzoVwKHXlICoIQFnoECCcQAQ&amp;url=https%3A%2F%2Fwww.gnu.org%2Fsoftware%2Fgsl%2Fdoc%2Fhtml%2Fmontecarlo.html&amp;usg=AOvVaw1hwqYmuh7FrfQpeIVU6efr" rel="nofollow noreferrer">this one</a>. The examples look very appealing but I see that all test functions use very barebones C arrays.</p> <p>Do you have any recommendations of libraries for numerical integration over functions of multi-variate arguments that bind easily to Eigen?</p> <p>I suppose in a pinch, one could always use Eigen's .data() methods and element accessors to emulate the C arrays in their examples.</p> <p>I just think it would be nice to interface with Eigen in a slightly more streamlined way.</p> <p>EDIT: To make things more precise consider the following problem:</p> <p>We have a known function <span class="math-container">$f(\lambda,a)$</span> and want to numerically integrate over <span class="math-container">$\lambda$</span>, where <span class="math-container">$\lambda \in R^{m}$</span>, so we are looking at a high-dimensional integral and possibly need to use Monte Carlo integration or high-dimensional quadrature methods.</p> <p>Any leads would be highly appreciated!</p> https://scicomp.stackexchange.com/q/42599 2 Solving 2D Poisson equation with nonhomogeneous boundary conditions (Dirichlet) and a source Omer Paz https://scicomp.stackexchange.com/users/45185 2023-03-12T14:06:20Z 2023-03-17T03:25:14Z <p>I am a physicist who is fairly new to numerical analysis, currently, I am trying to simulate a non-linear paraxial equation, and part of my calculation involves solving a 2D Poisson equation with Dirichlet boundary conditions and a source function.</p> <p>In order to solve the 2D Poisson equation, I tried using the iterative relaxation method, meaning I repeated the process:</p> <p><span class="math-container">$$\varphi_\mathrm{new}(x) = \frac{1}{4}(\varphi(x+h,y) + \varphi(x-h,y) + \varphi(x,y+h) + \varphi(x,y-h)-h^2\cdot I(x,y))$$</span></p> <p>With a uniformly spaced grid and <span class="math-container">$N=2048$</span> points on each axis. <span class="math-container">$I$</span> is a source function given (in my case a Gaussian distribution), and there are non-homogeneous boundary conditions. However, this method converges very slowly at this resolution.</p> <p>Are there any other recommended methods that will be applicable to my case? It will be extremely helpful if they are either relatively simple or are available in Python libraries.</p> https://scicomp.stackexchange.com/q/42597 6 How to optimize an approximated matrix multiplication? Zuba Tupaki https://scicomp.stackexchange.com/users/44961 2023-03-12T07:16:33Z 2023-03-15T18:38:54Z <p>[<strong>UPDATING</strong>]</p> <p><strong>The old one is a simplified version of the current one.</strong> Here is a solution based on the answer proposed by professor Bangerth down below. To describe what I am trying to do, first rewrite the question into <span class="math-container">$$\max_X\|(I-\alpha X)^{-1}X\circ AE\|_F^2$$</span> where <span class="math-container">$X,A$</span> still are <span class="math-container">$n$</span> by <span class="math-container">$n$</span> square matrices, <span class="math-container">$E$</span> is a vector with all the elements are one such that <span class="math-container">$E=(1,\cdots,1)^\top$</span>. According to professor Bangerth's transformation, the objective becomes <span class="math-container">$$\max_B\left\|\frac{1}{\alpha}(B-I)\circ AE\right\|_F^2$$</span> where <span class="math-container">$B=(I-\alpha X)^{-1}$</span>. Now the gradient of objective with respect to <span class="math-container">$B$</span> is <span class="math-container">$$\nabla_B\left(\frac{1}{\alpha^2}\text{Tr}(B-I)\circ(AEE^\top A^\top)\circ(B-I)^\top\right) =\frac{1}{\alpha^2}(B-I)\circ AE(I\circ AE)^\top$$</span> By the stochastic gradient descent (SGD) algorithm, <span class="math-container">$B$</span> follows <span class="math-container">$$B^{t+1} = B^t + \eta \nabla_B$$</span> So initialize <span class="math-container">$B$</span> as <span class="math-container">$B^0 = (I-\alpha X^0)^{-1}$</span> for a guess <span class="math-container">$X^0$</span>, <span class="math-container">$B$</span> can keep iterating till convergence.</p> <p>The constraint of <span class="math-container">$B$</span> is the specification of matrix <span class="math-container">$X$</span>, which is a sparse matrix and each column has and only has one one. To update <span class="math-container">$X$</span> for each iteration, compute the objective, denoted by <span class="math-container">$J$</span>, with the <span class="math-container">$B^t$</span> from SGD, find the location of maxima in each row of <span class="math-container">$$A + \alpha J$$</span> then mark the same location in <span class="math-container">$X$</span> as one, the rest keeps being zero, this called <span class="math-container">$X'$</span>. In Python, this can be done by</p> <pre><code>ini = np.argmax(A + alpha * J, axis = 1) X = np.zeros([n,n]) X[np.arange(n),ini] = 1 </code></pre> <p>The updated <span class="math-container">$B'$</span> can be obtained by <span class="math-container">$X'$</span> such that <span class="math-container">$B'=(I-\alpha X')^{-1}$</span>. We can also put <span class="math-container">$B'$</span> into SGD to get <span class="math-container">$B^{t+1}$</span>, these iterations stop when <span class="math-container">$\|B'-B^t\|&lt;\epsilon$</span>.</p> <p>The whole thing can be implemented theoretically and pratically, I think. But I still have one pivotal question:</p> <p><strong>How to add the constraint to the objective or express it as a panelty function?</strong></p> <hr /> <p>[<strong>OLD</strong>]</p> <p>Suppose the objective I try to maximize is <span class="math-container">$$\max_{X} \|(I - \alpha X)^{-1}XA\|_F$$</span> where <span class="math-container">$X$</span> is the matrix needs to be pinned down, <span class="math-container">$\alpha$</span> is a scalar, and <span class="math-container">$\|\cdot\|_F$</span> is the Frobenius norm. <strong>Note that <span class="math-container">$(I-\alpha X)^{-1}$</span> is invertible only when <span class="math-container">$\alpha\neq1$</span>, so assume <span class="math-container">$\alpha&lt;1$</span></strong>. All matrices are <span class="math-container">$n$</span> by <span class="math-container">$n$</span> square matrices. An extra condition of matrix <span class="math-container">$X$</span> is that the vectors <span class="math-container">$X(:,k)$</span> for <span class="math-container">$k=1,\cdots,n$</span> is a standard basis vector selected by the optimization, and <span class="math-container">$X$</span> may have repeated columns, for example a 3 by 3 <span class="math-container">$X$</span> could be</p> <p><span class="math-container">$$X = \left( \begin{matrix} 1&amp;0&amp;0\\ 0&amp;1&amp;1\\ 0&amp;0&amp;0 \end{matrix} \right)$$</span></p> <p>One possible way to obtain the maxima of the objective is to apply the matrix multiplication via random sampling, where the objective is rewritten as</p> <p><span class="math-container">$$\max_{p_k}\left\| \frac{1}{n}\sum_{k=1}^n\frac{1}{p_k}B(:,k)X(k,:)A\right\|_F$$</span></p> <p>where <span class="math-container">$B$</span> is the inverse and is approximated by Neumann series such that <span class="math-container">$$B = (I - \alpha X)^{-1} = \sum_{i=0}^n (\alpha X)^k$$</span></p> <p><strong>Note that we do not know the form or the specification of <span class="math-container">$X$</span> in advance but choose each column of <span class="math-container">$X$</span> by random sampling with probability <span class="math-container">$p_k$</span></strong>.</p> <p>Now I have 3 questions:</p> <ol> <li><p>How to compute <span class="math-container">$B$</span> when it is a random sampling approximation?</p> </li> <li><p>Is this method too complex to implement numerically?</p> </li> <li><p>Is there any better way to solve this puzzle?</p> </li> </ol> https://scicomp.stackexchange.com/q/42585 4 Suggestions for libraries that can numerically compute geodesics from a given Riemannian metric? Spencer Kraisler https://scicomp.stackexchange.com/users/45162 2023-03-09T00:52:06Z 2023-03-19T22:45:52Z <p>I am dealing with a non-trivial Riemannian metric <span class="math-container">$H$</span> defined on a particular subset of Euclidean space (<span class="math-container">$E \subset \mathbb{R}^n$</span>). I was able to show the Riemannian manifold <span class="math-container">$(E,H)$</span> is geodesically complete. I am interested in finding a geodesic that connects arbitrary points. I think this is accomplished via the shooting method...</p> <p>I was wondering if there were any libraries, preferably in Matlab or Python, dedicated to numerically computing geodesics on manifolds. I rather not re-invent the wheel and instead take advantage of some (hopefully) library with a dedicated community for fixing bugs and adding cool features I haven't thought of.</p> <p>ManOpt is a great library, but the Manifold object assumes that you have the exponential map readily available. In my case, I do not. I only have a closed-form expression of the Riemannian metric.</p> <p>I attempted at doing this myself but ran into many difficulties. First, the resulting geodesic ODE seems to be <a href="https://en.wikipedia.org/wiki/Stiff_equation" rel="nofollow noreferrer">stiff</a>, and hence Ode45 isn't useful. I found that other ODE solvers for stiff equations work nice (to an extent...), but frankly, I'm dealing with stuff I just don't understand. I want a library written by people who <em>do</em> understand these issues. Not to mention, trying to implement the shooting method for solving 2-point boundary problems for the geodesic will be quite difficult.</p> <p>I would appreciate any help with this.</p> https://scicomp.stackexchange.com/q/42573 -2 Numerically solving for eigenvalues of Schrödinger equation in python [closed] Daniel Lima https://scicomp.stackexchange.com/users/35502 2023-03-05T17:56:38Z 2023-03-19T19:13:07Z <p>I've been strugling with this problem for a while. I'm trying to find and to plot the eigenvalues of a Schrödinger-like equation: <span class="math-container">$$(1-x²)u'' - 2xu' + (\lambda - \phi^{2}(1-x²)-\frac{m²}{1-x²})u,$$</span> <span class="math-container">$\lambda = -2m\phi + E$</span> as a function of a magnetic potential <span class="math-container">$\phi$</span> for some values of angular and magnetic quantum numbers: <span class="math-container">$l = (1...7)$</span> and <span class="math-container">$m = (-l...l)$</span>.</p> <p>As an attempt of a solution, I tried to discretize it using the finite difference method, but I'm not really sure if I did it correctly, since I haven't found examples of this method been applied in general cases. So using the finite difference approximations for the first and second order derivatives one has <span class="math-container">$$\frac{1-x²}{h²}\left(u_{i-1}-2u_i+u_{i+1}\right) +x\left(\frac{u_{i+1}-u_{i-1}}{h}\right) + V(x_i)u_i=\lambda u_i,$$</span> where <span class="math-container">$$V(x)=- \phi^{2}(1-x²)-\frac{m²}{1-x²}$$</span> And in matrix form <span class="math-container">$$\begin{bmatrix}(1-x²)\frac{1}{h^2}+V_1 &amp; \frac{x}{h}+\frac{1-x²}{h²} &amp; 0 &amp; 0...\\ \frac{x}{h}-\frac{1-x²}{h²} &amp; (1-x²)\frac{1}{h^2}+V_2 &amp; \frac{x}{h}+\frac{1-x²}{h²} &amp; 0... \\ ...&amp; ... &amp; ... &amp; \frac{x}{h}+\frac{1-x²}{h²}\\...0 &amp; 0 &amp; \frac{x}{h}-\frac{1-x²}{h²} &amp; (1-x²)\frac{1}{h^2}+V_{N-1} \\ \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ ... \\ u_{N-1} \end{bmatrix} = \lambda \begin{bmatrix} u_1 \\ u_2 \\ ... \\ u_{N-1} \end{bmatrix}$$</span></p> <p>Based on some examples on line, I tried to implement this method in python like this:</p> <pre class="lang-py prettyprint-override"><code> import numpy as np from scipy import sparse as sparse from scipy.sparse.linalg import eigs from matplotlib import pyplot as plt N = 10000 x = np.linspace(-0.99,0.99, N) dx = np.diff(x) phi = np.linspace(0,20,N) m=0 def V(phi,m,x): potential=(phi**2)*(1-x**2)-m**2/(1-x**2) potential_matrix=sparse.diags((potential)) return potential_matrix def kin(phi, m, x): main_diag = (1-x**2)*2*np.ones(N)/dx**2 + V(phi,m)*np.ones(N) off_diag = -x*np.ones(N-1)/dx**2-(1-x**2)/dx**2*np.ones(N-1) kinetic_term = sparse.diags([main_diag, off_diag],(0,1,-1)) return kinetic_term def matrix(phi,m,x): H = kin(phi, m, x) + V(phi,m,x) return H E, v = eigs(matrix,k=100,which='SM') E = np.sort(E) for l in range(1,8): for m in range(-l,l,1): plt.plot(E,phi) plt.grid() plt.xlim([0,20]) plt.ylim([0,60]) plt.tight_layout() plt.show() </code></pre> <p>But I just get an error: AttributeError: 'function' object has no attribute 'shape'</p> <p>So my question is twofold: Did I make the discretization correctly? If not, what is the right way? How can I code this in Python?</p> https://scicomp.stackexchange.com/q/42510 1 What are the benefits of cutting by half the number of multiplications needed to calculate n? user25406 https://scicomp.stackexchange.com/users/45068 2023-02-21T10:39:09Z 2023-03-21T14:36:18Z <p>First we need to present the details of what makes that cut possible. The issue is linked to the sum of squares <span class="math-container">$1^2+2^2+3^2+4^2+...+n^2$</span>.<br /> We consider the two cases:<br /> 1-the case of a sum of even squares<br /> 2-the case of a sum of odd squares.<br /> However both sums lead to practically the same result which is a cut of the number of multiplications by half.</p> <p>Here we consider the two cases, odd and even squares, with a new way of calculating the sum of squares. The sum in both cases involves summing up numbers that are not squares and are located in diagonals perpendicular to the main diagonal of the squares <span class="math-container">$1^2,2^2,3^2,4^2,5^2...$</span>.</p> <p><strong>1- case of the sum of odd squares:</strong> <span class="math-container">$1^2$</span> is just <span class="math-container">$1^2$</span>. However the first sum we consider is <span class="math-container">$1^2+3^2=10$</span> but also <span class="math-container">$1^2+3^2=3+4+3=10$</span>. Here <span class="math-container">$3+4+3$</span> are the elements of the perpendicular to the main diagonal and passing through <span class="math-container">$3,4,3$</span> to add up to <span class="math-container">$10$</span>. The example of <span class="math-container">$1^2+3^2+5^2=35=5+8+9+8+5$</span> is another example. Here also <span class="math-container">$5+8+9+8+5$</span> are the elements of the perpendicular to the main diagonal and going through <span class="math-container">$5,8,9,8,5$</span>. We can give other examples but they all add up to the sum of squares calculated with squares <span class="math-container">$1^2+3^2+5^2+...(2n+1)^2$</span>.</p> <p><strong>2- case of the sum of even squares:</strong> In this case, the first element is always even but all the elements of the diagonal perpendicular to the main diagonal are even as opposed to the odd case where some elements can be even or odd, except of course the first one which is always odd. To give an example, <span class="math-container">$0^2+2^2+4^2=4+16=20=4+6+6+4$</span>. So <span class="math-container">$4,6,6,4$</span> are the elements of the diagonal that add up to <span class="math-container">$20$</span>. Like in the odd case, we can give more examples but they will confirm the rules for sums of odd or even squares. Formulas for the sums are given in the OEIS but no mention is made that they are calculated as shown in this post.</p> <p><strong>Calculating factorial <span class="math-container">$n$</span> using the sum of odd and even squares:</strong></p> <p><strong>1- the odd case</strong><br /> factorial(5) is given by <span class="math-container">$5!=1⋅2⋅3⋅4⋅5=120$</span>. Multiplying the elements given above for <span class="math-container">$1^2+3^2+5^2=35=5+8+9+8+5=35$</span> gives <span class="math-container">$5⋅8⋅9⋅8⋅5=14400=120^2$</span> which is exactly the square of <span class="math-container">$5!$</span>. So it is clear that we need to eliminate some elements of the product <span class="math-container">$5⋅8⋅9⋅8⋅5$</span> and leave only <span class="math-container">$5⋅8⋅3=120$</span> which is basically half the number of multiplications needed to get <span class="math-container">$5!$</span>. This is a general rule that applies to all cases involving factorial of odd numbers. The question of which numbers are involved is shown below. We will just copy the numbers from the multiplication table for factorial of odd numbers to show how the numbers to be multiplied are selected.<br /> <span class="math-container">$3!=3⋅2$</span><br /> <span class="math-container">$5!=5⋅8⋅3$</span><br /> <span class="math-container">$7!=7⋅12⋅15⋅4$</span><br /> <span class="math-container">$9!=9⋅16⋅21⋅24⋅5$</span><br /> These few examples show that the element coming after the starting point <span class="math-container">$5$</span> is: <span class="math-container">$8=5+3$</span>. The element to add to the number whose factorial we want is always the previous odd element in the first column, in this case <span class="math-container">$3$</span>. The next element is <span class="math-container">$3$</span> which is the square root of <span class="math-container">$9$</span>. The number of elements that need to be multiplied is <span class="math-container">$(n+1)/2=(5+1)/2=3$</span>. A reminder for the case of a factorial of an odd number is that each element is always a square away from the element at the intersection of the main diagonal and the perpendicular diagonal. In the case of <span class="math-container">$5$</span> we see that <span class="math-container">$9−5=4=2^2$</span> and <span class="math-container">$8−9=1^2=1$</span>. The last element is always the square root of the square on the diagonal of squares, in this case <span class="math-container">$3$</span>. So now we have the 3 numbers needed to calculate <span class="math-container">$5!$</span>. The last element is simply given by <span class="math-container">$l=(5+1)/2=3$</span>, so in fact, there is no need to take any square root. This has to do with the fact that the diagonal starting with <span class="math-container">$5$</span> is aligned with the square <span class="math-container">$3^2=9$</span> so both <span class="math-container">$5$</span> and <span class="math-container">$9$</span> have the same sum of factors <span class="math-container">$6$</span>. The last element is just the factor <span class="math-container">$q$</span> of <span class="math-container">$q^2$</span>, which in this case is <span class="math-container">$3$</span> or <span class="math-container">$(5+1)/2$</span>.</p> <p>The element coming after the starting point 7 is: <span class="math-container">$7+5=12$</span>. The element to add to the number whose factorial we want is always the previous odd element, in this case <span class="math-container">$5$</span>. This can be easily checked by looking at the multiplication table. The next number to add to <span class="math-container">$12$</span> to get the next element is <span class="math-container">$3$</span> so that <span class="math-container">$7$</span> is distant from the square <span class="math-container">$16$</span> by <span class="math-container">$5+3+1=9$</span>. The number of elements to multiply is <span class="math-container">$(7+1)/2=4$</span> so we have <span class="math-container">$7!=7⋅12⋅15⋅4=5040$</span>. The element before the last number is always <span class="math-container">$1$</span> less than the square of the last element, here <span class="math-container">$15=4^2−1$</span>.</p> <p><strong>2 the even case</strong> The even case is a bit different since the numbers involved are on a line perpendicular to the main diagonal but passing between two squares of the main diagonal. A quick look at the multiplication table will confirm that. In this case there are as many elements on each side of the main diagonal so there is no need to consider the square root or the factors of a square. We give few examples to show how the numbers involved in the multiplication are determined.<br /> <span class="math-container">$2!=1⋅2=2$</span><br /> <span class="math-container">$4!=4⋅6=24$</span><br /> <span class="math-container">$6!=6⋅10⋅12=720$</span><br /> <span class="math-container">$8!=8⋅14⋅18⋅20=40320$</span><br /> Here too we see that the number of multiplications to get <span class="math-container">$n!$</span> has been halved. The same rule to determine the next element applies here. if we look at <span class="math-container">$6!$</span>, we see that we need to add <span class="math-container">$4$</span> to get <span class="math-container">$10$</span> and <span class="math-container">$4$</span> is just the previous even number below <span class="math-container">$6$</span>. The next number is <span class="math-container">$10+2=12$</span> and this will be the last number needed to calculate <span class="math-container">$6!$</span>. The rule is to decrease the number by <span class="math-container">$2$</span> until we get to the last number. We stop when we have added <span class="math-container">$2$</span> to whatever previous number is. In summary, for the factorial of odd numbers <span class="math-container">$1,3,5,...2n+1$</span>, we add odd numbers to the starting number whose factorial we want. For factorial of even numbers, we add even numbers <span class="math-container">$2,4,6...2n$</span>. The numbers are given in reverse order as explained previously.</p> <p>One last thing to show the equivalence between the two methods of calculating a factorial.<br /> <span class="math-container">$8!=1\cdot2\cdot3\cdot4\cdot5\cdot6\cdot7\cdot8=40320$</span>.</p> <p>The above method is equivalent to eliminating some of the factors, two at a time, but without having to do the multiplications. <span class="math-container">$7\cdot2$</span> was replaced by its value with the addition <span class="math-container">$8+6$</span> as explained above. <span class="math-container">$6\cdot3$</span> was replaced by <span class="math-container">$14+4$</span>. Latsly, <span class="math-container">$5\cdot4$</span> was replaced by <span class="math-container">$18+2$</span> so we end up with <span class="math-container">$8!=8\cdot14\cdot18\cdot20=40320$</span>. For the odd case the same principle of replacing factors in the classical way by their values calculated as shown above. The only difference is the handling of the middle term which is always a square and whose value is given by <span class="math-container">$(n+1)/2$</span> with <span class="math-container">$n$</span> the number whose factorial is to be calculated.</p> <p><strong>Edit March 21 2023</strong> A further reduction of the number of multiply for odd numbers.</p> <p>It turned out that we can calculate an odd factorial using numbers from an even factorial. This trick will decrease the number of multiplications needed by one.</p> <p><span class="math-container">$7!=1\cdot2\cdot3\cdot4\cdot5\cdot6\cdot7=5040$</span>.<br /> <span class="math-container">$7!=7\cdot12\cdot15\cdot4=5040$</span></p> <p>Instead of using <span class="math-container">$7!$</span> and getting the numbers needed <span class="math-container">$7,12,15,4$</span>, we consider the even number above <span class="math-container">$7$</span> and get the numbers needed as explained above (for even numbers).</p> <p><span class="math-container">$8!=8\cdot14\cdot18\cdot20$</span>. Then to get <span class="math-container">$7!$</span>, we just drop the first number from <span class="math-container">$8!$</span> to get:</p> <p><span class="math-container">$7!=14\cdot18\cdot20=5040$</span>.</p> <p>So we went from <span class="math-container">$5$</span> (the classical method) multiply to <span class="math-container">$3$</span> (the method above) to <span class="math-container">$2$</span> (with using the numbers from the even factorial) which is more than half the multiplications needed.</p> <p>Another potential improvment would be the use of pre-stored squares and convert the factorial into a product of squares. For example, <span class="math-container">$7!=14\cdot18\cdot20$</span> becomes <span class="math-container">$7!=(16^2-2^2)\cdot20=14\cdot(19^2-1)$</span>.<br /> There are rules that provide the square to be subtracted which depend on the position of the numbers used in the factorial and their distance from the main diagonal of the multiplication table. They are given by <span class="math-container">$1^2,2^2,3^2...$</span></p> <p><strong>Note:</strong> I am not a computer scientist so I am not familiar with the tags needed for this method. I already posted the above in the math stack exchange but it didn't get a good reception. I thought posting it here will get more interest since faster methods of doing calculations are always welcomed. Feel free to add tags, make corrections (if needed).</p> https://scicomp.stackexchange.com/q/40286 3 Calculations on discontinous grids Maxim Umansky https://scicomp.stackexchange.com/users/4325 2021-10-28T04:04:05Z 2023-03-24T06:04:55Z <p>Suppose for a grid-based calculation a grid is used such that the grid Jacobian is discontinuous. For example, in 1D, for a domain <span class="math-container">$x \in$</span> [0,1], one half of the domain is covered uniformly by twice as many grid points as the other side of the domain (and this always remains so as the number of grid points is increased). This grid can be illustrated by this plot, showing the x coordinate of the grid point vs. its index normalized to the total number of grid points.</p> <p><a href="https://i.stack.imgur.com/d9o6nm.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/d9o6nm.jpg" alt="enter image description here" /></a></p> <p>The Jacobian for this grid, representing the transformation of the x coordinate to the grid index is clearly discontinuous at x=1/2. What would be the implications of using such grids for finite-difference of finite-volume calculations such as solving ODEs or PDEs? Would it reduce the grid convergence to the first order? Or not necessarily? Can it lead to losing grid convergence altogether? Would a finite-element (or spectral element) method have an advantage on discontinuous grids, compared to finite difference or finite volume based methods?</p> https://scicomp.stackexchange.com/q/40104 3 Doubt regarding GMRES(m) and preconditioned GMRES Manuel Oliveira https://scicomp.stackexchange.com/users/41264 2021-09-23T22:09:22Z 2023-03-19T04:10:41Z <p>I have the two following algorithms for GMRES(m) and left preconditioned GMRES.</p> <p><strong>GMRES(m)</strong></p> <p><a href="https://i.stack.imgur.com/iNU9s.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/iNU9s.png" alt="enter image description here" /></a></p> <p><strong>Left preconditioning</strong></p> <p><a href="https://i.stack.imgur.com/rRi3n.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/rRi3n.png" alt="enter image description here" /></a></p> <p>I would like to know if anyone could explain why steps 10 through 12 are not used in the left preconditioned GMRES. Are they lost, or is it implicit that these need to happen?</p> <p>If not, am I correct to assume that left-preconditioned works in blocks of m-vectors (number of vectors chosen before the restart of the Krylov subspace). So, if I choose 20 vectors to store before the restart, it will assemble all 20 Krylov subspace basis, and only afterward move to step 12?</p> https://scicomp.stackexchange.com/q/37674 0 deal.ii - ParaView "warp by scalar" of my output is not continuous FEGirl https://scicomp.stackexchange.com/users/33120 2021-06-28T15:42:32Z 2023-03-21T01:06:40Z <p>During our finite element course, we've solved the linear elasticity problem in 2D on a square (<code>GridGenerator::hyper_cube</code>) with <span class="math-container">$Q_1$</span> bilinear finite elements in each component. We imposed neumann homogeneous boundary conditions on one face, and homogeneous Dirichlet on the other three faces.</p> <p>As outputs, we chose:</p> <ul> <li>magnitude of the solution <span class="math-container">$u$</span></li> <li><span class="math-container">$u_x$</span> (x-displacement)</li> <li><span class="math-container">$u_y$</span> (y-displacement)</li> </ul> <p>The output of the magnitude of <span class="math-container">$u$</span> is the following: <a href="https://i.stack.imgur.com/Dvypt.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/Dvypt.png" alt="enter image description here" /></a></p> <p>So far so good. Now, I select <span class="math-container">$u_x$</span>, and I'd like to warp it by scalar, as it is a scalar valued function. So first let's see <span class="math-container">$u_x$</span>: <a href="https://i.stack.imgur.com/F8aMe.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/F8aMe.png" alt="enter image description here" /></a></p> <p><strong>Now, I warp this <span class="math-container">$u_x$</span> by scalar</strong>, and the plot is the following:</p> <p><a href="https://i.stack.imgur.com/1Ngnk.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/1Ngnk.png" alt="enter image description here" /></a></p> <p>i.e. it seems that the solution is flat, which is absolutely non-sense. Also, if I increase the scale factor, I got something which to me doesn't make any sense at all:</p> <p><a href="https://i.stack.imgur.com/itZ2d.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/itZ2d.png" alt="enter image description here" /></a></p> <p>Does anyone know if this is normal, or is there something wrong in my finite element solver? If the latter, this would really surprise me</p> https://scicomp.stackexchange.com/q/34024 3 Estimation of viscosity from critical properties kedarb https://scicomp.stackexchange.com/users/33553 2019-12-12T11:12:29Z 2023-03-17T18:41:31Z <p><a href="https://i.stack.imgur.com/AyHD8.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/AyHD8.png" alt="O. A. Uyehara and K. M. Wastson, *Nat. Petroleum News, Tech. Section*, 36, 764(Oct. 4, 1944); revised"></a></p> <p>The above graph represents <strong><code>reduced viscosity</code></strong> as a function of <strong><code>reduced temperature</code></strong> for several values of the <strong><code>reduced pressure</code></strong>. </p> <hr> <p>I am writing a code which will estimate the viscosity, in the following steps :</p> <ol> <li>Calculating <strong><code>critical viscosity</code></strong>(<span class="math-container">$\mu_c$</span>) by using the formula, </li> </ol> <p><span class="math-container">$$\mu_c = 7.70M^{0.5}p_c^{2/3}T_c^{-1/6}$$</span></p> <ol start="2"> <li>Calculating <strong><code>reduced temperature and pressure</code></strong> as,</li> </ol> <p><span class="math-container">$$T_r = \frac{T}{T_c} \\ p_r = \frac{p}{p_c}$$</span></p> <p>Now, from the graph at the top, estimating <span class="math-container">$\mu_r$</span> by using <span class="math-container">$T_r$</span> and <span class="math-container">$p_r$</span> values calculated from the previous step.</p> <ol start="3"> <li>Finally, calculating the predicted value of <span class="math-container">$\mu$</span> as,</li> </ol> <p><span class="math-container">$$\mu = \mu_r\mu_c$$</span></p> <p>(this value of <span class="math-container">$\mu$</span> is unusually a good agreement with the measured value)</p> <hr> <h2>Question: How do I feed/extract the data/equation of the plot (top) which is experimentally generated so that I can also plot/test it in my code?</h2> <p>P.S - All other parameters <span class="math-container">$T_c$</span> , <span class="math-container">$p_c$</span> , <span class="math-container">$T$</span> , <span class="math-container">$p$</span>, <span class="math-container">$M$</span> will be input by the user</p> <hr> <p><strong>REFERED</strong> : O. A. Uyehara and K. M. Wastson, <em>Nat. Petroleum News, Tech. Section</em>, 36, 764(Oct. 4, 1944); revised | Transport Phenomena, 2nd edition, R. Byron Bird, Warren E. Stewart, Edwin N. Lightfoot</p> https://scicomp.stackexchange.com/q/32949 1 Library for Discontinuous Galerkin method: FEniCS vs deal.ii Zxcvasdf https://scicomp.stackexchange.com/users/24703 2019-06-25T05:04:20Z 2023-03-24T14:04:55Z <p>I am aware that both FEniCS and deal.ii are capable of solving problems with Discontinuous Galerkin (DG) method. I would like to specifically know if any of these two softwares can cater these requirements. Other software suggestions are also welcome; I am aware of these two because they are actively developed. I am specifically interested in solving hyperbolic (wave-dominated) problems.</p> <ol> <li>Local DG implementation</li> <li>User defined numerical flux function.</li> <li>Access and modify nodal/modal basis function coefficients. This is required because I want to implement limiters (for shock capturing).</li> <li>Support for both structured and unstructured meshes.</li> </ol> <p>I request FEniCS/deal.ii users to kindly answer.</p> https://scicomp.stackexchange.com/q/32495 0 Specifying mesh spacing for DFT in numpy user30058 https://scicomp.stackexchange.com/users/0 2019-04-24T19:35:46Z 2023-03-20T09:07:10Z <p>I was testing the <a href="https://docs.scipy.org/doc/numpy/reference/routines.fft.html" rel="nofollow noreferrer">.fft package</a> of numpy 1.16.1 in Python 3.7.2. In particular I was trying to verify that the transform resembles the analytical one for: <span class="math-container">$$f(x) = \mathrm{exp}\left[-\left(\frac{x-5}{2}\right)^{2}\right]$$</span></p> <p>I get from <a href="https://www.wolframalpha.com/input/?i=fourier%20transform%20of%20e%5E%7B-(x-2)%5E%7B2%7D%2F4%7D" rel="nofollow noreferrer">Wolfram Alpha</a> that <span class="math-container">$\hat{f} = \mathcal{F}[f]$</span> looks like this:</p> <p><a href="https://i.stack.imgur.com/xfLJU.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/xfLJU.png" alt="FT by Wolfram"></a></p> <p>Then I tried to replicate this plot with numpy and matplotlib, with the following code:</p> <pre><code>import numpy as np import matplotlib.pyplot as plt x = np.arange(0, 10, 1/1000) y = np.exp(-((x-5)**2)/4) y_hat = np.fft.fftshift(np.fft.fft(y)) re_y_hat = np.real(y_hat) im_y_hat = np.imag(y_hat) fig, ax = plt.subplots() ax.plot(x, re_y_hat, "b-", x, im_y_hat, "r-") plt.show() plt.close() </code></pre> <p>But the image I obtain differs greatly from the one Wolfram gives:</p> <p><a href="https://i.stack.imgur.com/Xh3LK.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/Xh3LK.png" alt="DFT in Python"></a></p> <p>In the last image the zero frequecy was shifted to the center by using <code>np.fft.fftshift()</code> so the spike corresponds to frequency zero.</p> <p>I already figured out that the problem is that nowhere in <code>np.fft.fft()</code> is the <span class="math-container">$\Delta x$</span> being specified, so what numpy is interpreting is that I have data that varies very slowly, almost constant<span class="math-container">$^{1}$</span>, and thus the transform is close to that of a constant function.</p> <p>I looked at the numpy documentation and other SE posts to see how this can be fixed but found nothing. <strong>Does anyone know how to fix this?</strong></p> <hr> <p><span class="math-container">$^{1}$</span> We can calculate the average slope of the function numpy <em>sees</em> by <span class="math-container">$\frac{\mathrm{max}\{f\}-\mathrm{min}\{f\}}{x_{f_{\mathrm{max}}}-x_{f_{\mathrm{min}}}} = \frac{f(5)-f(0)}{n\Delta x} \approx \frac{1}{n\Delta x}$</span> where <span class="math-container">$n$</span> is the number of nodes separating the maximum from the minumum. In this case, since numpy takes <span class="math-container">$\Delta x = 1$</span> by default, the slope is about 1/5000=0.0002</p> https://scicomp.stackexchange.com/q/1273 9 Suggestions for numerical integral over Pólya Distribution yep https://scicomp.stackexchange.com/users/960 2012-02-12T06:13:42Z 2023-03-23T03:02:20Z <p>This problem arises from a Bayesian statistical modeling project. In order to compute with my model, I need to perform an integration in which part of the integrand is the "Pólya" or "Dirichlet-Multinomial" Distribution,</p> <p>$$p(n\mid \alpha) = \frac{(N!) \Gamma(K\alpha)}{\Gamma(\alpha)^K \Gamma ( N + K\alpha)} \prod_{i=1}^K \frac{\Gamma(n_i + \alpha)}{ n_i!}$$</p> <p>where $n_i$ and $N = \sum_{i=1}^K n_i$ are integers, $n = \left(n_1, n_2, \dots, n_K\right)$, and $\alpha &gt; 0$. The integral I wish to compute, $\int_0^\infty (\text{other terms})p(n|\alpha) d\alpha$, works well for small $N$, but the quadrature methods I've attempted (in MATLAB) break down as $N$ becomes large. I haven't tried Monte Carlo; an accurate, fast quadrature method would be very nice for my project. </p> <p>Currently, the "best" method when $N$ is large is to compute $\log[p(n|\alpha)]$ over a grid in alpha, normalize, and exponentiate. This is inaccurate (I lose essentially all detail about the distribution except its peaks), but at least produces a number. </p> <p>I would appreciate any advice on improving this computation, or pointers to different algorithms/methods or existing software. </p> <p>EDIT: I should maybe add that that my evaluation of $p(n|\alpha)$, performed by computing $\log p(n|\alpha)$ using some carefully-written code to compute $\log \Gamma(x)$ for large $x$, does not appear to be causing any problems. </p> <p>EDIT 2: Additionally, "large" values would be on the order of $N\sim 10^8$, with the largest $n_i\sim 10^5$, along with many small values of $n_i$. The other terms are numerically well-behaved. As a simplification with roughly the appropriate tail behavior, you could take </p> <p>$(\text{other terms}) = \exp(-\alpha)$</p>