Hot answers tagged

8

Since your figure is a closed loop, its parametric curves $x(t)$ and $y(t)$ must be periodic functions. This suggests one way to generate such figures, by constructing random smooth periodic functions $x(t)$ and $y(t)$ via summation of sinusoids/harmonics with randomized amplitude and phase. Unfortunately, it would be difficult to guarantee such a figure ...


7

This is my comment expanded into an answer. MATLAB is a column-major programming language; it is not very hard to find it in the documentation once you know what to look for: https://www.mathworks.com/help/matlab/matlab_external/matlab-data.html#f22019 The layout of data structures (most commonly row-major or column-major) used in a programming language is ...


7

Julia's DifferentialEquations.jl has a lot of tooling for automatically deriving (sparse) matrices. For more information, see the JuliaCon 2020 video on Auto-Optimization and Parallelism in DifferentialEquations.jl. Combined with the orders of magnitude acceleration commonly seen over the MATLAB solvers, this might be a good option for you and is a quick ...


6

Q1: No. Here's a counter-example: >> A = eye(4)*1e-300 A = 1.0e-300 * 1.0000 0 0 0 0 1.0000 0 0 0 0 1.0000 0 0 0 0 1.0000 >> rank(A) ans = 4 >> >> rank([A, ones(4, 1)]) ans = 1 >> so, if the added column ...


6

As Federico has mentioned, you probably don't want to deprive yourself of the learning experience. I'll just give you a small nudge in the right direction. You will never be able to store $A$. You also won't be able to store $(A+uv^T)^{-1}$. However, you don't really need to. You can easily write down a formula for each of the entries in $A$. Instead of ...


6

I got your point: you have the ODE $y''=y'+x$, and we can see that $y(x)=-\frac{x^2}{2} - x$ solves the problem. As you want to integrate it numerically, you want to set two initial conditions, so you imposed $y(0)=0$ (and that's fine) and also another one $y(h)$, which is wrong because it's not the required information you need to solve an ODE, mainly ...


5

You are starting from a uniform temperature and you have insulated boundary conditions; so there is no heat conduction occurring. Likewise, your initial mole fraction is also constant in $x$ so that the heat input is uniform along the length. So the fact that the temperature and mole fraction don't change as a function of $x$ is exactly what you should ...


5

Pseudoinverse can be computed using the SVD $A = USV^\top$ by: $$ A^+ = V\Sigma^+ U^\top $$ where $\Sigma^+$ is formed from $\Sigma$ by taking the reciprocal of all the non-zero elements. With that in mind, you could use MATLAB's svds function as follows: [U,S,V] = svds(A,k); Ainv = V*diag(1./diag(S))*U'; Here k refers to the rank and svds computes only a ...


4

The reason that this particular mesh does not give the correct, uniform displacement solution to this problem is that it is "non-conforming." Specifically, at the intersection of the two cubes in the model, the two element edges that cross that face don't align with each other but instead cross each other. A typical face in a non-conforming mesh ...


4

I see that you are doing some redundant computation, for example y.^2+z.^2 is computed 4 times, with y.^2 and z.^2 are computed 9-10 times. You can define a set of variables y2=y.^2, z2=y.^2, y2pz2=y2+z2; and push some of the computation cost to memory -given that you have enough memory. That would save you a good amount of time. MATLAB is a column-major ...


4

This is easy to formulate in CVX, under MATLAB. A CVXPY solution, under Python, is similar. CVX code: cvx_begin sdp variable X(n,n) hermitian semidefinite minimize(norm_nuc(X-A)) X <= B cvx_end or alternatively cvx_begin variable X(n,n) hermitian semidefinite minimize(norm_nuc(X-A)) B - X == semidefinite(n) cvx_end Edit 2: CVX is very fussy about ...


4

You simply have a bug in your code. The flux is $\frac{1}{2} u^2$ and not $\frac{1}{4} u^2$.


4

Some things I can think of: use sparse matrices for Matrix1 and Matrix2 to speed up the computations of dZand dY use larger integration tolerances reltoland abstol, especially if your are searching for steady-state solutions and/or do not need a precise resolution of the transient dynamics of your system. you are solving with ode15s, which is an implicit ...


4

T = 1.0; f = @(x,t) (0 <= t & t <= T/2)*20*exp(-4*x.^2) + 0.0; x = linspace(0,1,100); plot(x, f(x,T/4), x, f(x,T)) legend('t=T/4','t=T')


4

The result you see is due to the long integration interval and the accumulation of truncation errors. Each step with its truncation error will switch to a slightly different exact solution. Now the exact solutions are in general $$ y(x)=-\frac12x^2-x+C+De^x $$ Due to the first order approximation of the first derivative, the coefficients $C$, $D$ will grow ...


4

It's not overcomplicated: That's how $x^y$ is actually defined -- via $$ x^y = \exp(y \log x). $$ The point is that for non-integer $y$, it's not at all obvious what $x^y$ actually means: It is not just a product of $x$ with itself, repeated $y$ times. The definition of $x^y$ then is done by pulling the expression back to functions we have previously ...


