Recent Questions - Computational Science Stack Exchange most recent 30 from scicomp.stackexchange.com 2020-11-26T02:18:41Z https://scicomp.stackexchange.com/feeds https://creativecommons.org/licenses/by-sa/4.0/rdf https://scicomp.stackexchange.com/q/36375 0 Converting for loop from matlab to python New2Python https://scicomp.stackexchange.com/users/36389 2020-11-25T17:51:04Z 2020-11-25T17:55:21Z <p>I am converting some MATLAB code in to python and have the encountered the error &quot;ValueError Traceback (most recent call last) in 1 for ig in range(nbas): ----&gt; 2 psi[:,ig] = np.polyval(np.array(pp[ig,:ig]),nodes)</p> <p>ValueError: could not broadcast input array from shape (56,1) into shape (56) &quot; I cannot find the issue with following code:</p> <pre><code>import numpy as np import scipy.special as scl ##basis parameters (nbas&lt;ngrid+1) nbas = 54 #Basis set size in one dimension ngrid = 56 #quadrature size #Gauss-Hermite-Quadrature #Hermite matrix def hermipol(n): p = np.zeros((n+1,n+1)) p = 1 if n == 0: p = np.array([[1,0],[2,0]], dtype = float) if n &gt; 0: p[range(0,2)] = np.array([2,0]) if n &gt;=1: for k in range(2,n+1): p[k][range(n)] = 2*p[k-1][range(0,n)] p[k][range(2,n+1)] += -2*(k-1)*p[k-2][range(0,n-1)] for i in range(0,n+1): p[i,:] /= np.sqrt(np.sqrt(np.pi)*2**(i)*scl.factorial(i)) return(p) #Generation Gauss-Hermite Quadrature nodes and weights def ghquad(n): return(np.polynomial.hermite.hermgauss(n)) [nodes,weights] = ghquad(ngrid) pp = hermipol(ngrid) #Evaluation and store psi_i()(x,y)_k) nodes = np.array([nodes]).T weights = np.array([weights]).T print(weights) print(nodes) np.shape(nodes) psi = np.zeros((ngrid,nbas)) np.shape(psi) for ig in range(nbas): psi[:,ig] = np.polyval(np.array(pp[ig,:ig]),nodes) </code></pre> <p>Any guidance would be greatly appreciated.</p> <p>Thanks!</p> https://scicomp.stackexchange.com/q/36367 0 performance comparison between PETSc and SLATE user2348209 https://scicomp.stackexchange.com/users/13113 2020-11-25T05:08:52Z 2020-11-25T05:08:52Z <p>We want to start a new project to solve a large-scale inverse problem (O(10^6) number of parameters) to invert for subsurface wave speeds. We will use FEM to solve forward and adjoint PDEs. In our project, we will need to have both sparse and dense solvers. Previously in similar projects, we were using PETSc library but today when I was searching for new libraries, <a href="https://bitbucket.org/icl/slate/src/master/#markdown-header-documentation" rel="nofollow noreferrer">SLATE library</a> came to my attention (<a href="https://xsdk.info/packages/" rel="nofollow noreferrer">link</a>). I was just wondering if there is any comparison between these two libraries from the performance point of view. Perhaps to make this question richer, I know Trilinos is another competitive library, but based on my personal experience and research that I did, I found PETSc is faster while Trilinos is more reliable. Please comment on this if you have any objections.</p> https://scicomp.stackexchange.com/q/36366 -1 Integer $u$ is known to divide integer $v$. Division algorithms for computing $v÷u$? CopyPasteIt https://scicomp.stackexchange.com/users/37556 2020-11-25T00:35:27Z 2020-11-25T00:35:27Z <p>Integer <span class="math-container">$u$</span> is known to divide integer <span class="math-container">$v$</span>. Division algorithms for computing <span class="math-container">$v÷u$</span>?</p> <p>Let <span class="math-container">$0 \lt u \lt v$</span> are numbers such that <span class="math-container">$u \mid v$</span>.</p> <blockquote> <p>Are there algorithms (computer technology) that can use this information to perform the quotient calculation <span class="math-container">$\Large \frac{v}{u}$</span>?</p> </blockquote> <p><strong>My work</strong></p> <p>I've been thinking about (as a novice) 'residue machines' and can imagine calculating residues</p> <p><span class="math-container">$\quad a \pmod{n}$</span></p> <p>in a cryptographic application without calculating, at first, the quotient.</p> <p>If later we need to divide, say <span class="math-container">$m$</span> into, say <span class="math-container">$b$</span>, and it can be determined that</p> <p><span class="math-container">$\quad b = 0 \pmod{m}$</span></p> <p>can computing cycles be reduced to find the integer <span class="math-container">$k$</span> solving <span class="math-container">$b = km$</span>?</p> https://scicomp.stackexchange.com/q/36365 1 Asymptotic complexity of fixed-rank SVD mana https://scicomp.stackexchange.com/users/37002 2020-11-24T22:48:59Z 2020-11-25T11:27:51Z <p>According to the W<a href="https://en.wikipedia.org/wiki/Singular_value_decomposition#Calculating_the_SVD" rel="nofollow noreferrer">ikipedia article on Singular Value Decomposition</a>, the asymptotic complexity of computing the SVD of an arbitrary m×n matrix M with m&gt;n by the popular Householder QR methods is O(mn2). Are there any algorithms (perhaps Householder QR) that provide better asymptotic guarantees for fixed-rank matrices?</p> <p>In other words: let Sn,k be the collection of n×n matrix of rank k. Are there algorithms that provide a better asymptotic complexity of computing the SVD of elements of Sn,k as n→∞ than O(n3)?</p> https://scicomp.stackexchange.com/q/36362 0 Tool to compare if two logical expressions are same! user50121 https://scicomp.stackexchange.com/users/37554 2020-11-24T17:46:09Z 2020-11-25T22:56:06Z <p>We have challenge in my current assignment where we need to modify/minimize an existing logical expression to another new logical expression. But the result should be the same.</p> <p>For eg: the ask to convert</p> <p>(1 AND (2 OR 9 OR 10) AND (3 OR 4 OR 5) AND (6 OR 7) AND 8) OR (1 AND (2 OR 9 OR 10) AND (11 OR 12) AND 8)</p> <p>into</p> <p>1 AND (2 OR 9 OR 10) AND 8 AND (((3 OR 4 OR 5) AND (6 OR 7)) OR 11 OR 12)</p> <p>and see if both of them evaluate to same, across a various combinations of ( T or F) for each number.</p> <p>I need to the same activity for almost 100+ expressions. hence i am looking for a tool to help expedite.</p> https://scicomp.stackexchange.com/q/36369 0 A better word to indicate slowness/high latency? bd3lk https://scicomp.stackexchange.com/users/0 2020-11-24T09:49:41Z 2020-11-25T09:43:18Z <p>We are comparing two techniques in computer science. We want to say X has <strong>&quot;significantly high latency&quot;</strong> when executed on system Y.</p> <p>Is there a better one-word term we can use for the above to mean the same as 'significantly high latency'? e.g., slow &lt;-- but a better word than slow?</p> https://scicomp.stackexchange.com/q/36360 2 Does incomplete LU preconditioning improve the asymptotic scaling of Krylov subspace methods? gTcV https://scicomp.stackexchange.com/users/21342 2020-11-24T02:46:40Z 2020-11-24T02:46:40Z <p>It is well known that unpreconditioned Krylov subspace methods applied to the finite-difference-discretised Poisson equation with <span class="math-container">$n$</span> grid points per direction require <span class="math-container">$O(n \, |\log(\varepsilon)|)$</span> iterations to reduce the error in the initial guess by a factor <span class="math-container">$\varepsilon \in (0,1)$</span>. This is typically a prohibitively large number of iterations; hence it is often suggested to reduce the number of iterations through incomplete LU (ILU) preconditioning.</p> <p>My question is: does ILU reduce the <span class="math-container">$\text{# iterations} = O(n \,|\log(\varepsilon)|)$</span> estimate asymptotically (e.g. to <span class="math-container">$\text{# iterations} = O(\sqrt{n} \,|\log(\varepsilon)|)$</span>), or does it merely reduce the prefactor implied by the big-O notation?</p> <p>My own experiments seem to suggest that ILU preconditioning only improves the prefactor, but maybe I am not using the method properly?</p> https://scicomp.stackexchange.com/q/36359 1 Which 2D PDE with an exact solution can I use to test/verify my FEM-PDE code? Marten https://scicomp.stackexchange.com/users/20541 2020-11-23T23:37:58Z 2020-11-24T21:57:29Z <p>I have created a program to solve 2D, time-dependent PDEs with the finite element method and get reasonable looking results for the 2D acoustic wave equation. Now I would like to go further and solve a PDE with a known exact/analytic solution to compare against. However, I have a lot of trouble finding a suitable equation. It seems that there are no 2D, time-dependent equations with an analytic solution that does not involve infinite sums or the like. I thought about taking a 1D equation and extend it to two dimensions by just solving it on a 2D domain without changing anything, in effect replicating the equation along the y-axis. However, boundary conditions seem to mess things up and these 1D equations work on infinite domains, i.e. the whole real axis.</p> <p>Is there such an equation for me to use? Is there a different established way of testing the correctness of a PDE solver implementation?</p> https://scicomp.stackexchange.com/q/36357 0 System of second order ODEs Runge Kutta 4th order J Wright https://scicomp.stackexchange.com/users/37349 2020-11-23T14:22:18Z 2020-11-23T16:21:13Z <p>I have two masses, and two springs which are atattched to a wall on one end. For context I attached the system of equations. <a href="https://i.stack.imgur.com/ctOnE.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/ctOnE.png" alt="enter image description here" /></a></p> <p>How can I add a command to see my analytic solution for my system of equations? I want to be able to compare the numeric vs analytic. I think it would be something like (t, answer, t, answer). Thank you.</p> <pre><code>import numpy as np import matplotlib.pyplot as plt def f(x,t): k1=20 k2=30 m1=3 m2=5 return np.array([x, (-k1*x-k2*(x - x))/m1, x, (-k2*(x-x)/m2)]) h=.01 t=np.arange(0,15+h,h) y = np.zeros((len(t), 4)) y[0, :] = [1, 0, 0, 0] for i in range(0,len(t)-1): k1 = h * f( y[i,:], t[i] ) k2 = h * f( y[i,:] + 0.5 * k1, t[i] + 0.5 * h ) k3 = h * f( y[i,:] + 0.5 * k2, t[i] + 0.5 * h ) k4 = h * f( y[i,:] + k3, t[i+1] ) y[i+1,:] = y[i,:] + ( k1 + 2.0 * ( k2 + k3 ) + k4 ) / 6.0 plt.plot(t,y[:,0],t,y[:,2]) plt.gca().legend(('x_1','x_2')) plt.show() <span class="math-container"></span> </code></pre> https://scicomp.stackexchange.com/q/36355 0 Are spatial boundary conditions required for PDEs discretized with Method of Lines? Marten https://scicomp.stackexchange.com/users/20541 2020-11-23T11:14:52Z 2020-11-23T12:39:07Z <p>As far as I understand, you need to define boundary conditions in time and space to select a unique solution to a PDE and make it solvable. However, in ODEs I only need to specific the initial value in time which I think of as a boundary condition. If I now use the method of lines to transform a PDE into a system of ODEs, do I need the spatial boundary conditions? If yes, are they imposed by fixing some of the &quot;lines&quot; to given functions? If no, how can it be that they are first required but then somehow the need for them disappears?</p> <p>I got this question because I did this with the 2D acoustic wave equation and the result looks reasonable, yet I did not impose any spatial boundary conditions. I have also read somewhere (that I sadly cannot find anymore) that there are implicit boundary conditions. Does that mean that me not specifying any conditions is equivalent to choosing some implicitly?</p> https://scicomp.stackexchange.com/q/36354 0 Help Installing MFEM WizardGeneral https://scicomp.stackexchange.com/users/37546 2020-11-23T06:39:40Z 2020-11-23T06:39:40Z <p>I am fairly new to both FEM and C++ (I have used FEM in SolidWorks and C++ for a couple of projects) and I am hoping to use MFEM for an existing Computational Astrophysics Project.</p> <p>I am unable to install even using the simple &quot;straightforward&quot; serial build instructions available on the MFEM webpage.</p> <p>I have downloaded mfem-4.2 and glvis3.4 and run the command:</p> <p><code>make serial -j</code></p> <p>which seems to work, but when I run:</p> <p><code>make MFEM_DIR=../mfem-4.2 -j</code></p> <p>I get a series of fatal errors, such as <code>openglvis.cpp:12:10: fatal error: 'GL/gl.h' file not found</code></p> <p>simply typing make or any other <code>make</code> command results in <code>makefile:192: *** The MFEM library is not built. Stop.</code></p> <p>I have no idea what I'm doing wrong. I have hit many dead ends trying to get a usable FEM library for C++ operational.</p> <p>Can I somehow configure the MFEM files for convenient use in Xcode?</p> <p>Any advice or direction appreciated.</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 2020-11-23T16:43:07Z <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/36350 -3 Scientific calculator device malfunctioning [closed] Prashant Akerkar https://scicomp.stackexchange.com/users/37431 2020-11-22T18:41:37Z 2020-11-22T18:41:37Z <p><a href="https://en.wikipedia.org/wiki/Scientific_calculator" rel="nofollow noreferrer">https://en.wikipedia.org/wiki/Scientific_calculator</a></p> <p>Can there be scenario/s that a scientific calculator gives incorrect answers in calculations to the end user?<img src="https://i.stack.imgur.com/wsPXs.jpg" alt="enter image description here" /></p> https://scicomp.stackexchange.com/q/36349 0 Temporal/spectral conversion for large fields - best approach? arc_lupus https://scicomp.stackexchange.com/users/20545 2020-11-22T18:18:02Z 2020-11-22T18:18:02Z <p>I am currently working on a more efficient implementation of a pulse propagation algorithm. The propagation is done in the spectral domain, but several evaluations (such as energy calculations) are done in the time domain. The field is either defined in cylindrical coordinates (i.e. <span class="math-container">$E=E(r, t)$</span>) or in cartesian coordinates (i.e. <span class="math-container">$E=E(x, y, t)$</span>).</p> <p>To convert the field from the temporal domain into the spectral domain, I can either apply the hankel transformation for converting the field from <span class="math-container">$r$</span>-space to <span class="math-container">$k$</span>-space, and afterwards a 1d-FFT for converting the field from the <span class="math-container">$t$</span>-space into the <span class="math-container">$\omega$</span>-space (for the cylindrical coordinate system), or I can first apply a 2d-FFT for converting the field from the <span class="math-container">$x, y$</span>-space into the <span class="math-container">$k_x,\,k_y$</span>-space, and afterwards a 1d-FFT for going from the <span class="math-container">$t$</span>-space into the <span class="math-container">$\omega$</span>-space. The Hankel transformation can be reduced into a matrix-matrix-multiplication.</p> <p>In my current implementation I am calculating everything on a GPU, while using <code>ArrayFire</code> as library in the backend. Still, for higher resolutions and longer pulses the memory of my GPU is not enough, and I have to switch to using the CPU instead. According to benchmarks, running the calculations on a single node this might increase the runtime by a factor of 5-20, depending on parameters. Therefore, I was intending to use several nodes together (using MPI) to reduce that factor at least a bit. Would that make sense?</p> <p>Unfortunately, <code>ArrayFire</code> does not support such calculations (zgemm/n-FFT) via MPI, but I also do not want to use several different libraries just depending on the target device. Are there libraries which support those operations both for GPUs, local threads on a CPU node and MPI for several nodes? I have seen that PETSc might support that, but I am not sure about that. Are there others? Or would another approach be easier here?</p> https://scicomp.stackexchange.com/q/36348 0 How to make a directed graph symmetric? IPribec https://scicomp.stackexchange.com/users/37438 2020-11-22T17:22:59Z 2020-11-22T17:22:59Z <p>Say I have a directed graph given as an adjacency matrix <span class="math-container">$A$</span> in <a href="https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_row_(CSR,_CRS_or_Yale_format)" rel="nofollow noreferrer">CSR format</a> represented by the arrays <code>ia</code> (row indexes) and <code>ja</code> (column indexes). In my application the graph serves as the underlying spatial connectivity for an <a href="https://github.com/treverhines/RBF" rel="nofollow noreferrer">RBF-FD</a> discretization of a PDE. The graph originates from finding the <span class="math-container">$N$</span>-nearest neighbor stencil for each node <span class="math-container">$\boldsymbol{x}_i$</span> of a <span class="math-container">$d$</span>-dimensional point cloud.</p> <p>For some purposes it is benefical to use symmetric stencils <span class="math-container">$S$</span>, meaning that <span class="math-container">$\boldsymbol{x}_j \in S(\boldsymbol{x}_i)$</span> implies <span class="math-container">$\boldsymbol{x}_i \in S(\boldsymbol{x}_j)$</span> for all <span class="math-container">$i,j$</span> (or at least for interior nodes). This corresponds to making the directed graph symmetric (undirected). One example of where this is needed is in the METIS partitioning library, which expects graphs to be symmetric. Similarly, the Reverse Cuthill McKee ordering algorithm also expects a symmetric matrix.</p> <p><strong>Question:</strong> Given an adjacency matrix in CSR format as arrays <code>ia</code> and <code>ja</code>, how can I find the symmetric graph adjacency matrix arrays <code>ias</code>, <code>jas</code>?</p> <p>I have noted in Scipy, the way they achieve this is by forming the matrix <span class="math-container">$A + A^T$</span> (see Scipy source <a href="https://github.com/scipy/scipy/blob/master/scipy/sparse/csgraph/_reordering.pyx" rel="nofollow noreferrer">here</a>). Note that I don't want to compute the actual matrix values, but only the structure of the resulting matrix.</p> <p>Is there any other way to achieve this, apart from the naive way of of initializing <span class="math-container">$M$</span> empty adjacency lists (for all <span class="math-container">$M$</span> nodes of my point cloud), iterating through all the nodes, and pushing each connection <span class="math-container">$(i,j)$</span> into the right list? It seems kind of unneccesary to process each connection twice.</p> https://scicomp.stackexchange.com/q/36347 1 Solution of Coupled Differential equation for a 2d linear flow using RK4 method in python 3 Anirban Majumdar https://scicomp.stackexchange.com/users/37541 2020-11-22T16:20:46Z 2020-11-22T19:01:40Z <p>I want to study the dynamics of a 2d linear flow, whose dynamical equation is- <span class="math-container">$\begin{pmatrix} \dot{x_1}\\ \dot{x_2}\\ \end{pmatrix}=\begin{pmatrix} 1 &amp; 1\\ 4 &amp; -2\\ \end{pmatrix}\begin{pmatrix} x_1\\ x_2\\ \end{pmatrix}$</span>. Now I have tried to solve and plot y vs. x of this coupled differential equation using RK4 in python for the initial condition <span class="math-container">$(y_0=2, x_0=-1)$</span>. My code is following, but the graph is not correct [In the graph origin should be a saddle point, <span class="math-container">$\begin{pmatrix} -0.25\\ 1\\ \end{pmatrix}$</span> axis should be stable manifold and <span class="math-container">$\begin{pmatrix} 1\\ 1\\ \end{pmatrix}$</span> axis should be unstable manifold]-</p> <pre><code>import numpy as np from math import sqrt import matplotlib.pyplot as plt # Equations: def V(u,t): x1, x2, v1, v2 = u return np.array([ v1, v2, (x1+x2), -(4*x1-2*x2)]) def rk4(f, u0, t0, tf , n): t = np.linspace(t0, tf, n+1) u = np.array((n+1)*[u0]) h = t-t for i in range(n): k1 = h * f(u[i], t[i]) k2 = h * f(u[i] + 0.5 * k1, t[i] + 0.5*h) k3 = h * f(u[i] + 0.5 * k2, t[i] + 0.5*h) k4 = h * f(u[i] + k3, t[i] + h) u[i+1] = u[i] + (k1 + 2*(k2 + k3) + k4) / 6 return u, t u, t = rk4(V, np.array([2.0, 0., -1, 1.]) , 0. , 1. , 1000) x1, x2, v1, v2 = u.T plt.plot(x1,x2) plt.show() </code></pre> <p><em><strong>Can anyone helps me to write this code.</strong></em></p> https://scicomp.stackexchange.com/q/36344 1 Understanding why scipy.fft.fft (fast Fourier transform) doesn't work as expected PiKindOfGuy https://scicomp.stackexchange.com/users/37540 2020-11-22T10:45:21Z 2020-11-22T11:42:05Z <p>I write the following fast Fourier transform code into my Python notebook expecting to see a plot wherein there's a spike at <span class="math-container">$1/2\pi$</span> since that's the frequency of the sin function, but instead I get the plot below. Any ideas as to what's going on?</p> <pre><code>import math import matplotlib.pyplot as plt import numpy as np from scipy.fft import fft x = np.arange(1000) y = np.sin(x) freq = fft(y) plt.plot(np.real(freq)) plt.show() </code></pre> <p><a href="https://i.stack.imgur.com/2UcYi.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/2UcYi.png" alt="fft plot" /></a></p> https://scicomp.stackexchange.com/q/36342 3 Advantage of diagonal "jitter" for numerical stability? Algo https://scicomp.stackexchange.com/users/13154 2020-11-22T06:50:22Z 2020-11-24T13:52:17Z <p>In a machine learning <a href="https://github.com/mml-book/mml-book.github.io/blob/master/tutorials/tutorial_linear_regression.solution.ipynb" rel="nofollow noreferrer">code</a>, that computes optimum parameters <span class="math-container">$\theta _{MLE}$</span> of a linear regression model, by maximum likelihood estimation:</p> <p><span class="math-container">$$\boldsymbol \theta^\text{ML} = (\boldsymbol\Phi^T\boldsymbol\Phi )^{-1}\boldsymbol\Phi^T\boldsymbol y$$</span></p> <p>Where <span class="math-container">$y$</span> is the target vector and <span class="math-container">$\Phi$</span> is the polynomial feature matrix. In the linked notebook we can find:</p> <blockquote> <p>For reasons of numerical stability, we often add a small diagonal &quot;jitter&quot; <span class="math-container">$\kappa$</span> to <span class="math-container">$\boldsymbol\Phi^T\boldsymbol\Phi$</span> so that we can invert the matrix without significant problems so that the maximum likelihood estimate becomes <span class="math-container">$$\boldsymbol \theta^\text{ML} = (\boldsymbol\Phi^T\boldsymbol\Phi + \kappa\boldsymbol I)^{-1}\boldsymbol\Phi^T\boldsymbol y$$</span></p> </blockquote> <p>In the code, <span class="math-container">$\kappa$</span> is very small value of 1e-08.</p> <p>So, how does the diagonal &quot;jitter&quot; <span class="math-container">$\kappa$</span> affects stability?</p> https://scicomp.stackexchange.com/q/36337 1 Optimization on the multinomial manifolds of stochastic non-square matrices Hu jiaran https://scicomp.stackexchange.com/users/37534 2020-11-21T12:41:07Z 2020-11-22T07:20:10Z <p>Thanks for note! So I have an optimization problem with simple form but the decision variable is a large-scale matrix. My problem is similar to a existing problem <a href="https://scicomp.stackexchange.com/questions/23940/optimization-on-the-manifold-of-stochastic-matrices/36335#36335">here</a> about multinomial manifolds and my objective function is simpler than the counterpart in that problem. The decision variable <span class="math-container">$\mathbf T \in {\bf R}^{1300\times 3100}$</span> is a stochastic matrix (each row sums to 1) with non-negative elements. <span class="math-container">$\mathbf S \in {\bf R}^{600\times 1300}$</span> is a known constant matrix with non-negative elements. <span class="math-container">$\mathbf L \in {\bf R}^{3100\times 9000}$</span> is a known constant sparse matrix. There are exactly one <span class="math-container">$1$</span> and one <span class="math-container">$-1$</span> in each column of <span class="math-container">$\mathbf L$</span>, so the column sum of <span class="math-container">$\mathbf L$</span> is equal to <span class="math-container">$0$</span> for every column. Let <span class="math-container">$\mathbf X$</span> denote the product of above three matrices, i.e. <span class="math-container">$\mathbf X=\mathbf S\mathbf T\mathbf L$</span>. The optimization problem is below. <span class="math-container">$$\text{minimize}\hspace{3mm}f(\mathbf T)=\sum_{i=1}^{600}\sum_{j=1}^{9000}|\mathrm X_{ij}|$$</span> <span class="math-container">$$\text{subject to}\hspace{19mm}\mathbf T{\bf 1}=\bf{1}$$</span> <span class="math-container">$$\hspace{39mm}\mathrm T_{ij}&gt;0$$</span> The optimization function means I expect to obtain a matrix <span class="math-container">$\mathbf X$</span> as similar to zero matrix as possible, so it can be replaced by minimizing Frobenius norm of <span class="math-container">$\mathbf X$</span> or other similar forms. As far as I know, my problem resembles optimal transport in some ways. In addition, there may be some trick to reduce the amount of computation and cost of storage given the special nature of matrix <span class="math-container">$\mathbf L$</span>. I'm doubting whether the manopt toolbox in Matlab <a href="https://www.manopt.org/reference/manopt/manifolds/multinomial/multinomialfactory.html" rel="nofollow noreferrer">here</a> is capable of dealing with this large-scale optimization problem (speed is not so important). I have thought that a python package named pymanopt <a href="https://www.pymanopt.org/" rel="nofollow noreferrer">here</a> might work but the multinomial manifolds is not currently supported by pymanopt. I'm planning to employ a feasible algorithm or a computation environment to solve this optimization task. If there are ideas to make my mission possible, what would be the best approach here? Thanks very much for any advice or comments.</p> https://scicomp.stackexchange.com/q/36334 3 Equations that are easier to verify than to solve? Yaroslav Bulatov https://scicomp.stackexchange.com/users/18786 2020-11-21T01:31:20Z 2020-11-22T02:27:20Z <p>Are there interesting examples of (systems) of equations where it is known to be harder to find a solution (in terms of scaling with respect to problem size) than verifying a provided solution for correctness?</p> <p>A non-expert may find it surprising that Stein, Riccati, Sylvester matrix equations with <span class="math-container">$d\times d$</span> matrices all have the same <span class="math-container">$O(d^3)$</span> complexity for solving as for verifying, wondering if this is a rule that holds more generally.</p> https://scicomp.stackexchange.com/q/36333 2 In Eigen, can a sparse matrix contain vectors/objects instead of simple scalar values? Pietro https://scicomp.stackexchange.com/users/10441 2020-11-20T23:59:20Z 2020-11-21T11:36:36Z <p>I need to have a sparse matrix whose elements are not simple numbers, but objects, e.g. a couple of floating point values and a bunch of integer indices.</p> <p>I am wondering if Eigen has something similar, beyond its <code>Eigen::SparseMatrix&lt; _Scalar, _Options, _StorageIndex &gt;</code> class template (the only sparse matrix I have found so far), which only accepts a scalar template parameter.</p> <p>I know I can use a SoA (Structure of Arrays) approach, but in this specific case, for data locality reasons, I am trying with an AoS (Array of Structures) one.</p> https://scicomp.stackexchange.com/q/36330 0 How to determine guitar tones played as early as possible? Emil https://scicomp.stackexchange.com/users/19947 2020-11-20T19:05:20Z 2020-11-22T08:06:39Z <p>I want to detect chords on a guitar as early as possible, but my approach with a sliding window and a filter bank seems to introduce too much lag.</p> <p>Would required observation time decrease by using a model where there are only a finite number of tones possible and only a finite number of simultaneous tones? (I.e. the different strings of the guitar).</p> <p>I would suppose that for the system to not be underdetermined the number of samples would have to be at least as many as the number of guitar strings, and the time window would have to be at least on the same time scale as the tone with the shortest period. Or maybe the number of samples would have to be at least the same as the dimension of the model space (something like ~6 * 20)? And probably the amplitude resolution of the microphone together with the slowest frequency would set a constraint too?</p> https://scicomp.stackexchange.com/q/36292 0 How to insert a(x) function in non homogeneous parabolic pde for implicit method in Python? Mr.Podilatis https://scicomp.stackexchange.com/users/37062 2020-11-15T20:50:30Z 2020-11-22T18:13:53Z <p>I have the following inhomogeneous parabolic initial/boundary value problem: <span class="math-container">$$u_{t}(t,x) = (1-x^{2})u_{xx}(t,x)+u(t,x),$$</span> for <span class="math-container">$t \in [0,1]$</span> and <span class="math-container">$x \in [-1,1]$</span> <span class="math-container">$$u(0,x) = \sin(\pi x),$$</span> for <span class="math-container">$x \in [-1,1]$</span> initial condition <span class="math-container">$$u(t,-1)=u(t,1)=0,$$</span> <span class="math-container">$t \in [0,1]$</span> Dirichlet boundary conditions.</p> <p>I want to construct a Backward Euler method with <span class="math-container">$N_{x} =39$</span> and <span class="math-container">$Nt=400$</span></p> <p>and a Crank Nicolson method with <span class="math-container">$N_{x} =39$</span> and <span class="math-container">$Nt=20$</span> but I don't know how to put <span class="math-container">$$a = (1-x^2)$$</span> inside my method in the script below.Any help?</p> <p>Backward Euler</p> <pre><code> def g(x): return(np.sin(np.pi*x)) Nx = 39 Nt = 400 L = 1 dx = (L - (-L))/(Nx - 1) t0 = 0 Tf = 1 dt = (Tf - t0)/(Nt - 1) h = (L - (-L))/(Nx+1) t = Tf / Nt m = t/h**2 print(&quot;m =&quot;, round(m)) x = np.linspace(-L, L, Nx+1) t = np.linspace(t0, Tf, Nt+1) a = np.array([1-x**2]).reshape(Nx+1) u = np.zeros(Nx+1) u_n = np.zeros(Nx+1) A = np.zeros((Nx+1, Nx+1)) b = np.zeros(Nx+1) for i in range(1, Nx): A[i,i-1] = -m A[i,i+1] = -m A[i,i] = 1 + 2*m A[0,0] = A[Nx,Nx] = 1 A = A*a #--- initial condition u(x,0) = g(x) for i in range(0, Nx+1): u_n[i] = g(x[i]) for n in range(0, Nt): # Compute b and solve linear system for i in range(1, Nx): b[i] = -u_n[i] b = b[Nx] = 0 u[:] = scipy.linalg.solve(A, b) # Update u_n before next step u_n[:] = u plt.plot(u) plt.show() </code></pre> <p>Crank - Nicolson</p> <pre><code>from scipy.sparse.linalg import spsolve Nx = 39 Nt = 20 L = 1 dx = (L - (-L))/(Nx - 1) t0 = 0 Tf = 1 dt = (Tf - t0)/(Nt - 1) h = (L - (-L))/(Nx+1) t = Tf / Nt m = t/h**2 print(&quot;m =&quot;, round(m)) x = np.linspace(-L, L, Nx+1) t = np.linspace(t0, Tf, Nt+1) # Representation of sparse matrix and right-hand side main = np.zeros(Nx+1) lower = np.zeros(Nx) upper = np.zeros(Nx) b = np.zeros(Nx+1) # Precompute sparse matrix main[:] = 1+m lower[:] = -1/2*m upper[:] = -1/2*m # Insert boundary conditions main = 0 main[Nx] = 0 A = scipy.sparse.diags( diagonals=[main, lower, upper], offsets=[0, -1, 1], shape=(Nx+1, Nx+1), format='csr') A = A*a print(A) # Set initial condition for i in range(0,Nx+1): u_n[i] = g(x[i]) for n in range(0, Nt): b = u_n b = b[-1] = 0.0 # boundary conditions u[:] = scipy.sparse.linalg.spsolve(A, b) u_n[:] = u plt.plot(u) plt.show() <span class="math-container"></span> </code></pre> https://scicomp.stackexchange.com/q/36136 1 Calculating Error for Poisson Equation using Successive Over-Relaxation technique, Python Drishika Nadella https://scicomp.stackexchange.com/users/37306 2020-10-19T20:53:56Z 2020-11-21T13:03:25Z <p>I am trying to solve the Poisson Equation</p> <p><span class="math-container">$\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} = 32(x(x-1) + y(y-1))$</span></p> <p>for a 61x61 grid using Python3 with boundary conditions being <span class="math-container">$T=0$</span> on all four boundaries and taking an initial guess of <span class="math-container">$T=50$</span>. I am considering a relaxation parameter <span class="math-container">$w=1.6$</span>.</p> <p>To solve the Poisson Equation using SOR technique, I discretize using the Finite Difference method and do:</p> <p><span class="math-container">$T^{(n+1)}_{i,j} = 0.25w{(T^{n}_{i+1, j} + T^{n}_{i-1, j} + T^{n}_{i, j+1} +T^{n}_{i, j-1} - 32(x_i(x_i - 1) + y_i(y_i-1))dx^2)} + T^n_{i,j}(1-w)$</span></p> <p>I would like to find the error between the iterations by simply finding the difference between the older iteration and the newer one (by finding the maximum difference in the array), and checking if that error is lesser than my permissible error (<span class="math-container">$10^{-7}$</span>). If it is lesser, the iterations stop, otherwise they continue.</p> <p>However, this method doesn't seem to converge even after 2000 iterations. I would like to know where I am making the mistake (I suspect it is in the error calculation). Why doesn't this simple error calculation work? What changes can I implement to make it work?</p> <p>Here is my code:</p> <pre><code>import numpy as np import matplotlib.pyplot as plt # Set Dimension and delta nx = 61 #grid size my = 61 x = 1.0 #total x length y = 1.0 #total y length dx = x/(nx-1) dy = y/(my-1) xarr = np.linspace(0,x,nx) yarr = np.linspace(0,y,my) print(xarr) &quot;&quot;&quot; Output: [0. 0.01666667 0.03333333 0.05 0.06666667 0.08333333 0.1 0.11666667 0.13333333 0.15 0.16666667 0.18333333 0.2 0.21666667 0.23333333 0.25 0.26666667 0.28333333 0.3 0.31666667 0.33333333 0.35 0.36666667 0.38333333 0.4 0.41666667 0.43333333 0.45 0.46666667 0.48333333 0.5 0.51666667 0.53333333 0.55 0.56666667 0.58333333 0.6 0.61666667 0.63333333 0.65 0.66666667 0.68333333 0.7 0.71666667 0.73333333 0.75 0.76666667 0.78333333 0.8 0.81666667 0.83333333 0.85 0.86666667 0.88333333 0.9 0.91666667 0.93333333 0.95 0.96666667 0.98333333 1. ] &quot;&quot;&quot; # Boundary condition Ttop = 0 Tbottom = 0 Tleft = 0 Tright = 0 # Initial guess of interior grid Tguess = 50 # Set colour interpolation and colour map colorinterpolation = 50 colourMap = plt.cm.jet # Set meshgrid X, Y = np.meshgrid(np.arange(0, nx), np.arange(0, my)) # Set array size and set the interior value with Tguess T = np.empty((nx, my)) T.fill(Tguess) #Boundary conditions T[(my-1):, :] = Ttop T[:1, :] = Tbottom T[:, (nx-1):] = Tright T[:, :1] = Tleft T_init=T print(&quot;The initial matrix is: \n&quot;, T_init) &quot;&quot;&quot; Output: The initial matrix is: [[ 0. 0. 0. ... 0. 0. 0.] [ 0. 50. 50. ... 50. 50. 0.] [ 0. 50. 50. ... 50. 50. 0.] ... [ 0. 50. 50. ... 50. 50. 0.] [ 0. 50. 50. ... 50. 50. 0.] [ 0. 0. 0. ... 0. 0. 0.]] &quot;&quot;&quot; #SOR Technique def SORAlgo(error, w, T, MaxIter): for n in range(MaxIter): Tn=T.copy() #Solving the Poisson Equation using array operations T[1:-1, 1:-1] = w*0.25*((Tn[2:, 1:-1] + Tn[:-2, 1:-1])*dy**2 + (Tn[1:-1, 2:] + Tn[1:-1, :-2])*dx**2 - 32*((xarr[1:-1]*(xarr[1:-1]-1) + yarr[1:-1]*(yarr[1:-1]-1)))*dx**2*dy**2)/(2*(dx**2 + dy**2)) + (1-w)*Tn[1:-1, 1:-1] #Tn will be the older value, T will be the newer value. Finding the max difference in corresponding values of both arrays max_error = (abs(T-Tn)).max() if max_error&lt;error or n==MaxIter-2: print(&quot;The relaxation parameter is: &quot;, w) print(&quot;The number of iterations taken is: &quot;, n) print(&quot;The error is: &quot;, max_error) break return Tn #Taking error = 10^-7, relaxation parameter=1.6 and maximum iterations=2000 T_final = SORAlgo(0.0000001, 1.6, T_init, 2000) print(T_final) #print(&quot;The final matrix is: &quot;, T_final) cp = plt.contourf(X, Y, T_final, colorinterpolation, cmap=colourMap) plt.colorbar() plt.show() </code></pre> <p>And this is my output:</p> <pre><code>The relaxation parameter is: 1.6 The number of iterations taken is: 1998 The error is: 2.35862095815448e-05 [[0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00 0.00000000e+00 0.00000000e+00] [0.00000000e+00 2.25244388e-05 4.41640095e-05 ... 4.41640095e-05 2.25244388e-05 0.00000000e+00] [0.00000000e+00 2.40057069e-05 4.75951601e-05 ... 4.75951601e-05 2.40057069e-05 0.00000000e+00] ... [0.00000000e+00 2.40057069e-05 4.75951601e-05 ... 4.75951601e-05 2.40057069e-05 0.00000000e+00] [0.00000000e+00 2.25244388e-05 4.41640095e-05 ... 4.41640095e-05 2.25244388e-05 0.00000000e+00] [0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00 0.00000000e+00 0.00000000e+00]] </code></pre> <p>Notice how it goes until nearly 2000 iterations and stops only because I explicitly ask the loop to break at 1999 iterations and the error is higher than specified (<span class="math-container">$10^{-7}$</span>). This is a link to the plot since I am unable to directly paste it here:</p> <p><a href="https://drive.google.com/file/d/1xLbVSp9XA92saZr26b0gMXYZZ_YsR2GQ/view?usp=sharing" rel="nofollow noreferrer">https://drive.google.com/file/d/1xLbVSp9XA92saZr26b0gMXYZZ_YsR2GQ/view?usp=sharing</a></p> <p>Thank you.</p> <p>Edit: I took the suggestions in the comments.</p> <p>I tweaked the source term and simplified the equation, so instead of</p> <pre><code>Tn[1:-1, 1:-1] = w*0.25*(Tn[2:, 1:-1] + Tn[:-2, 1:-1] + Tn[1:-1, 2:] + Tn[1:-1, :-2])/(dx**2) - 32*((xarr[1:-1]*(xarr[1:-1]-1) + yarr[1:-1]*(yarr[1:-1]-1)))*dx**2*dy**2 + (1-w)*Tn[1:-1, 1:-1] </code></pre> <p>I did:</p> <pre><code>b = np.zeros((my, nx)) #The new source term for i,j in zip(range(nx), range(my)): b[i, j] = 32*(xarr[i]*(xarr[i]-1) + yarr[j]*(yarr[j]-1)) ... T[1:-1,1:-1] = (1-w)*Tn[1:-1, 1:-1] + w*0.25*(T[1:-1, 2:] + T[1:-1, :-2] + T[2:, 1:-1] + T[:-2, 1:-1] - dx**2*b[1:-1, 1:-1]) </code></pre> <p>When I make the source term 0 and take the relaxation parameter as 1, i.e.</p> <pre><code>T[1:-1,1:-1] = (1-w)*Tn[1:-1, 1:-1] + w*0.25*(T[1:-1, 2:] + T[1:-1, :-2] + T[2:, 1:-1] + T[:-2, 1:-1]) #w=1 </code></pre> <p>I get the correct plot and the solution converges at 349 iterations: <a href="https://drive.google.com/file/d/1IFqzcjkCxmlQjJwGdfwp7ACZVsP9Dd-Q/view?usp=sharing" rel="nofollow noreferrer">https://drive.google.com/file/d/1IFqzcjkCxmlQjJwGdfwp7ACZVsP9Dd-Q/view?usp=sharing</a></p> <p>When I use the source term b and keep the relaxation parameter w=1, i.e.</p> <pre><code>T[1:-1,1:-1] = (1-w)*Tn[1:-1, 1:-1] + w*0.25*(T[1:-1, 2:] + T[1:-1, :-2] + T[2:, 1:-1] + T[:-2, 1:-1] - dx**2*b[1:-1, 1:-1]) </code></pre> <p>I get this plot, which takes a long time to converge: <a href="https://drive.google.com/file/d/1diOWXXY1lnP9IEgJsz7Ih7mIBAhxdtBk/view?usp=sharing" rel="nofollow noreferrer">https://drive.google.com/file/d/1diOWXXY1lnP9IEgJsz7Ih7mIBAhxdtBk/view?usp=sharing</a></p> <p>Finally, when I introduce the relaxation parameter w as any value in (1.1, 1.2, 1.3....1,9) the values and the errors shoot up really high and the temperature matrix values shows 'inf'.</p> <p>It seems like the problem occurs when I use the relaxation parameter, but I cannot figure out why it's happening.</p> https://scicomp.stackexchange.com/q/35427 0 Integrating Matrix Elements TypeError: f() takes 1 positional argument but 3 were given New2Python https://scicomp.stackexchange.com/users/36389 2020-06-25T03:03:20Z 2020-11-22T04:06:40Z <p>I'm working on a linear variational problem for a general PIB and I keep encountering the same problem, and I know its a rather simple solution. Any suggestions?</p> <pre><code>import numpy as np import scipy.integrate as sp #Parameters alpha = 0.1 L = 10.0 N = 5 x = np.linspace(0,L,100) #Matrix elements linear combination components def f(x,n,L): return(np.sqrt(2./L)*np.sin(n*np.pi*x/L)) def f1(x,m,L): return( np.sqrt(2./L)*np.sin(m*np.pi*x/L)) #Linear Potential def V(x): return( alpha*x) #Matrix element functions to integrate def Int_1(x,n,m,L): return(f1(x,m,L)*f(x,n,L)+f1(x,m,L)*V(x)*f(x,n,L)) def Int_2(x,m,n,L): return(f1(x,m,L)*V(x)*f(x,n,L)) # Trying to integrate the matrix components def c(m,n,L): return(sp.quad(Int_1,0,L, args = (m,n,L), limit = 100)) def c1(m,n,L): return(sp.quad(Int_2,0,L, args = (m,n,L), limit = 100)) # Generating the matrix def cal_Hmn(m,n): if m == n: c(m,n,L) elif(m+n)%2 ==1: c1(m,n,L) else: return(0) #Filling the Matrix H = np.zeros((N,N), float) for i in range(N): for j in range(N): H[i,j] = cal_Hmn(i+1, j+1) </code></pre> <p>I am positive the issue arises in the <span class="math-container">$&lt;{\psi}|{\hat{H}_{mn}}|{\psi}&gt;$</span> integrals, Thanks!</p> https://scicomp.stackexchange.com/q/35202 1 How to write a code of 2D ADI method in matlab? I love math https://scicomp.stackexchange.com/users/36100 2020-05-22T14:03:41Z 2020-11-21T17:06:00Z <p>I tried to write a code for the <a href="https://en.wikipedia.org/wiki/Alternating_direction_implicit_method" rel="nofollow noreferrer">alternating direction implicit (ADI) method</a> in 2D, but I got stuck.</p> <p>My equation is:</p> <p><span class="math-container">$$\frac{\partial U(t,x,y)}{\partial t} = 2\Delta U(t,x,y) -10(\frac{\partial U(t,x,y)}{\partial x} + \frac{\partial U(t,x,y)}{\partial x}) + [(x^2 + y^2 + 4t(5x + 5y -2)] \; \forall \;(x,y) \in (0,1)\times(0,1),\; t\in(0,T]$$</span></p> <p>Initial and boundary conditions are:</p> <p><span class="math-container">\begin{align} &amp;U(0,x,y) = 0 \;\forall\;(x,y)\in[0,1]\times[0,1]\\ &amp;U(t,x,y) =t (x^2 + y^2) \; \forall \; (x,y) \in \partial[0,1]\times[0,1] \end{align}</span></p> <p>I solved this equation (I used ADI method) and I found:</p> <p><strong>1.step</strong> <span class="math-container">$(-15)u_{i-1,j}^{k+1/2} + (26)u_{i,j}^{k+1/2} + (-10)u_{i+1,j}^{k+1/2} = (15)u_{i,-1j}^{k} + (-24)u_{i,j}^{k} + (10)u_{i,j+1}^{k} + [\frac{x_i^2}{2} + \frac{y_j^2}{2} + 2t(5x + 5y -2)]$</span> <strong>from Dirichlet condition:</strong> <span class="math-container">$i=1 \qquad \qquad u_{0,j}^{k+1/2}= t_{k+1/2}(0^2 + y_j^2) \\ i=n-1 \; \qquad u_{n,j}^{k+1/2}=t_{k+1/2}(1^2+y_j^2)$</span></p> <p><strong>2.step</strong> <span class="math-container">$(-15)u_{i-1,j}^{k+1} + (26)u_{i,j}^{k+1} + (-10)u_{i+1,j}^{k+1} = (15)u_{i,-1j}^{k+1/2} + (-24)u_{i,j}^{k+1/2} + (10)u_{i,j+1}^{k+1/2} + [\frac{x_i^2}{2} + \frac{y_j^2}{2} + 2t(5x + 5y -2)]$</span> <strong>from Dirichlet condition:</strong> <span class="math-container">$j=1 \qquad \qquad u_{i,0}^{k+1}= t_{k+1}(x_i^2 + 0^2) \\ j=n-1 \; \qquad u_{i,n}^{k+1}=t_{k+1}(x_i^2+1^2)$</span></p> <p>Although I found both steps, but I do not know, hot to write it in Matlab. Could somebody help me, please?</p> <p>My attempt is:</p> <pre><code>n = 10; h = 1/n; tau = 1/10; mi = tau/h/h; g0x=@(t,y)0.*x; g1x=@(t,y)0.*x; g0y=@(t,y)0.*y; g1y=@(t,y)0.*y; u0=@(x,y)0,*+0,*y; x = 0 : h : 1; y = 0 : h : 1; t = 0 : tau:1; [X,Y]=meshgrid(x,y) u = zeros(J+1,Nt); for i = 1:n+1 u(i,n)=(-15)*u(i-1,n) + 26*u(i,n) + (-10)*u(i+1,n)+ (-15)*u(i,n-1) + 24*u(i,n)+(-10)u(i,n+1) end u(0,n)=t*u(y*y) u(n,n)=t*u(1+y*y) for j = 1:n+1 u(n,j)=(-15)*u(n,j-1) + 26*u(n,j) + (-10)*u(n,j+1)+ (-15)*u(n,j-1) + 24*u(n,j)+(-10)u(n,j+1) end u(0,n)=t*u(y*y) u(n,n)=t*u(1+y*y) </code></pre> https://scicomp.stackexchange.com/q/33619 1 Robin Boundary Condition with Implicit Upwind - Finite Difference Method for 2D Convection-Diffusion Equation Tan Phan https://scicomp.stackexchange.com/users/33088 2019-10-20T03:43:33Z 2020-11-21T09:08:30Z <p>I am trying to solve a problem with 2D Convection-Diffusion equation with U = Concentration (<span class="math-container">$mg/m^{2}$</span>) using Implicit Upwind Finite Difference Method like this</p> <p><span class="math-container">$$\frac{\partial U}{\partial t} + v_{y}\frac{\partial U}{\partial y} = D(\frac{\partial^{2} U}{\partial x^{2}}+ \frac{\partial^{2} U}{\partial y^{2}})$$</span></p> <p>Set <span class="math-container">$\sigma_D = D\frac{\Delta t}{\Delta x^{2}} = D\frac{\Delta t}{\Delta y^{2}}$</span> and <span class="math-container">$\sigma_{C} = \frac{v_{y}\Delta t}{\Delta y}$</span></p> <p>I have dicretized the equation above using Implicit Upwind scheme as I have vy = -0.01 m/s going down towards the ground:</p> <p><span class="math-container">$$U_{i,j}^{m+1}(1+4\sigma_D+\sigma_{C})-\sigma_{D}(U_{i+1,j}^{m+1}+U_{i-1,j}^{m+1} + U_{i,j+1}^{m+1})-(\sigma_{C}+\sigma_{D})U_{i,j-1}^{m+1} = U_{i,j}^{m}$$</span> with domain [x,y] = [10km,4km]. Dirichlet BCs are at the top, left and right of the domain and Robin BC is at the bottom (modelled as ground):</p> <p><span class="math-container">$$D\frac{\partial U}{\partial y} - v_{y}U = 0$$</span></p> <p>Set <span class="math-container">$\sigma_{R} = \frac{v_{y}\Delta y}{D}$</span></p> <p>Dicretized Robin BC using downwind method (since the wind is going down):</p> <p><span class="math-container">$$v_{y}U_{i,1}^{m+1}-D(\frac{U_{i,1}^{m+1} - U_{i,0}^{m+1}}{\Delta y})$$</span></p> <p>Rearrange the equation above gives:</p> <p><span class="math-container">$$U_{i,1}^{m+1}(\sigma_{R}-1)+U_{i,0}^{m+1}=0$$</span></p> <p>Initial Conditon: 1000kg of chemical is dropped in the middle of domain [5km,2km] ( U = Concentration [mg/m2] ):</p> <p><span class="math-container">$$U_{init} = \frac{1000E6}{\Delta x \Delta y}$$</span></p> <p>This is my sample code to run with this method:</p> <pre class="lang-py prettyprint-override"><code>Lx = 10E3 Ly = 4E3 dx, dy = 20, 20 nx = int(Lx/dx + 1) ny = int(Ly/dy + 1) D = 0.5 x = np.linspace(0.0, Lx, nx) y = np.linspace(0.0, Ly, ny) mid_x = int(nx/2) mid_y = int(ny/2) mass = 1000E6 #mg vy = -0.01 sigma_D = 1 dt = sigma_D * min(dx, dy)**2 / D sigma_c = vy*dt/dy day = 4 time = 3600*24*day nt = int(time/dt) def IC(nx, ny, mass, dx, dy, mid_x, mid_y): Dis = np.zeros((ny, nx)) Con = mass / (dy*dx) Dis[mid_y, mid_x] = Con return Dis, Con U0, C = IC(nx, ny, mass, dx, dy, mid_x, mid_y) def Implicit(U0, nt, dx, dy, D, vx, vy, frn, dt, nx, ny): sigma_D = D*dt/(dx)**2 sigma_c = vy*dt/dy sigma_R = vy*dy/D AA = csr_matrix((nx*ny, nx*ny)).tolil(copy = True) #Boundary Condition #Dirichlet for j in range(ny): for i in range(nx): AA[i + j*nx, i + j*nx] = 1 #Inner nodes for j in range(1, ny - 1): for i in range(1, nx - 1): #n + 1 side AA[i + j*nx, i + j*nx] = 1 + 4*sigma_D + sigma_c AA[i + j*nx, i+1 + j*nx] = -sigma_D AA[i + j*nx, i-1 + j*nx] = -sigma_D AA[i + j*nx, i + (j+1)*nx] = -sigma_D AA[i + j*nx, i + (j-1)*nx] = -(sigma_D +sigma_c) #Robin BC for j in range(0, 1): for i in range(nx): AA[i + j*nx, i + j*nx] = -1 AA[i + j*nx, i + (j+1)*nx] = -sigma_R + 1 Matrix_1 = AA.tocsr() for n in tqdm(range(nt)): U0 = 0 U0 = csr_matrix(U0.reshape(ny*nx, 1)) U1 = scipy.sparse.linalg.bicgstab(Matrix_1, U0.todense()) U0 = U1.copy().reshape((ny, nx)) return U0 </code></pre> <p>This is what I got from the code after 4 days (3600*24*4 seconds):<br> Concentration Plot in 3D <a href="https://i.stack.imgur.com/cai3f.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/cai3f.png" alt="Concentration Plot"></a></p> <p>Total concentration over time<br> <a href="https://i.stack.imgur.com/kPuHL.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/kPuHL.png" alt="Total Concentration"></a></p> <p>I am not sure why my total concentration keeps losing over time. Theoretically, it should stay the same over time since when it hits the ground, it is stuck there. </p> <p>I believe my boundary condition for Robin BC might be wrong but I am not sure where I went wrong. Any suggestions?</p> https://scicomp.stackexchange.com/q/30406 3 GPGPU/FPGA programming for Combinatorial Analysis dogman_1234 https://scicomp.stackexchange.com/users/29034 2018-10-22T23:07:07Z 2020-11-24T09:01:30Z <p>Recently, I have taken an interest in performing combinatorial analysis for the game of 21 (blackjack) and attempted to use my AMD APU to try and thread the program via the 4 cores on the chip. Generating the proper strategy and expectations has proven to be a success, as well as taking very few seconds (around 12) to compute the optimal strategy for any given rule(s) and deck composition. However, when it comes to computing other deck compositions, it would take 12.63 years to compute all deck states...for a single deck. So, another approach would be to simulate a select limit of random deck subsets via a Monte-Carlo method. This will work! However, the 12 seconds per analysis is a <em>huge</em> bottleneck! Randomly selecting just 10E6 subsets would take around 139 days. Not feasible with my current CPU, even with threading.</p> <p><em>A primer on how I am analysing the strategy chart(s):</em></p> <p>1.) I first start by computing what the overall expectation for standing is by taking a set of n_cards for the player and enumerating all possible dealer hand sets. I compare the outcome of each specific dealer total against the current hand total of the player hand.</p> <p>2.) I then compute the hitting expectation by taking a player hand set {23} and finding matching hand sets for x:2->A, so that we can compute the weighed expectation of hitting to {232}, {233}, {234}, ...,{23A}. We do this by computing all hard 21's first, then hard 20's...down to all hard 12's; then do the same by computing hitting all soft 21's to soft 13's; then compute all hard 11's to hard 5; and finally, we compute hitting all pair cards ({22}, {33}, ..., {TT}, {AA}.)</p> <p>3.) Doubling is the same as hitting, except we only draw one card and stand for any 2_card player hand. We return the expectation of drawing to a specific 3_card hand and multiply the expectation by 2 times the probability of that drawn card.</p> <p>4.) For splitting, we will take steps 1.) to 3.), we will repeat 1->3 up to 11 times. Why? We are computing the expectation of splitting to a new hand by removing the rank for which we are splitting. (For example: We want to split {22}. We will compute standing, hitting, doubling for a given deck subset, minus one of the ranks we are splitting. If we originally split {22} for a deck subset of {4 4 4 4 4 4 4 4 16 4: 2-A} or {3 1 4 3 4 4 3 1 14 3: 2-A}, our new subsets would be {3 4 4 4 4 4 4 4 16 4: 2-A} or {2 1 4 3 4 4 3 1 14 3: 2-A}.)</p> <p>5.) After all of this, we should get back out optimal strategy with respect to the computed expectation(s) for each possible decision for each player hand.</p> <p>Now, I have been looking into seeing if there is a quicker way of analysing the strategy using either GPGPU or FPGA/ASIC. I first looked at GPU programming as a possible route, however, I am not sure if it is feasible as there are 3072 unique player hand states, meaning there will be 3072 unique data structures to compute for standing, hitting, doubling, and splitting.</p> <p>With all this said, I was hoping there was a way to perform some form of simultaneous processing, where the state of each expectation changes dependent on the change in the rules and/or the deck state. That is, rather than recompute each expectation for every change in rules/deck_state, there would be some dependence on the output for each data point. For example: We have three data points named A, B, and C. We also have state data {1 2 3}. For each data point, there is a state equation that we possess:</p> <p>A = data</p> <p>B = A + data</p> <p>C = B +data</p> <p>Under serial processing, we would evaluate A first for 1, B for 2 plus A, and C for 3 plus B. If we change data to be {2 2 3}, then A would be 2, B would be 4, and C would be 7. Under serial processing, this would take several cycles per data point, for each data point. Is there a method under FPGA where the change in state is instantaneous? That is if we change any point of data, that A, B, and C will change after data changes, simultaneously?</p> <p>If not, what benefits does FPGA/GPGPU processing offer to a programmer who wants to accurately and quickly compute large volumes of floating point data in a near instant? Basically, is there a way to rapidly compute the optimal strategy of a blackjack game that is faster than multi-threaded CPU processing that is not 12 seconds long?</p> https://scicomp.stackexchange.com/q/28388 10 Resources on mesh generation for finite element methods philm https://scicomp.stackexchange.com/users/18230 2017-12-05T15:56:30Z 2020-11-26T00:28:05Z <p>I know that this is not really apart of the rules as this is a recommendation question and these don't really have an answer per say. But, like this forum posting: <a href="https://stackoverflow.com/questions/388242/the-definitive-c-book-guide-and-list">https://stackoverflow.com/questions/388242/the-definitive-c-book-guide-and-list</a>. To be honest, I am not really sure where else to ask this question.</p> <p>I would like to start a book guide list on Finite Element and meshing. The thing is, I am not doing this as university research nor as job research. I am interested in this topic and have been for some years. I would like to learn this on my own time. I have been working on building my own simulator for non-linear electromagnetic simulations in 2-D and I am currently using gmsh as the mesher. I am currently working on integrating gmsh into my source. Progress has been positive. I am able to create a mesh with gmsh where the source code is directly integrated into my project. I would like to modify the current workflow for meshing which means that I need to code to work around GMSH limitations. THis also means, I need to have an understanding of creating a numerical grid for simulations so that I can navigate the source of GMSH and have a better understanding of what is going on. </p> <p>I am finding that I have a lack on this aspect (numerical grid generation) since the source is referring to alot of terms that I am not familiar with. I apologize if this is breaking any forum rules. </p> <p>But, I was wondering if anyone can point me to any resources about numerical grid generation? What would be some good references to use for introductory, beginner, intermediate, and advanced? Right now, I feel like I am at beginner. I currently have the resource "Numerical Grid Generation: Foundations and Applications by Joe F. Thompson". Would this resource be good as a starting point?</p> <p>(I will definitely be going through the deal.II tutorials/manual as the documentation is very detailed. I like that it doesn't quite overload the reader with technical details but explains it with simplicity but to the point)</p> <p>As a side note, I think that there should be a section for videos. Sometime, those can be helpful.</p> https://scicomp.stackexchange.com/q/23940 5 Optimization on the manifold of stochastic matrices Thoth https://scicomp.stackexchange.com/users/6849 2016-05-09T05:57:42Z 2020-11-21T12:51:40Z <p>So I have an optimization problem of the form $$\text{maximize}\hspace{3mm}f(A):{\bf R}^{K\times K}\rightarrow{\bf R}$$ $$\text{subject to}\hspace{19mm}A^T{\bf 1}=\bf{1}$$ $$\hspace{33mm}A\geq 0$$ where the constraint set is thus the manifold of left stochastic matrices (each column sums to 1). The objective function $f$ is not concave and thus my goal is to very quickly find a local maxima (speed is important).</p> <p>The function $f$ is a nasty piece of work, nevertheless the gradient $\nabla f$ is analytically computable. However even without the constraint, setting the gradient equal to zero and solving for the extrema is intractable, and thus I can't imagine that back solving for Lagrange multipliers is doable.</p> <p>The method I've tried is manifold gradient ascent using some code found <a href="http://www.manopt.org/reference/manopt/manifolds/multinomial/multinomialfactory.html" rel="nofollow noreferrer">here</a> from the Manopt package. They refer to the manifold of left stochastic matrices as the <em>multinomial manifold</em>, and the code provides methods for basically projecting the gradient onto the tangent space, and then retracting this new gradient onto the manifold itself.</p> <p>Unfortunately as far as retractions go it's fairly primitive I think, and it involves entry-wise dividing the projected gradient by the current point and then exponentiating it, which causes overflow unless you make your step-sizes quite small (however I don't fully understand the intricacies of manifold optimization so I could be wrong about this).</p> <p>I could compute the Hessian numerically and try to run a 2nd order optimization on the multinomial manifold, but I'd have to do some reading to figure out exactly how to do this. Basically I need something that is <em>fast</em>, speed is more important than accuracy, as I have to do this thousands of times and it's only 1 step in a larger coordinate ascent algorithm. </p> <p>The gradient $\nabla f$ is complicated and thus each evaluation of it is quite costly, however the one saving grace is that the dimension of $A$ is relatively small, almost certainly $&lt;100$.</p> <p>What would be the best approach here?</p> <p><strong>Update</strong></p> <p>As requested, here is the function and its gradient: \begin{align} f(A )&amp;= \Big[-\psi\big(\sum_{i=1}^KC_{ij}\big)+\sum_{i=1}^KA_{ij}(\psi(C_{ij})-\log C_{ij})\Big]_{j=1:K} \cdot \big(\sum_{n=2}^NA^{n-2}\big)\cdot {\bf x}\\ &amp;+ \sum_{n=1}^N[\log B_{in}]_{i=1:K}^T\cdot A^{n-1}\cdot{\bf x} \end{align}</p> <p>\begin{align} \nabla f(A) &amp;= [\psi(C_{ij}) - \log A_{ij} -1]_{i=1:K,\hspace{1mm}j=1:K}\cdot \text{diagm}\Big(\big(\sum_{n=2}^NA^{n-2}\big)\cdot{\bf x}\Big)\\ &amp;+ \sum_{n=3}^{N}\sum_{r=0}^{n-3}\Big((A^r)^T\cdot\Big[-\psi\big(\sum_{i=1}^KC_{ij}\big)+\sum_{i=1}^KA_{ij}(\psi(C_{ij})-\log C_{ij})\Big]_{j=1:K}^T \cdot {\bf x}^T\cdot (A^{n-3-r})^T\Big)\\ &amp;+ \sum_{n=2}^N\sum_{r=0}^{n-2}(A^r)^T\cdot[\log B_{in}]_{i=1:K}\cdot{\bf x}^T\cdot(A^{n-2-r})^T \end{align} Where $\psi$ is the digamma function, $A$ and $C$ are $K\times K$ matrices, $B$ a $K\times N$ matrix and ${\bf x}$ is a length $K$ vector, all real-valued. The function $\text{diagm}$ converts a vector to a diagonal matrix. Also I checked my derivation of the gradient numerically, so you can be sure its correct.</p> <p>I'm writing this in Julia, and I tried using the NLopt package both using a gradient-based method and a derivative free method, but they were even slower than my original manifold approach.</p>