Recent Questions - Computational Science Stack Exchange most recent 30 from scicomp.stackexchange.com 2021-10-20T06:53:09Z https://scicomp.stackexchange.com/feeds https://creativecommons.org/licenses/by-sa/4.0/rdf https://scicomp.stackexchange.com/q/40250 0 2D DFT for lower frequencies only; is there something significantly faster than numpy.fft.fft2 (throwing away high frequencies)? uhoh https://scicomp.stackexchange.com/users/17869 2021-10-20T03:14:20Z 2021-10-20T04:52:30Z <p>I do a lot of 2D discrete FFT in python using <code>np.fft.fftshift(np.fft.fft2(y))</code>, then throw away 90% <em>or more</em> of the array, keeping only the central low-frequency area.</p> <p>I understand that there's <code>.rfft2()</code> for cases where the input is real.</p> <p>I'm wondering if there are significantly faster ways than numpy's fft to do this, especially if I can specify the low frequency range of interest.</p> <p>I have heard of <a href="https://en.wikipedia.org/wiki/FFTW" rel="nofollow noreferrer">FFTW</a> (as <a href="https://pypi.org/project/pyFFTW/" rel="nofollow noreferrer">pyFFTW</a>) and its ability to choose the best algorithm for a given input, but don't know if there's any extra benefit available for a limited frequency ROI.</p> <hr /> <p>Example of full FT and my ROI:</p> <p><a href="https://i.stack.imgur.com/F0sxc.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/F0sxc.png" alt="enter image description here" /></a></p> <p><a href="https://i.stack.imgur.com/6brT2.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/6brT2.png" alt="enter image description here" /></a></p> https://scicomp.stackexchange.com/q/40248 1 When do not use preconditioners for sparse linear system of equations? yemino https://scicomp.stackexchange.com/users/17630 2021-10-19T22:59:49Z 2021-10-19T23:44:08Z <p>I'm implementing a solver of Finite Element Method, and to solve the linear system of equations I'm using gmres from MKL of Intel. Exists the option with and without a preconditioning. In what case it is not necessary the preconditioning?</p> <p>The matrix is sparse and not symmetric.</p> https://scicomp.stackexchange.com/q/40247 2 Regularisation of ill-conditioned matrix-vector problem Joce https://scicomp.stackexchange.com/users/17850 2021-10-19T12:19:55Z 2021-10-19T15:04:13Z <p>I have a linear* problem which arises from an integro-differential system, and writes: <span class="math-container">$$(\mathbf{I}+\lambda \mathbf{A})x = b$$</span> where <span class="math-container">$\mathbf{A}$</span> is a real full matrix, size <span class="math-container">$n\times n$</span>, but is not symmetric and has some zero eigenvalues. I know what these eigenvalues correspond to.</p> <p>When <span class="math-container">$\lambda$</span> increases (but is still much smaller than 1) the condition number of <span class="math-container">$\mathbf{I}+\lambda \mathbf{A}$</span> (let's call it <span class="math-container">$\mathbf{K}$</span>) deteriorates extremely rapidly (to <span class="math-container">$10^{10}$</span> and more). Even with a direct solver (luckily <span class="math-container">$n=100$</span> is enough for me) the results become meaningless.</p> <p>I have no certainty, but I believe this is linked with the zero eigenvalues of <span class="math-container">$\mathbf{A}$</span>. Is there a way I could exploit my knowledge of the corresponding eigenvectors to regularize the numerical resolution?</p> <p>(* note that this is actually the fixed-point linear problem that I am solving, the original problem has <span class="math-container">$\lambda=\lambda(x)$</span>)</p> <p><strong>EDIT:</strong> <span class="math-container">$\mathbf{A}$</span> itself has indeed a very large condition number (numerically evaluated by numpy with L2-norm), <span class="math-container">$10^{19}$</span>, with <span class="math-container">$\sigma_{\min} = 1.8e-07$</span> and <span class="math-container">$\sigma_{\max} = 3.2e11$</span> (<span class="math-container">$\sigma_{\min}$</span> should analytically be 0 of course, since as said above, <span class="math-container">$\mathbf{A}$</span> has a 0 eigenvalue). Its largest eigenvalue has modulus about 600, here is the complete spectrum for <span class="math-container">$n=102$</span></p> <p><a href="https://i.stack.imgur.com/ozEPg.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/ozEPg.png" alt="enter image description here" /></a></p> https://scicomp.stackexchange.com/q/40239 2 Is there a Python version of the ODE tool pplane? An Ignorant Wanderer https://scicomp.stackexchange.com/users/41323 2021-10-18T00:08:04Z 2021-10-19T11:21:56Z <p>This is the same question as <a href="https://mathematica.stackexchange.com/questions/142076/is-there-a-mathematica-version-of-ode-tools-pplane-and-dfield">this one</a>, except for Python instead of Mathematica. Basically, the MATLAB software <a href="https://math.rice.edu/%7Epolking/odesoft/dfpp.html" rel="nofollow noreferrer">pplane</a> is a staple in ODE courses. Is there a Python equivalent?</p> <p>Sample outputs from pplane: <a href="https://i.stack.imgur.com/oSKc5.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/oSKc5.png" alt="enter image description here" /></a></p> <p><a href="https://i.stack.imgur.com/OvIsC.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/OvIsC.png" alt="enter image description here" /></a></p> <p>EDIT: I don't know much about the software since my class just started using it, but I basically need it for plotting nullclines and trajectories, as well as finding fixed points and their stability. <code>pplane</code> can basically do all of this through a nice GUI without having to write any code.</p> https://scicomp.stackexchange.com/q/40236 0 Computing eigenvalues of Schrodinger equation with spin celerion https://scicomp.stackexchange.com/users/41383 2021-10-17T19:02:35Z 2021-10-19T08:16:48Z <p>I want to solve a 2-dimensional particle in box problem with two electrons in the quantum well.I would like to take into account spin of electrons and Coulomb interactions to compute singlet and triplet eigenstates.</p> <p>My question is it possible to solve this problem using finite difference method, if yes then please guide me through it.</p> https://scicomp.stackexchange.com/q/40235 -1 Null Christoffel symbols associated to the FLRW metric obtained via Mathematica [migrated] Frédéric Laurent https://scicomp.stackexchange.com/users/41544 2021-10-17T18:55:54Z 2021-10-18T17:12:01Z <p>I'm trying to make a mathematica notebook that computes the Christoffel symbols associated to the Friedmann-Lemaître-Robertson-Walker metric (noted FLRW in the notebook) describing an homogenous and isotrope universe that is expanding (or contracting).</p> <p>My issue is that my code fails to recover the Christoffel symbols associated to this metric as they all vanish, which is of course impossible as the described spacetime is curved.</p> <p>Here is presented my code with the various outputs :</p> <p><a href="https://i.stack.imgur.com/7rYOI.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/7rYOI.png" alt="first part" /></a></p> <p><a href="https://i.stack.imgur.com/HE767.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/HE767.png" alt="second part" /></a></p> <p>If my explanations are not clear, do not hesitate to let me know. Thank you in advance for your help ! ^^</p> https://scicomp.stackexchange.com/q/40234 3 What determines the order of a finite volume scheme? CuteCompute https://scicomp.stackexchange.com/users/38580 2021-10-17T13:09:31Z 2021-10-17T13:09:31Z <p>I often hear that cell centred finite volume is second order accurate but at the same time I come across notions of high order FVM flux schemes. Is there a distinction between the two? If I were to use something like a first order upwind referencing vs something like a third order MUSCL would that make my entire FVM scheme first or third order accurate?</p> https://scicomp.stackexchange.com/q/40233 7 What can a computational scientist do in the fourth industrial revolution? Dude https://scicomp.stackexchange.com/users/33505 2021-10-17T10:28:33Z 2021-10-19T13:41:44Z <p>This question is neither scientific nor technical but more career related.</p> <p>I am at a junction in my professional life where I need to make a decision with regard to the future of my career. At the moment I am an R&amp;D manager and my job is very management-oriented (supervising, problem identification) so I asked my own manager about changing my functionality to go on a more technically oriented quest.</p> <p>Following my request, he asked me a simple question that made me understand how unaware I am about the computational science market. He told me: &quot;OK, no problem, but tell me precisely, what is it that you wanna do?&quot;. It's a simple question, right? But it's so difficult for me to find an answer, not because I don't know what I would like to do, but because in my mind I have an image of what I want to do and I can't put it in words. I'll tell you why.</p> <p>I have a master's degree in computational mechanics and a PhD in computational fracture mechanics. What I want to do must include modeling and simulation and programming at the same time. But I can't describe it this way to my manager, I have to be specific about the job in industrial jargon. I'll give you an example. If I tell him: &quot;Hey, listen, I have a background in CFD (Computational Fluid Dynamics)&quot; he would understand and he will easily find me a mission in that area.</p> <p>I am aware of the &quot;fourth industrial revolution&quot; and new automation such as &quot;automation in engineering simulation&quot;, so I am looking for careers in line with future demands. Is automation in engineering simulation a thing? What are the jobs for computational mechanics engineers in the context of industry 4.0?</p> <p><strong>P.S.</strong></p> <p>For those who might say: &quot;Hey dude listen, what you are looking for is software development&quot; I would reply that's right, but, because of my computational background, most software development teams remain skeptical about my software development skills, although I am a skilled programmer.</p> https://scicomp.stackexchange.com/q/40230 0 Help with debugging block GMRES anantdevi https://scicomp.stackexchange.com/users/41476 2021-10-16T19:07:38Z 2021-10-16T19:07:38Z <p>I have written block version of GMRES by referring  and MATLAB implementation of gmres. I need to write it for complex matrices. My block implementation when run on single RHS is giving correct answer but when run with block RHS is converging but LHS is not matching. Can anyone please help me to debug or point in direction how can this be debugged? Are there any standard set of matrix and vectors where intermediate vectors like Krylov subspace basis or H matrix are know (or can be easily computed)?</p> <p> Iterative Methods for Sparse Linear Systems, Saad Yousef</p> <p>Pasting here my blockGMRES and wrapper testcode I am using.</p> <pre><code> function [x, error, iter, flag] = myblockgmres1( A, x, b, max_it, tol ) % input A complex nonsymmetric positive definite matrix % x complex initial guess vector block % b complex right hand side vector block % max_it INTEGER maximum number of iterations % tol complex error tolerance % % output x complex solution vector block % error real error norm % iter INTEGER number of iterations performed % flag INTEGER: 0 = solution found to tolerance % 1 = no convergence given max_it iter = 0; flag = 0; numRHS = min(size(b)); bnrm2 = norm( b ); if ( bnrm2 == 0.0 ), bnrm2 = 1.0; end r = b-A*x ; error = norm( r ) / bnrm2; if ( error &lt; tol ) return, end [n,n] = size(A); m = max_it; p = numRHS; V(1:n, 1:numRHS) = zeros(n,numRHS); H = []; QH = []; cs = []; sn = []; e1 = zeros(n, numRHS); for i=1:numRHS e1(i,i) = 1 + 0i; end r = b-A*x; %calculate QR of initial residue r0 [r0,R] = qr(r, 0); V(:,1:numRHS) = r0(:,1:numRHS); g = e1*R(1:numRHS, 1:numRHS); for j = 1:max_it, Avj = A*V(:,(j-1)*p+1:j*p); for i = 1:j H((i-1)*p+1:i*p, (j-1)*p+1:j*p) = V(:,(i-1)*p+1:i*p)'*Avj; Avj = Avj - V(:,(i-1)*p+1:i*p)*H((i-1)*p+1:i*p, (j-1)*p+1:j*p); end [Q2, R2] = qr(Avj, 0); H(j*p+1:(j+1)*p, (j-1)*p+1:j*p) = R2(1:numRHS, 1:numRHS); V(:, j*p+1:(j+1)*p) = Q2(:, 1:numRHS); %QH is H which gets reduced through Given's QH = [QH; zeros(p, (j-1)*p)]; QH = [QH zeros((j+1)*p, p)]; QH(:, (j-1)*p+1:j*p) = H(:, (j-1)*p+1:j*p); %apply Givens rotation in all previous iterations %to new block of columns for k = (j-1)*p+1:j*p %columns-wise for l = 1:k-1 %reverse order of rows in that block for i = p:-1:1 % apply Givens rotation col = k; row1 = l + i - 1; row2 = l + i; rr = QH(row1, col); hh = QH(row2, col); QH(row1, col) = cs(l, i)*rr - sn(l, i)*hh; QH(row2, col) = conj(sn(l, i))*rr + cs(l, i)*hh; end end %new column reduction for i = p:-1:1 col = k; l = k; row1 = l + i - 1; row2 = l + i; rr = QH(row1, col); hh = QH(row2, col); rnorm = norm(rr); hnorm = norm(hh); rhnorm = sqrt(rnorm*rnorm + hnorm*hnorm); sn(l, i) = -1*(hh*rr) / (rnorm*rhnorm); cs(l, i) = real(rnorm/rhnorm); temp = conj(sn(l, i))*g(row1, :) + cs(l, i)*g(row2, :); g(row1, :) = g(row1, :)*cs(l, i) + sn(l, i)*g(row2, :); g(row2, :) = temp; QH(row1, col) = cs(l, i)*rr - sn(l, i)*hh; QH(row2, col) = 0.0; end end 'iter=', j noexit = 0; %F-norm for i=1:numRHS error(i) = norm(g(j*p+1:(j+1)*p, i)); end norm(error) if(norm(error) &gt; tol) noexit = 1; end %solve triangular systtem if ( noexit == 0 || j == max_it), y = zeros(size(b)); y = QH(1:j*p,1:j*p) \ g(1:j*p,:); %y = y*norm(b); y1(1:j*p, :) = y(1:j*p, :); x = x + V(:,1:j*p)*y1; break; end end if ( error &gt; tol ) flag = 1; end; </code></pre> <p>Testcode to run blockgmres, here you can replace A, RHS and LHS with your test matrices and vectors. If required I can attach test matrix and vector I am using.</p> <pre><code>clear all; A = dlmread('SystemMatrix_1.00e+006_real.txt') + 1i*dlmread('SystemMatrix_1.00e+006_imag.txt'); RHS1 = dlmread('System_1.00e+006_RHS_block1_real.txt') + 1i*dlmread('System_1.00e+006_RHS_block1_imag.txt') ; RHS0 = dlmread('System_1.00e+006_RHS_block0_real.txt') + 1i*dlmread('System_1.00e+006_RHS_block0_imag.txt') ; LHS0 = zeros(size(RHS0)); LHS1 = zeros(size(RHS0)); LHS0_ = dlmread('System_1.00e+006_LHS_block0_real.txt') + 1i*dlmread('System_1.00e+006_LHS_block0_imag.txt') ; LHS1_ = dlmread('System_1.00e+006_LHS_block1_real.txt') + 1i*dlmread('System_1.00e+006_LHS_block1_imag.txt') ; err = 0; iter = 0; %[LHS0,flag,relres,iter,resvec] = gmres(A, RHS0, length(RHS0), 1e-4); %[LHS1,flag,relres,iter,resvec] = gmres(A, RHS1, length(RHS0), 1e-4); %[LHS0, err, iter] = mygmres(A, LHS0, RHS0, length(RHS0), 1e-4); %[LHS1, err, iter] = mygmres(A, LHS1, RHS1, length(RHS0), 1e-4); %call block GMRES LHS = [LHS0 LHS1]; RHS = [RHS0 RHS1]; %[LHS0, err, iter] = myblockgmres1(A, LHS(:,1), RHS(:,1), length(RHS0), 1e-4); %[LHS1, err, iter] = myblockgmres1(A, LHS(:,2), RHS(:,2), length(RHS0), 1e-4); [LHS, err, iter] = myblockgmres1(A, LHS, RHS, length(RHS0), 1e-4); LHS0 = LHS(:,1); LHS1 = LHS(:,2); close all plot(abs(LHS0_)) hold on plot(abs(LHS0)) legend('expected LHS0','LHS0') fprintf('relative error = %f\n', norm(LHS0-LHS0_)/norm(LHS0_)) fprintf('relative error = %f\n', norm(LHS1-LHS1_)/norm(LHS1_)) <span class="math-container">`</span> </code></pre> https://scicomp.stackexchange.com/q/40225 -1 If L is regular so is collapsing doubles [closed] Addem https://scicomp.stackexchange.com/users/34121 2021-10-15T15:49:51Z 2021-10-15T15:49:51Z <p>Suppose L is a regular language. I need to prove that</p> <p><span class="math-container">$$L'=\{x_0\cdots x_n:x_0x_0x_1x_1\cdots x_nx_n\in L\}$$</span></p> <p>I thought I could take a DFA which computes L, and take each accepting state, together with the state before it along an edge labeled 0, and the state before that along an edge labeled 0, and make an edge that skips over the middle state. If you either cut out the middle state or if you leave it and regard the result as an NFA, I can't seem to show that the resulting automaton to computes L'.</p> https://scicomp.stackexchange.com/q/40223 4 Backward Euler + Quasi Newton(Broyden) method fails to solve Van der Pol's equation(Stiff ODE) underdog https://scicomp.stackexchange.com/users/41523 2021-10-15T12:28:48Z 2021-10-15T19:44:18Z <p>The first guess is using the forward Euler approach. The first jacobian is using finite differences. Then NR method is used to solve for the next iteration and Broyden's method is used to update the jacobian until tolerance is met. The figure shows the results that I got compared with correct results. What improvements can I make to my code so that I am moving in the right direction(at the very least). Here is the MATLAB code.</p> <pre><code>function [t,y] = bdf1_qn_broyden(fcn,tspan,Y0,h,tol) t = (tspan(1):h:tspan(2))'; %timestep vector n = length(t); %no of steps y = zeros(n,length(Y0)); %creting output array y(1,:) = Y0; %first row of o/p is initial condition i = 2; jac = zeros(length(Y0)); unity = eye(length(Y0)); while i &lt;= n H = h*unity; for j = 1:1:length(Y0) jac(:,j) = fcn(t(i),y(i-1,:)' + H(:,j)) - fcn(t(i),y(i-1,:)'); end iterate = (unity - jac)\unity; jac = (1/h)*jac; ycurr = y(i-1,:)' + h*fcn(t(i-1),y(i-1,:)'); count = 0; diff = 1; while any(diff &gt; tol) &amp;&amp; (count &lt; 25) curreval = fcn(t(i),ycurr); ynext = ycurr - iterate*(ycurr - y(i-1,:)' - h*curreval); delta = ynext - ycurr; deltaeval = fcn(t(i),ynext) - curreval; jac = jac + (1/dot(delta,delta))*(deltaeval - jac*delta)*delta'; iterate = (unity - h*jac)\unity; diff = abs(delta); ycurr = ynext; count = count +1; end if count &gt;= 25 disp('Iterative method failed to converge within 25 steps'); end y(i,:) = ynext'; i = i+1; end end </code></pre> <p><a href="https://i.stack.imgur.com/exyoq.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/exyoq.jpg" alt="Results" /></a></p> <p>As you might've guessed, I am new to this and my code is quite inefficient. However, I would like to atleast be accurate and slow if I'm going to be slow. Is this not possible without use of adaptive step size? Is BDF1 inept at solving such problems? Should I choose a different method? Or should I solve my implicit function better? Is there a mistake in my implementation? I don't know which way to proceed. Please point me in the right direction. Huge Thanks!</p> https://scicomp.stackexchange.com/q/40222 0 State change of input-output system outlaw https://scicomp.stackexchange.com/users/24800 2021-10-15T08:50:46Z 2021-10-16T08:20:14Z <p><em><strong>Edited</strong></em></p> <p>Given a computer model <span class="math-container">$F:\mathbb{R}^3 \to \mathbb{R}$</span>, with inputs <span class="math-container">$x, w$</span> and <span class="math-container">$z$</span>, and output <span class="math-container">$y=F(x,w,z)$</span>, where for any input, we are able to evaluate the output, my goal is to tune the inputs to achieve a certain state change (change in output).</p> <p>More precisely:</p> <p><strong>Goal</strong>: Given an input <span class="math-container">$x_0, w_0, z_0$</span> and an output <span class="math-container">$y_0$</span>, I want to increase my output to a state <span class="math-container">$y_1$</span>.</p> <p><strong>Question</strong>: How to change my inputs to achieve the new state <span class="math-container">$y_1$</span>.</p> <p><strong>My approach</strong>: <span class="math-container">\begin{align} dy &amp;= \frac{\partial F}{\partial x}(x_0, w_0, z_0)dx + \frac{\partial F}{\partial w}(x_0, w_0, z_0) dw + \frac{\partial F}{\partial z}(x_0, w_0, z_0)dz \\ y_1-y_0 &amp;= \frac{\partial F}{\partial x}(x_0, w_0, z_0) (x_1 - x_0) + \frac{\partial F}{\partial w}(x_0, w_0, z_0) (w_1 - w_0)+ \frac{\partial F}{\partial z}(x_0, w_0, z_0) (z_1 - z_0) \end{align}</span> Since <span class="math-container">$x_0, w_0, z_0, y_0, y_1$</span> are given, and assuming we are able to compute from the model <span class="math-container">$F$</span> the numerical approximation of the derivative of <span class="math-container">$F$</span> with respect to each input at the point <span class="math-container">$(x_0, w_0, z_0)$</span>, the new values <span class="math-container">$x_1, w_1, z_1$</span> then should lie on the plane given by the second equation above.</p> <p>I am also aware that I made a linear approximation that is only exact when <span class="math-container">$F$</span> is linear with respect to the inputs.</p> <p>I am going in the right direction in formulation the problem ? If not, any suggestions?</p> https://scicomp.stackexchange.com/q/40219 2 Different sources of error in Finite Element computations Beni Bogosel https://scicomp.stackexchange.com/users/3548 2021-10-14T10:04:19Z 2021-10-15T02:04:50Z <p>Consider the problem <span class="math-container">$-\Delta u = f$</span> in <span class="math-container">$\Omega$</span>, with <span class="math-container">$u=0$</span> on <span class="math-container">$\partial \Omega$</span>. Suppose that <strong><span class="math-container">$\Omega$</span> is a polygon</strong> and that we approximate the solutions to the previous problem using <strong>Lagrange finite elements</strong> on a <strong>triangular mesh</strong>. When solving such a finite element problem there are multiple sources of error:</p> <ol> <li><p>The discrete mesh may not be an <em>exact</em> mesh for the polygon due to rounding errors (typically of the order of the machine <span class="math-container">$\varepsilon$</span>). This gives slight perturbations for the mass and rigidity matrices.</p> </li> <li><p>The discrete solution comes from solving a linear system. This linear system is solved numerically, so there is an additional error here.</p> </li> <li><p>The difference between the analytical solution and the finite element one is usually quantified with a priori estimates depending on the mesh size <span class="math-container">$h$</span>.</p> </li> </ol> <p>Usually only point 3. above is mentioned when dealing with error estimates from finite elements. I guess this is due to the fact that errors coming from points 1 and 2 above (rounding errors, linear systems) are usually orders of magnitude smaller than the discretization error in 3. However, I am interested in knowing more about how one can quantify the errors coming from the first two points and eventual references treating these aspects.</p> <blockquote> <p>Can you indicate references in the literature or the main ideas that can be used in order to better understand errors coming from points 1-2 above?</p> </blockquote> https://scicomp.stackexchange.com/q/40213 1 Why aren't face integrals for an element calculated in FEM but they show up in FVM? CuteCompute https://scicomp.stackexchange.com/users/38580 2021-10-12T23:18:54Z 2021-10-16T20:28:10Z <p>Consider the Laplace problem: <span class="math-container">\begin{align} -\nabla^2 u = f \qquad \text{in } \Omega \\ u = 0 \qquad \text{on } \Gamma \end{align}</span> The weak problem is find <span class="math-container">$u_h \in V \subset H^1$</span> such that <span class="math-container">$\forall v_h \in K \subset H^1_0$</span>, <span class="math-container">$u_h$</span> satisfies:<br /> <span class="math-container">\begin{equation} \int_\Omega \nabla v_h \cdot \nabla u_h \, d\Omega = \int_\Omega f v_h \, d\Omega \end{equation}</span> FVM is said to be a special case of FEM where the space of test functions is <span class="math-container">$1$</span>. In FVM we get to work with the integral <span class="math-container">\begin{equation} \oint_e \nabla u \, d\Omega \end{equation}</span> for each local element. Why does this integral show up in FVM but not in Galerkin FEM? This doesn't feel consistent to me.</p> https://scicomp.stackexchange.com/q/40207 1 Importance Sampling for multidimensional integrals and random numbers from multivariable pdf's pollux33 https://scicomp.stackexchange.com/users/41455 2021-10-12T15:22:15Z 2021-10-13T08:53:11Z <p>I am aiming to get a numerical value for a five-dimensional integral using Monte Carlo Integration. I am getting good results using the Mean Value Method, but I would like to try to use Importance Sampling to get even better results.</p> <p>My main source gives me an example of how to do this method using a one-dimensional integral, but it doesn't even mention how to replicate it for multivariable integrals. My main issue is generating a set of random numbers from a multivariable probability distribution function.</p> <p>If I understand correctly the process of doing it for a single variable is setting up the equation:</p> <p><span class="math-container">$$\int_0^{x}dx^{\prime}\,p(x^\prime) = z$$</span></p> <p>and then solving for <span class="math-container">$x$</span> (<span class="math-container">$z$</span> are the uniformly distributed random numbers generated).</p> <p>How can I replicate this procedure for an arbitrary 5-dimensional pdf <span class="math-container">$p(x_1,x_2,x_3,x_4,x_5)$</span>?</p> <p>I could find no good sources on this, most of the ones that I found dealt specifically with n-dimensional Gaussian distributions.</p> https://scicomp.stackexchange.com/q/40199 1 Fitting line to a staircase function adch99 https://scicomp.stackexchange.com/users/41477 2021-10-10T15:08:48Z 2021-10-16T04:16:52Z <p>I have a staircase/step function <span class="math-container">$n(E)$</span>. I know the points <span class="math-container">$\{E_i\}$</span> at which each &quot;step&quot; occurs and all steps are of constant height 1. I need to fit a line <span class="math-container">$a + bE$</span> to this function and find the least-squares deviation. In particular, I have to calculate the quantity</p> <p><span class="math-container">$$\Delta = \min_{a,b} \int_{E_i}^{E_f} [n(E) - a - bE]^2 dE$$</span></p> <p>In Python/NumPy, I could try and recreate the function <span class="math-container">$n(E)$</span> with <code>np.heaviside</code> and then try to fit a line to it, but that feels inefficient. Is there a better way to fit a line to a staircase function?</p> <p>One approach might be to break up the integral into parts between each &quot;step&quot; and then optimize the resultant expression wrt <span class="math-container">$a,b$</span>. But I'd like to know if there is a cleaner, more efficient way to fit a line to a step function numerically.</p> https://scicomp.stackexchange.com/q/40197 6 General approach to infinite sums Lewkrr https://scicomp.stackexchange.com/users/41465 2021-10-08T18:54:35Z 2021-10-13T05:01:32Z <p>My question is specific to algorithms and models of computation.</p> <p>I would like to write code to evaluate the following expression quickly and accurately:</p> <p><span class="math-container">$$\log \left( \sum_{i=1}^{\infty}{I_{\nu+i}(2\lambda)} \right)$$</span></p> <p>Where <span class="math-container">$I_\nu(x)$</span> is the modified bessel function of the first kind.</p> <p>I have limited experience with numerical summation of infinite series. Is there a general approach/algorithm for this? <a href="https://en.wikipedia.org/wiki/Kahan_summation_algorithm" rel="nofollow noreferrer">Kahan summation</a> algorithm? Some other general purpose approach? Any reference material to help me learn are greatly appreciated.</p> <p>So far my only ideas are to either iteratively sum to N and stop when <span class="math-container">$\log \left( \sum_{i=1}^{N}{I_{\nu+i}(2\lambda)} \right)$</span> meets some convergence criteria (though I am not certain how precise it ought to be) or I might try summing <span class="math-container">$\log \left( \sum_{i=N}^{\infty}{I_{\nu+i}(2\lambda)} \right)$</span> using a recurrence relation <span class="math-container">$I_{\nu + 1}(x) = I_{\nu-1}(x) - \frac{2\nu}{x} I_{\nu}(x)$</span> until the tail value is small enough, and then calculating the partial sum up to <span class="math-container">$N$</span>.</p> <p>For those interested, I am working in R.</p> <hr /> <p>UPDATE: Based on the comments, I have a few followup points.</p> <ol> <li><strong>Convergence:</strong> I am working with the SDF of the absolute value of a Skellam random variable, i.e <span class="math-container">$P(|W|&gt;\nu)$</span>. In probability theory, the SDF is always between 0 and 1, and it can be shown that the SDF is equal to:</li> </ol> <p><span class="math-container">$$P(|W|&gt;\nu) = 2\exp(-2\lambda)\sum_{i=1}^{\infty}{I_{\nu+i}(2\lambda)}, \ \ \lambda &gt; 0, \nu \in \{0,1,2,3, \dots , \infty\}$$</span></p> <p>so</p> <p><span class="math-container">$$0 &lt; P(|W|&gt;\nu) \le 1 \\ 0 &lt; 2\exp(-2\lambda)\sum_{i=1}^{\infty}{I_{\nu+i}(2\lambda)} \le 1 \\ 0 &lt; \sum_{i=1}^{\infty}{I_{\nu+i}(2\lambda)} \le \frac{\exp(2\lambda)}{2}$$</span></p> <p>demonstrating that the infinite sum does indeed converge.</p> <ol start="2"> <li><p><strong>Range for <span class="math-container">$\lambda$</span> and <span class="math-container">$\nu$</span></strong>: As mentioned in part (1), <span class="math-container">$\lambda$</span> can be any real number greater than zero and <span class="math-container">$\nu$</span> (the random variable) can take on any value in the non-negative integers. Since we are dealing with probabilities, we will assume <span class="math-container">$\lambda$</span> is fixed and attempt to compute the infinite sum for any value of <span class="math-container">$\nu$</span>. I would PREFER an approach to calculating the value that is general for any fixed value of <span class="math-container">$\lambda$</span> and multiple <span class="math-container">$\nu$</span>, but for argument's sake, anyone can take <span class="math-container">$\lambda = 2000$</span> and <span class="math-container">$\nu = 1700$</span></p> </li> <li><p><strong>Convergence rate</strong>: I am not a computer scientist (I am a statistician) so I am not sure that I have a <em>target</em> convergence rate, I only hope to find an accurate approximation of the infinite sum.</p> </li> </ol> https://scicomp.stackexchange.com/q/40145 1 How to compute the Eigenvalue and Eigenstates of Quantum well with Effective mass using finite difference method in Python? celerion https://scicomp.stackexchange.com/users/41383 2021-09-29T11:52:20Z 2021-10-18T19:35:50Z <p>I want to compute the eigenvalues and eigenstates of a quantum well with different effective masses of electron in the barrier and in the quantum well. As can be seen : <a href="https://github.com/mholtrop/QMPython/blob/master/Finite%20Well%20Bound%20States.ipynb" rel="nofollow noreferrer">https://github.com/mholtrop/QMPython/blob/master/Finite%20Well%20Bound%20States.ipynb</a> , the mass of the electron is constant, therefore the Hamiltonian can be easily diagonalized and eigenvalues/eigenstates can be evaluated.</p> <p>However, in my situation, the mass of the electron takes different values in the barrier and the well. Could any one suggest how I could obtain the matrix which satisfies the effective mass condition, and hence can be used to compute eigenvalues and states using similar method.</p> <p><span class="math-container">$$a_i\psi_{i-1}+b_i\psi_{i}+c_i\psi_{i+1}=E\psi_{i}$$</span></p> <p><span class="math-container">$$a_{i+1}=c_i=\frac{-\hbar^2}{2m^*(\delta z)^2} \text{ and } b_i=\frac{\hbar^2}{m^*(\delta z)^2}+V_i$$</span></p> <p>In the above equations, it is apparent that the Hamiltonian can be diagonalized as mentioned in the link above.</p> <p><span class="math-container">$$a_{i+1}=c_i=\frac{-\hbar^2}{2m_{i+\frac{1}{2}}^*(\delta z)^2} \text{ and } b_i=\frac{\hbar^2}{2(\delta z)^2}\bigg(\frac{1}{m^*_{i+\frac{1}{2}}}+\frac{1}{m^*_{i-\frac{1}{2}}}\bigg)+V_i$$</span></p> <p>My question is: how can I diagonalize my Hamiltonian with <span class="math-container">$m^*$</span>(z) being different in the barrier and the well using the finite difference.</p> <p>The complete Schrodinger equation can be written as,</p> <p><span class="math-container">$$\frac{-\hbar}{2}\frac{\partial}{\partial z}\frac{1}{m^*(z)}\frac{\partial}{\partial z}\psi(z)+V(z)\psi(z)=E\psi(z)$$</span></p> <p><span class="math-container">$\qquad\qquad\qquad\qquad\qquad\quad$</span><a href="https://i.stack.imgur.com/NIaaB.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/NIaaB.png" alt="Quantum square-well diagram" /></a></p> <p>Here, mass depends on z growth axis of quantum well (in picture its x instead of z) <span class="math-container">$m^*(z)$</span> in case of finite quantum well the value mass in barrier remains same in both regions 1 and 3 and is different in the well region 2.</p> <p>I would like to know how can I use the same approach as mentioned  to solve the Schrodinger equation with a spatially dependent mass.</p> https://scicomp.stackexchange.com/q/40020 1 Two-dimensional ordering issue – alternate sort order ascending/descending to reduce fluctuations - trivial? Rien https://scicomp.stackexchange.com/users/41177 2021-09-06T14:18:03Z 2021-10-17T23:11:14Z <p>I have a solution in search of a problem that some of you could perhaps help me with.</p> <p>Let <span class="math-container">$L$</span> be a list of elements. Each element has two inherent properties/attributes (<span class="math-container">$a$</span>, <span class="math-container">$b$</span>) that can each be expressed as real numbers. The values of the two attributes of each element are completely uncorrelated. The goal is to re-order list <span class="math-container">$L$</span>, such that for successive elements in the list, the differences in values of both attributes <span class="math-container">$a$</span> and <span class="math-container">$b$</span> are kept relatively small – i.e. we want to avoid big jumps in values from one element to the next.</p> <p>Now, we could minimize fluctuations in attribute a by sorting the list by attribute <span class="math-container">$a$</span>. However, since <span class="math-container">$a$</span> and <span class="math-container">$b$</span> are uncorrelated, the values for <span class="math-container">$b$</span> will then fluctuate randomly between minimum and maximum. If we instead sort by <span class="math-container">$b$</span>, then <span class="math-container">$a$</span> is all over the place.</p> <p>Instead, we opt to create a number of groups/bins based on attribute <span class="math-container">$a$</span>, arrange those bins in order, and then order the elements within these bins by attribute <span class="math-container">$b$</span>. This way, fluctuations in attribute <span class="math-container">$a$</span> can be no larger than two times the range of the bins (worst case: when traversing from one bin into the next), and no larger than one time the range of the bins when traversing within a bin. Fluctuations in attribute <span class="math-container">$b$</span> are also relatively small, as within the bins the values are sorted by <span class="math-container">$b$</span>. The result looks something like this (attribute <span class="math-container">$a$</span> in blue, <span class="math-container">$b$</span> in orange):</p> <p><a href="https://i.stack.imgur.com/MTgsA.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/MTgsA.png" alt="enter image description here" /></a></p> <p>A remaining issue then is that, when traversing from one bin into the next, there remains a big fluctuation in attribute <span class="math-container">$b$</span> – from approximately its maximum value to approximately its minimum (the “sawtooth wave” above). To improve on that, we alternate the sorting order between subsequent bins, to create a sort of “triangle wave”. This ensures that fluctuations in attribute b are small everywhere in the list:</p> <p><a href="https://i.stack.imgur.com/94KDS.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/94KDS.png" alt="enter image description here" /></a></p> <p>I would like to ascertain how trivial it is for an engineer to come up with this ordering strategy – particularly the sorting within bins in an alternating fashion – if the objective is to reduce/minimize fluctuations in <span class="math-container">$a$</span> and <span class="math-container">$b$</span>. To that end, I am looking for problems or applications (or perhaps an entire discipline?) where one is likely to encounter this ordering strategy.</p> <p>While I have found an occasional post about someone trying to <em>implement</em> this particular order (e.g. <a href="https://codereview.stackexchange.com/questions/88710/alternating-sort-order-within-a-sorted-group">https://codereview.stackexchange.com/questions/88710/alternating-sort-order-within-a-sorted-group</a>), it has proven surprisingly difficult to find any concrete examples showing <em>why</em> one would want to do this.</p> <p>Does anything come to mind? Ultimately I aim to find (text)books showing this kind of problem and solution. Anything in that direction would be much appreciated.</p> https://scicomp.stackexchange.com/q/36394 3 Time Reversibility of Velocity Verlet Algorithm Shrey https://scicomp.stackexchange.com/users/37593 2020-11-28T23:17:07Z 2021-10-17T20:09:03Z <p>I'm very new to computational Physics and am finding conflicting statements on whether the velocity Verlet algorithm, defined as:</p> <p><span class="math-container">\begin{align} x_{n+1} &amp;= x_n + v_n \Delta t + \frac{1}{2} a_n \Delta t^2 \tag{1}\\ v_{n+1} &amp;= v_n + \frac{1}{2}(a_n + a_{n+1}) \Delta t \tag{2} \end{align}</span></p> <p>is time-reversible.</p> <p>From my understanding, a method is time-reversible if it is invariant under <span class="math-container">$\Delta t \to -\Delta t$</span>. Applying this to the velocity Verlet algorithm, equation 1 becomes:</p> <p><span class="math-container">\begin{align} x_{n-1} &amp;= x_n - v_n \Delta t + \frac{1}{2} a_n \Delta t^2 \\ x_n &amp;= x_{n-1} + v_n \Delta t - \frac{1}{2} a_n \Delta t^2 \\ x_{n+1} &amp;= x_{n} + v_{n+1} \Delta t - \frac{1}{2} a_{n+1} \Delta t^2 \tag{1'} \end{align}</span></p> <p>and equation 2 becomes:</p> <p><span class="math-container">\begin{align} v_{n-1} &amp;= v_n - \frac{1}{2}(a_n + a_{n-1}) \Delta t \\ v_n &amp;= v_{n-1} + \frac{1}{2}(a_n + a_{n-1}) \Delta t \\ v_{n+1} &amp;= v_{n} + \frac{1}{2}(a_{n+1} + a_{n}) \Delta t \tag{2'} \end{align}</span></p> <p>So far, we see equation 2' is equivalent to equation 2 (so the velocity part is time reversible for step <span class="math-container">$n$</span>), but equation 1' has the wrong times for acceleration and velocity, as well as wrong sign for the acceleration part. However, if we substitute in equation 2 (or 2') into equation 1', this recovers back equation 1 again. Therefore, I would have thought velocity Verlet <em>is</em> time reversible.</p> <p>On the contrary, in <em>Basic Concepts in Computational Physics</em> by Stickler and Schachinger, they state on page 108:</p> <blockquote> <p>The Stormer-Verlet algorithm<span class="math-container">$^{\dagger}$</span> is time-reversal symmetric (invariant under the transformation <span class="math-container">$\Delta t \to - \Delta t$</span>), hence reversible. This is a direct consequence of its relation to the symplectic Euler method. [...] <strong>The leap-frog algorithm or the velocity Verlet algorithm [methods] are not time-reversal invariant.</strong></p> </blockquote> <p>They repeatedly state that leapfrog and velocity verlet are not time-reversible throughout the rest of the book. They define leapfrog as:</p> <p><span class="math-container">\begin{align} x_{n+1} &amp;= x_{n} + v_{n+\frac{1}{2}} \Delta t \tag{3} \\ v_{n+\frac{1}{2}} &amp;= v_{n-\frac{1}{2}} + a_n \Delta t \tag{4} \\ v_{\frac{1}{2}} &amp;= v_0 + \frac{1}{2}a_0 \Delta t \tag{5} \end{align}</span></p> <p>Using the same procedure as for velocity Verlet, I found that equations 3 and 4 are invariant under time reversal, but the initialising condition does change to:</p> <p><span class="math-container">$$v_{\frac{1}{2}} = v_0 + \frac{1}{2}a_{\frac{1}{2}} \Delta t \tag{5'}$$</span></p> <p>where 5' differs from 5 since the acceleration at <span class="math-container">$t = \frac{1}{2} \Delta t$</span> is used instead of at <span class="math-container">$t = 0$</span>.</p> <p>Looking at the half time step implementation of velocity Verlet, defined as:</p> <p><span class="math-container">\begin{align} v_{n + \frac{1}{2}} &amp;= v_n + \frac{1}{2} a_n \Delta t \tag{6} \\ x_{n+1} &amp;= x_n + v_{n + \frac{1}{2}} \Delta t \tag{7} \\ v_{n+1} &amp;= v_{n + \frac{1}{2}} + \frac{1}{2} a_{n+1} \Delta t \tag{8} \end{align}</span></p> <p>I have found under time-reversal, these transform to:</p> <p><span class="math-container">\begin{align} v_{n + \frac{1}{2}} &amp;= v_n + \frac{1}{2} a_{n+\frac{1}{2}} \Delta t \tag{6'} \\ x_{n+1} &amp;= x_n + v_{n + \frac{1}{2}} \Delta t \tag{7'} \\ v_{n+1} &amp;= v_{n + \frac{1}{2}} + \frac{1}{2} a_{n+\frac{1}{2}} \Delta t \tag{8'} \end{align}</span></p> <p>So equation 7 is invariant, but the times for the accelerations at equations 6 and 8 are changed. Does this suggest the substitution of 2' into 1' was invalid and velocity Verlet isn't time reversible? Actually, equation 6' is equivalent to equation 8 and equation 8' is equivalent to equation 6 (as seen by adding or subtracting <span class="math-container">$\frac{1}{2}\Delta t$</span> to the times the equations are evaluated at). So does this mean velocity Verlet is still time reversible?</p> <p>Lastly, in <em>Computational Physics</em> by Thijssen, it's stated that:</p> <blockquote> <p>There exist two alternative formulations of the Verlet algorithm<span class="math-container">$^{\dagger}$</span>, which are exactly equivalent to it in exact arithmetic but which are less susceptible to errors resulting from finite numerical precision in the computer than the original version. The first of these [is] the leap-frog form [and the second is] the so-called velocity-Verlet algorithm which is also more stable than the original Verlet form.</p> </blockquote> <p>To summarise, my questions are:</p> <ol> <li>Are Stickler and Schachinger incorrect in stating that the velocity Verlet and leapfrog algorithms are not time reversible?</li> <li>Is the substitution of 2' into 1' valid to show that velocity Verlet is time reversible?</li> <li>Does the fact that the initialising condition of leapfrog is not time reversible mean that the leapfrog method itself is not time reversible? Other resources, including answers on this stack exchange, state leapfrog <em>is</em> time reversible.</li> <li>Since Thijssen states the velocity verlet and leapfrog algorithm are &quot;exactly equivalent&quot; in terms of the arithmetic, does this mean they inherit the time-reversibility (and possibly even the symplectic nature) of the Stormer-Verlet algorithm?</li> </ol> <hr /> <p><span class="math-container">$^{\dagger}$</span>For reference, the Stormer-Verlet algorithm is outlined <a href="https://en.wikipedia.org/wiki/Verlet_integration#Verlet_integration_(without_velocities)" rel="nofollow noreferrer">here</a>. This Wikipedia page also entails how to obtain velocities within this framework.</p> <p><a href="https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.6496" rel="nofollow noreferrer">https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.6496</a> and <a href="https://math.stackexchange.com/questions/1448005/what-does-the-time-reversibility-of-verlet-or-other-integration-mean">https://math.stackexchange.com/questions/1448005/what-does-the-time-reversibility-of-verlet-or-other-integration-mean</a> also both state velocity Verlet is time reversible</p> <p><a href="https://www5.in.tum.de/lehre/vorlesungen/sci_compII/ss13/uebungen/blatt9solution.pdf" rel="nofollow noreferrer">https://www5.in.tum.de/lehre/vorlesungen/sci_compII/ss13/uebungen/blatt9solution.pdf</a> uses the substitution outlined earlier to show velocity Verlet is time reversible</p> https://scicomp.stackexchange.com/q/36352 0 Solve non-linear equation in R camibc https://scicomp.stackexchange.com/users/37545 2020-11-23T05:50:42Z 2021-10-20T04:30:23Z <p>I need to solve the following equation for <span class="math-container">$x$</span> in [0, 1]. Assume <span class="math-container">$0&lt;\alpha&lt;1$</span> and <span class="math-container">$0&lt;\lambda$</span>.</p> <p><span class="math-container">$$(1 - x)^{\alpha+1} - \lambda (x+1)^{\alpha+1} = -2\lambda (\alpha + 1) x^\alpha$$</span></p> <p>Would very much appreciate any kind of help!</p> https://scicomp.stackexchange.com/q/35211 2 Determining the voxels between two boundary surfaces dimst23 https://scicomp.stackexchange.com/users/36042 2020-05-23T18:16:14Z 2021-10-17T08:02:30Z <h1>Issue description</h1> <p>I am working on human brain tACS simulations where I have the models of the skin, skull, csf, brain and ventricles in STL format. The shape does not matter and there are no intersections. I want to run FE analysis and use SfePy but the only problem I have so far is how to define the volume regions. Using pymesh I am able to generate a pretty detailed and good mesh using the following code:</p> <pre class="lang-py prettyprint-override"><code>import pymesh skin_stl = pymesh.load_mesh('skin.stl') skull_stl = pymesh.load_mesh('skull.stl') csf_stl = pymesh.load_mesh('csf.stl') brain_stl = pymesh.load_mesh('brain.stl') model = pymesh.merge_meshes( (skin_stl, skull_stl, csf_stl, brain_stl) ) model_tet = pymesh.tetrahedralize(model, 10) </code></pre> <p>The meshing I get from the code above can be seen in the following picture. <a href="https://i.stack.imgur.com/fr0rU.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/fr0rU.png" alt="enter image description here" /></a></p> <p>As you can see the boundaries are clear and the meshing algorithm took that into consideration. If I get the number of components of this mesh it is 4 as it should be since I have 4 bounding surfaces. Using the following code I am able to separate and get the vertices and faces of the above tessellated mesh.</p> <pre class="lang-py prettyprint-override"><code>new_mesh = pymesh.form_mesh(model_tet.vertices, model_tet.faces) sep_mesh = pymesh.separate_mesh(new_mesh) </code></pre> <p>I can not in any way separate the mesh and get the voxels included in between two consecutive bounding surfaces. Is there a way to do that using pymesh or something else?</p> https://scicomp.stackexchange.com/q/34952 2 Implementation of Z^2 error estimator in Abaqus for adaptive mesh refinement fruitiest Punch https://scicomp.stackexchange.com/users/32245 2020-04-27T22:11:46Z 2021-10-16T13:23:39Z <p>Currently, I am working on a remeshing routine for my simulations (Abaqus 6.14-1) using python scripts. The simulation deals with the Brinell indentation test and as the remeshing software I use Gmsh 4.5.2 (although the newest probably also works fine).</p> <p>The main script is given by</p> <pre><code>main.py: createIndentationTest.py runJob.py #(only one small step) odbExport.py writeGeo.py #(.geo files are Gmsh readable geometry files) remeshGeo.py #(here occurs the change from quadrangular elements to triangular ones) importInAbaqus.py runJob.py #(only one small step) while totalTime-stepTime&gt;0: odbExport.py writeGeo.py remeshGeo.py importInAbaqus.py runJob.py #(only one small step) totalTime-=stepTime </code></pre> <p>As an example, my last routine created 6 separate jobs, where after every step the mesh was exported to Gmsh, remeshed and imported back in Abaqus to continue the simulation. The first job starts with quadrangular elements (element size of 0.25) for one step. Then the remaining 5 jobs use triangular elements with a custom density function, as can be seen in the following figure (The reason for switching from quadrangular elements to triangular ones during the remeshing process is that it was easier for me to work with in Gmsh).</p> <p><a href="https://i.stack.imgur.com/c3UE1.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/c3UE1.png" alt="Figure 1"></a></p> <p>Now, part of my thesis is to determine the error during the whole simulation and remeshing process and I have been given the task to use the <span class="math-container">$Z^2$</span> error estimator. When compared to my reference simulation, where I ran a whole simulation from start to finish with quadrangular elements with the element size of 0.5, I can see there is a slight discrepancy when analyzing the stresses, as shown in Figure 2.</p> <p><a href="https://i.stack.imgur.com/uETbg.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/uETbg.png" alt="Figure 2"></a></p> <p>For now, I assume that the results on the left-hand side are more accurate, therefore I can assume that there is an error (or deviation) in my simulation results on the right-hand side where the remeshing was implemented. For the past 2 weeks, I have studied some works on the <span class="math-container">$Z^2$</span> error estimator, but I fail to see how to implement this in my script. Therefore, I wanted to ask, if anyone knows how I can apply the error estimator in my python code. </p> <p>I know that I have to export some data from the COORD <code>fieldOutput</code> and use these to solve this equation</p> <p><span class="math-container">$F(\underline{a})=\sum_{i=0}^n\,\left(\sigma_h(x_i,y_i)-\sigma_p^*(x_i,y_i)\right)^2=\sum_{i=0}^n\,\left(\sigma_h(x_i,y_i)-\underline{P}(x_i,y_i)\,\underline{a}\right)^2$</span></p> <p>where <span class="math-container">$\underline{a}$</span> is a set of unknowns, <span class="math-container">$\sigma_h$</span> is the approximate solution to <span class="math-container">$\sigma$</span>, <span class="math-container">$\underline{P}$</span> is a set of appropriate polynomial terms and <span class="math-container">$\sigma_p^*=\underline{P}\,\underline{a}$</span>.</p> <p>Another question I have stems from the fact that I change element types during the simulation process. To estimate the error, the <span class="math-container">$Z^2$</span> method uses so-called <em>superconvergent</em> points. But I have no clue if there will arise a problem when using these since I change my element type once during the whole simulation process.</p> <p>I hope, someone can help me with this.</p> https://scicomp.stackexchange.com/q/34299 0 Global to local transformation matrix in 2D frame structures user32251 https://scicomp.stackexchange.com/users/32251 2020-01-27T00:48:29Z 2021-10-19T19:06:05Z <p>In section 3.2 of <a href="https://doi.org/10.1007/s11081-013-9225-7" rel="nofollow noreferrer">this paper</a> , where 2D planar frame structures are being analyzed, the authors mentioned a transformation matrix to be used in extracting the element displacement vector from the displacement vector of the ground structure in global coordinates.</p> <p><span class="math-container">$u_i^e = (u_x^{(1)},u_y^{(1)}, \theta^{(1)},u_x^{(2)},u_y^{(2)}, \theta^{(2)})$</span> is the element <span class="math-container">$i$</span>'s displacement vector in local coordinates. <span class="math-container">$$u_i^e=T_i u$$</span> where <span class="math-container">$T_i \in \mathbb{R}^ {6\times d}$</span>. What should be the matrix <span class="math-container">$T_i$</span>?</p> <p><strong>Note:</strong> The same matrices are used to assemble the global stiffness matrix from elements stiffness matrix in local coordinates using the following formula: <span class="math-container">$$K= \sum_{i\in E}\sum_{j=1}^3 k_{ij}b_{ij}b_{ij}^\top$$</span> where <span class="math-container">$b_{i1}, b_{i2},b_{i3}\in \mathbb{R}^d$</span> are defined by (<span class="math-container">$\hat{b_{ij}}$</span> are <span class="math-container">$6\times 1$</span> constant vectors):</p> <p><span class="math-container">$$b_{ij}=T_i^\top \hat{b_{ij}}$$</span></p> <p>I am not an expert in statics or mechanical concepts but I need these relations for the optimization problem that I am using for my research.</p> <p> Kureta, Rui, and Yoshihiro Kanno. &quot;A mixed integer programming approach to designing periodic frame structures with negative Poisson’s ratio.&quot; Optimization and Engineering 15.3 (2014): 773-800.</p> https://scicomp.stackexchange.com/q/33454 2 How to avoid gsl root finder evaluate function outside its domain HD189733b https://scicomp.stackexchange.com/users/32861 2019-09-20T08:24:32Z 2021-10-19T23:39:03Z <p>When I use the <a href="https://www.gnu.org/software/gsl/doc/html/roots.html#c.gsl_root_fdfsolver_newton" rel="nofollow noreferrer">newton's method</a> or <a href="https://www.gnu.org/software/gsl/doc/html/multiroots.html#c.gsl_multiroot_fsolver_hybrids" rel="nofollow noreferrer">hybrid solver</a> in the GSL package to deal with 1-D or multidimensional root solving problems, the code frequently crashes when the solver requests function value outside its domain of definition. For example, to solve <span class="math-container">$f(x)=0$</span> where <span class="math-container">$f(x)$</span> is only defined at <span class="math-container">$x\geq0$</span>, even when I start with an initial guess very close to the root, the solver may still ask to evaluate the <span class="math-container">$f(x)$</span> at a negative <span class="math-container">$x$</span>. Bisection method may be helpful in the 1-D problem, but it won't work for a multidimensional problem.</p> <p>I typically try to solve the problem by arbitrarily define the <span class="math-container">$f(x)$</span> in the whole domain. But in some situation, especially in the case of complicated multidimensional root solving, I feel it's hard to extend the domain of the function and make sure it roughly maintain the general trend of the original function.</p> <p>So I'm wondering if there is a method to restrict the region that solver would evaluate the function.</p> https://scicomp.stackexchange.com/q/32577 1 Split-step Fourier method applied on Schrodinger equation decarat https://scicomp.stackexchange.com/users/31666 2019-05-04T19:46:15Z 2021-10-18T02:07:25Z <p>I'm trying to solve a Schrodinger equation of the form <span class="math-container">$i\frac{\partial}{\partial t}\psi=-\frac{\partial^2}{\partial x^2}\psi + (V(x)+\alpha|\psi|^2)\psi$</span> using the split-step Fourier method implemented via MatLab code. In order to make sure it works, I'm testing my code with <span class="math-container">$V(x)=x^2,\alpha=0$</span> and <span class="math-container">$\psi=Ce^{-x^2/2}$</span> where <span class="math-container">$C$</span> is the normalization coefficient (it's set to 1 in this code because it shouldn't matter for testing purposes). My MatLab code is as follows: </p> <pre><code>N = 100000; % Number of Fourier mode dt = .001; % Time step tfinal = 5; % Final time M = round(tfinal/dt); % Total number of time steps L = 5; % Total space length h = L/N; % Space step n =( -N/2:N/2-1); % Indices x = 2*n*h; % Grid points u_i = exp(-((x-0).^2)/2); % Intial pulse u = u_i; % Make a duplicate of the initial pulse in order to compare at the end. k = 2*pi/N * n; % wavenumber values. epsilon = 0; % nonlinear coefficient optical_potential=.5*(x-0).^2; figure(1); % Plotting the probability function of the initial wavefunction (not normalised) plot(x,abs(u_i).^2); hold on; plot(x,optical_potential); hold off; for m = 1:1:M/2 % Start time loop (M/2 since each loop is 2 time steps) c = (fft(u)); % Half time-step linear propagation c = exp(-dt/2*1i*k.^2).*c; u = ifft(c); % Full time-step nonlinear propagation u = exp(dt*1i*(optical_potential + epsilon *(abs(u).^2))).*u; c = (fft(u)); % Half time-step linear propagation c = exp(-dt/2*1i*k.^2).*c; u = ifft(c); end u_out=u; figure(2); % Plotting the probability function of the final wavefunction (not normalised) plot(x,abs(u_out).^2); hold on; plot(x,optical_potential); hold off; error_percent=sum(abs(((abs(u_i) - abs(u_out)))))/sum(abs((abs(u_i)))) </code></pre> <p>The final and initial plots should be the same since <span class="math-container">$\psi$</span> is an eigenstate of the potential (harmonic oscillator). My question is: <span class="math-container">$\text{What is the correct value for k?}$</span> I'm fairly certain what I currently have as my k value is incorrect, despite it giving the correct answer within an error on the magnitude of 1e-13. I've tried a few, and a few give the correct answer within reasonable error, so I'm not sure which is the actual correct value.</p> https://scicomp.stackexchange.com/q/30023 0 Is it possible to partition 2D data into bins such that each bin contains the same number of samples? Zac https://scicomp.stackexchange.com/users/25837 2018-08-10T21:48:30Z 2021-10-15T20:05:23Z <p>I am trying to sort data following a bivariate distribution into a numpy histogramdd, where each bin should contain the same number of data points (to the nearest whole sample).</p> <p>I expect that some kind of quantile-approach is required, and have tried <code>scipy.stats.mstats.mquantiles</code>, which according to the documentation takes up to 2D data. However, it seems to take the dimensions independently, splitting each dimension into to equal marginal probabilities, which doesn't achieve the desired result of equal-probability bins in 2D.</p> <p>Is there a built-in way in scipy/numpy or another package to achieve this (in 2D or higher)? If not, are there algorithms designed to achieve this which I can implement myself directly? </p> https://scicomp.stackexchange.com/q/28188 2 Is there a simple way to avoid carbuncles for FD WENO methods? omican https://scicomp.stackexchange.com/users/25060 2017-11-02T11:27:48Z 2021-10-19T06:00:23Z <p>I have implemented finite-difference WENO scheme for Euler equations (with some variants - WENO-JS, WENO-Z, WENO-M, different flux splitting). It works well, but have problem with so-called carbuncles (shock instability in the areas where its front is almost aligned with the grid). This problem seems to be well-studied for finite-volume methods (quick search shows that there are some fixes for Riemann solvers), but I can't find almost any information regarding carbuncles in finite-difference methods. </p> <p>Are there any known fixes for finite-difference methods?</p> https://scicomp.stackexchange.com/q/27200 5 Algorithms for computing winding numbers of 2-sphere maps Michael Seifert https://scicomp.stackexchange.com/users/24606 2017-06-19T21:41:23Z 2021-10-13T20:06:04Z <p>I have a question concerning computational geometry which arises in the simulation of fields with topological defects, and I'd like to know whether there's an efficient algorithm (or any algorithm) to solve it. </p> <p>The problem is basically the following: consider a grid cell $\mathcal{C}$ in a 3-D cubical grid. On each of the eight vertices of this grid cell, we are given a triplet of numbers $(a_i, b_i, c_i)$ such that $a_i^2 + b_i^2 + c_i^2 = 1$ for each $i = 1 \dots 8$. By triangulating the surface of the cube (which we denote $\partial \mathcal{C}$) and linearly interpolating between the vertices of each triangle in $\partial \mathcal{C}$, we can define a map $f : \partial \mathcal{C} \to \mathbb{R}^3$. Let us assume that the image of this map does not contain 0, i.e., the map so defined is actually from $f: \partial \mathcal{C} \to \mathbb{R}^3 \setminus\{0\}.$</p> <p>We have thus defined a map from a space homeomorphic to $S^2$ ($\partial \mathcal{C}$) to a space that has a non-trivial second homotopy group $\pi_2$. This map may or may not be contractible, and in fact it should have some notion of a winding number associated with it. My questions are:</p> <ul> <li>Is there an algorithm that, given the vertex values $(a_i, b_i, c_i)$ for a triangulated cube, calculates a winding number for the map $f$ so defined?</li> <li>Is there an algorithm that, given the vertex values $(a_i, b_i, c_i)$ for a triangulated cube, calculates whether the map $f$ so defined is nulhomotopic? (This is obviously a weaker question than the first one, but if this is all that can be done I wouldn't be too disappointed.)</li> </ul> <p>I would not be surprised if this question has been addressed in the literature somewhere, but it's quite hard to google for it: every reference to "winding number" seems to be for maps of 1-D curves into some space rather than 2-D spheres.</p> https://scicomp.stackexchange.com/q/14271 7 Tanh-sinh quadrature numerical integration method converging to wrong value ThomasNicholas https://scicomp.stackexchange.com/users/9828 2014-07-27T23:11:13Z 2021-10-15T11:11:11Z <p>I'm trying to write a Python program to use Tanh-sinh quadrature to compute the value of \begin{equation} \int_{-1}^1 \frac{dx}{\sqrt{1-x^2}} \end{equation} but although the program converges to a sensible value with no errors in every case, it's not converging to the correct value (which is $\pi$ for this particular integral) and I can't find the problem. </p> <p>Instead of asking for a desired level of accuracy, the program asks for the number of function evaluations wanted, to make comparisons of convergence with simpler integration methods easier. The number of evaluations needs to be an odd number as the approximation used is \begin{equation} \int_{-1}^1 f(x) dx = \sum_{k=-n}^n w_k f(x_k) \end{equation} Can anyone suggest what I might have done wrong?</p> <pre><code>import math def func(x): # Function to be integrated, with singular points set = 0 if x == 1 or x == -1 : return 0 else: return 1 / math.sqrt(1 - x ** 2) # Input number of evaluations N = input("Please enter number of evaluations \n") if N % 2 == 0: print "The number of evaluations must be odd" else: print "N =", N # Set step size h = 2.0 / (N - 1) print "h =", h # k ranges from -(N-1)/2 to +(N-1)/2 k = -1 * ((N - 1) / 2.0) k_max = ((N - 1) / 2.0) sum = 0 # Loop across integration interval while k &lt; k_max + 1: # Compute abscissa x_k = math.tanh(math.pi * 0.5 * math.sinh(k * h)) # Compute weight numerator = 0.5 * h * math.pi * math.cosh(k * h) denominator = math.pow(math.cosh(0.5 * math.pi * math.sinh(k * h)),2) w_k = numerator / denominator sum += w_k * func(x_k) k += 1 print "Integral =", sum </code></pre>