4

While $a^{b}$ is usually computed via $\exp(b \log (a))$, in real-life high-quality implementations of functions like pow() the logarithm is computed with extra precision to counter the error magnification effect of $\exp$, as demonstrated in this answer for example. This suggests that we should favor computation via a built-in general exponentiation ...


4

It looks like you have a advection diffusion PDE discretized with finite differences. This gives an ODE of the form $$ y' = f(y) = A y + D y, $$ where $A$ is the discretized advection operator and $D$ is the discretized diffusion operator. In your case, it seems you have zero Dirchlet and Neumann boundary conditions which can be encoded in the $A$ and $D$ ...


4

Before even discussing any vectorization, there are gross inefficiencies in the code you provided. Just by employing basic techniques, I got 21-fold speedup: n = 1000000; y = zeros(1,n); tProbs = [0.1 0.5; 0.9 0.5]; iProbs = sum(tProbs^100000,2); iProbsSqr = iProbs.^2; V0 = [1;1]; for i=1:n V0 = tProbs*V0; y(i) = sqrt(dot(V0.^2,iProbsSqr)); end If the ...


3

The notation $$\frac{\partial U}{\partial \eta}$$ means usually $$\eta \cdot \nabla U$$. This is correct even if the domain is the interval $[a,b]$. The normal vector on the interval $[a,b]$ @a is $\eta=-1$ and @b $\eta= 1$ both pointing outwards of the domain. Hence in 1D $\frac{\partial U}{\partial \eta}$ means $$\eta\cdot\nabla U=\eta \frac{dU}{dx}$$.


3

As far as I understand the fist part of your question is regarding the options of the ode15s, and which ones you can use to make the solver more efficient/accurate. Providing the Jacobian to the solver, especially if its reasonably easy to calculate, is a good idea to make things faster and more accurate. Using sparse systems for small systems of equations ...


3

Suppose you want to minimize $$\Phi(x)=\frac{1}{2}||Ax-b||^2$$ The gradient is $$\frac{\partial \Phi}{\partial x} = A^T(Ax-b)$$ The step size to guarantee convergence is $$\alpha=||A^TA||^{-1}$$ Why? The direct solution to the problem is: $$x_{opt}=(A^TA)^{-1}A^Tb$$ This can be achieved iteratively if we look at the update on the estimate $x_k$. Suppose we ...


3

You lose convergence order, or in the worst case convergence altogether. You can try this out: Take $$ f(x) = \begin{cases} 0 & \text{if $x<0$} \\ x^2 & \text{if $x\ge 0$}.\end{cases} $$ The function is differentiable, but not twice differentiable at $x=0$. It's exact derivative is $f'(0)=0$. Now compute the (second-order) symmetric finite ...


3

This problem is known as joint diagonalization, and it has two variants: orthogonal, in which the basis vectors are orthonormal, and non-orthogonal, which is harder to solve, but which may be more appropriate to your application. The simplest method I know of seeks a unitary matrix $U$ that minimizes the sum of squares of the off-diagonal elements of $U^HAU$...


3

The effect is due to quantization noise and/or aliasing, you are not computing with the true minima of the radius along the orbit, but with the closest sample point. This means that the sample point for the half-step integration can lie in-between. Close to the minimum this results in an $O(e·dt^2)$ ($e$ the eccentricity) difference between the computed ...


3

This problem is too small to actually be sparse. Sparse handling has a big overhead because the indexing is not "direct", i.e. you don't necessarily know where the next value will be without branch checking. So you need it to be "sparse enough" that the O(n^3) dense LU-factorization cost shrinking to the purely non-zero terms overcomes ...


3

Since you already know Matlab and Mathematica, I would go with SciPy Lecture notes. It is a course on using Python for Scientific purposes. Gael Varoquaux, Valentin Haenel, Pierre de Buyl, Gert-Ludwig Ingold, Emmanuelle Gouillart, Michael Hartmann, … João Felipe Santos. (2017, October 4). scipy-lectures/scipy-lecture-notes: Release 2017.1 (Version 2017.1). ...


3

Probably not what OP was waiting for, but I think it could be pretty instructive and useful. FEM codes use a much different approach to build the so-called stiffness matrix. In practice, they loop over elements and compute for each element small matrices (in your case if you use linear elements 3by3) which are distributed to the right entries of the global ...


3

(This answer is valid for both MATLAB and Octave, even though I mainly refer to MATLAB) There are two beasts to slay; but let's first understand the underlying data structure. MATLAB and Octave store sparse matrices in the coordinate (COO) format, i.e. a sparse matrix S is a collection of three arrays of the length equal to the number of nonzero entries of S:...


2

The elemental stiffness matrix must be always singular because while deriving it we do not impose any constraints or boundary conditions. Thus inverting the stiffness matrix to solve for displacements/position-vectors/degrees-of-freedom should yield indeterminate results.


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