numerical-analysis

[dusk meadow]

Ah I see. Technically what im doing is changing the basis of a vector. I'm basically solving for z in Gz = x, where G(x) is a pos-def matrix and a function of x, and I know x. Though, I believe this implies I must invert G, no?

It is not sparse. I haven't heard of a dense matrix before though.

[plucky kayak]

no, all you need is to pass G and x to a linear system solver

you would be already doing one more unnecessary matrix-vector multiplication by computing G^{-1} * x, which could've been avoided by solving for z in Gz=x

plus spend more memory on storing the inverted matrix

all of this could've been avoided by solving Gz=x and only storing z

[cosmic karma]

as the other person said, you have a plethora of choices to solve for z in Gz=x because G is positive definite. Inverting matrices are much more expensive than solving a linear system. You have Gauss-Seidel, LU factorization, Richardson methods, etc.

[dusk meadow]

Ill look into them. Thank you.

[fathom rain]

\

[silk drum]

If you want to process a high number of these efficiently on a computer, the LAPACK and PLASMA libraries are worth looking into for implementations of the algorithms you've been pointed to. If you are just looking to solve these one-off, MATLAB or Octave might be more user-friendly. Many of the dense algorithms factor the matrix G and allow G to be inverted without the trouble of actually finding the explicit inverse of G and applying the explicit inverse of G.

[dusk meadow]

I'm using python, but it seems that I dont have to invert G, but simply solve a linear equation which could speed things up. So, solving Gx = z, where G is a non-sparse positive definite with dimension 100x100, what algorithm do you think is the best choice?

[silk drum]

DPOSV seems reasonable. Some Python library should provide or directly call down to LAPACK.

https://netlib.org/lapack/explore-html-3.6.1/dc/de9/group__double_p_osolve_ga9ce56acceb70eb6484a768eaa841f70d.html

[uneven thicket]

Late to the party but numpy.linalg.solve will do you fine

numpy is backended by LAPACK

[dusk meadow]

@silk drum @uneven thicket Thank you, I'll check these out 🙂

[heady mesa]

LU decomposition is generally the best if you want the theory behind them.

Mpmath is a good python module for Numerical Maths enthusiasts, it allows you to use basically infinite precision and so you won't need to worry about matrix conditioning when inverting or solving linear equations.

[silent moat]

Today i had my exam at uni

Passed numerical analysis with 10

Thanks guys

[heady mesa]

👑

[jaunty saddle]

copy and pasted from dynamical systems

any help is appreciated

[worn creek]

Hello, could you please tell me what this symbol is called and how it is used

https://thumbs.dreamstime.com/z/contour-integral-sign-contour-integral-sign-white-background-189550343.jpg?w=768

[plucky kayak]

contour/line integral

[errant wagon]

You might find advanced engineering mathematics by Dennis Zill interesting

[cinder sandal]

I am having trouble implementing Simpson's rule. Can anyone help?

 simpsons <- function(f, partitions) {
  # Calculate the width of each subinterval (assumes equal width)
  a <- partitions[[1]][[1]]
  b <- partitions[[length(partitions)]][[2]]
  n <- length(partitions)
  h <- (b - a) / n
  
  # Adding initial endpoints f(x_0) and f(x_n)
  total <- f(a) + f(b)
  
  for (i in 2:(length(partitions) - 1)) {
    x_mid <- (partitions[[i]][[1]] + partitions[[i]][[2]]) / 2
    if (i %% 2 == 0) {
      # Even index: 2 * f(x_i)
      total <- total + 2 * f(x_mid)
    } else {
      # Odd index: 4 * f(x_i)
      total <- total + 4 * f(x_mid)
    }
  }
  
  # Multiply by h/3 and return the result
  total * h / 3
}

[fathom rain]

is that R? they do not have slicing/views? this is an example in python def simpson(f,a,b,N): h=(b-a)/(N-1) xi=np.linspace(a,b,N) fxi=f(xi) return 2*h*(sum(fxi[:-2:2]+4*fxi[1:-1:2]+fxi[2::2]))/6

[cinder sandal]

Yeah it’s R. Here is an update to my code.

latest code. I'll figure this out later.

simpsons <- function(f, partitions) {
  endpoints <- sort(unique(unlist(partitions)))
  
  # Calculate the width of each subinterval (assumes equal width)
  a <- min(endpoints)
  b <- max(endpoints)
  n <- length(endpoints) - 1
  h <- (b - a) / n
  
  # Adding initial endpoints f(x_0) and f(x_n)
  total <- f(a) + f(b)
  
  for (i in 2:n) {
    if (endpoints[i] != a & endpoints[i] != b) {
      if (i %% 2 == 0) {
        total <- total + 2 * f(endpoints[i])
      } else {
        total <- total + 4 * f(endpoints[i])
      }
    }
  }
  
  total * h/3
}

[cinder sandal]

The question asked me to use three randomly generated numbers in [0, 1].

[cinder sandal]

Python slicing is super nice.

[vapid plume]

Hello. Is anyone here familiar with discrete-time waveform relaxation for parabolic PDEs? Trying to implement this method in MATLAB using red-black Gauss-Seidel but my results aren't correct. After discretizing in space, we use the trapezoidal rule to discretize in time, so we just get the Crank-Nicolson scheme. As far as I know, this is supposed to be different from the standard time-marching schemes is by first looping over iterates, then over time-steps. If this is correct, then I just have questions about the initial guesses for the iterates.

Currently, for the first iteration for the first time-step, I am using the initial condition as the initial guess and then for the remaining time-steps, I am using the previous time-step iterations as the initial guesses. For the remaining iterations, I am using the previous time-step iteration for the first time-step and then for the remaining time steps, I have tried two approaches: Use the previous time-step iterations as the initial guesses or the previous corresponding time-step iterations as the initial guesses.

If this is how it is supposed to be, then there is something else I am doing wrong. Convergence is quick with these approaches and I am expecting slower convergence.

[raw sand]

Hi, I'm hoping I'm in the correct channel. I'm writing a video game that resembles a simulation of particles(Discrete Element Method) interacting with themselves and some physical/rigid boundaries. I'm using a physics engine that can model gravity/collisions pretty well, but it doesn't do anything like fluid dynamics or electromagnetics/electrostatic modeling. I'm able to apply external forces to the particles(angular and linear). So I was thinking I could model charged particle interactions by overlaying a mesh over the volume containing the particles and applying the particle charges as (I think) boundary conditions.

However, I'm not able to find examples of multiple charges. Scikit-fem has an example that models one charge using greens theorem: https://github.com/kinnala/scikit-fem/blob/9.0.1/docs/examples/ex38.py

The problem I'm having with sklearn-fem is the models are too abstracted and the effort to understand the math from reading the code is going to take me a long time.

I'm wanting to adapt this to use fenris (Rust based FEM), multiple point charges and a generalized mesh. Fenris has a poisson example, but it's the standard non-point source textbook example: https://github.com/InteractiveComputerGraphics/fenris/blob/master/examples/poisson2d.rs

My gut instinct is that the setup for this would be really easy if I knew what the heck I was doing 😄

[fathom rain]

maybe check out FMM https://en.wikipedia.org/wiki/Fast_multipole_method

[raw sand]

That definitely looks useful

I think I'm still going to have a difficult time implementing this for my game 😄

I took multiple differential equations and linear algebra classes 20 years ago, but everything I'm reading on FEM or now in FMM is making me feel really dumb

[fathom rain]

using fem for a game seems like overkill 🙂

i just has to look good, not be accurate

or maybe you want it accurate

[raw sand]

It could be. I decided to spend 1 day looking into it, since I thought applying forces to particles might be simpler for coding

For example, I'm considering modeling some fluid flow through the volume that interacts with the particles

My thought was that if I can get a handle on defining the problem I could just swap out the linear system for different layers of interaction

Then, since it's a game, set loose tolerances for convergence and not have too fine of a mesh, etc.

[fathom rain]

no idea how to do that for that setting (I have only developed DG codes for airplane design...no interaction)

[raw sand]

What's DG code?

[fathom rain]

discontinuous galerkin, higher order finite volume

so it is fem but with discontinuous elements

https://en.wikipedia.org/wiki/Discontinuous_Galerkin_method

[raw sand]

My eyes are clouding over from the greek letters omelette 😵‍💫 I'm seriously doubting I'm gonna be able to figure this out without taking some courses

[hoary light]

Can you provide a more general description of what you are trying to achieve?

And keep in mind the XY problem as you pose your questions.

[wooden kiln]

I'm a little confused about the following problem:

I numerically integrated the function from a), but got a convergence rate of roughly order 2.5, which doesn't make sense

[vapid plume]

@wooden kiln Fourth derivative of the function is causing problems

Not continuous over [0,1]

[wooden kiln]

But even the second derivative is unbounded, so why isn't it 1st order?

[vapid plume]

First derivative is bounded

second derivative isnt though

@wooden kiln I'm looking at an example, integral of sqrt(x) over [0,1]. Even the first derivative isn't bounded but simpson's rule is O(h^1.5).

[vapid plume]

Book: Parallel multigrid waveform relaxation, Stefan Vandewalle. Chapter 3, section 3.2, numerical experiments, in case anyone wants to see what I'm doing.

[raw sand]

My game is 2D, I was considering adding in some interactions with particles that resemble real world physics, electrostatics, magnetics, fluid flow(where the particles have a force applied to them from a fluid).

I thought it'd be cool if, for example with the fluid flow, there is turbulence, etc. Writing code to "pretend" I'm doing fluid dynamics would be more complicated. Or, for example, with electrostatics, coding the effect of a wall as an insulator or a conductor wouldn't be trivial.

So, I thought if I can incorporate FEM modeling, even if inaccurate, I could get a lot by using pre-existing FEM libraries.

[hoary light]

What is interacting with the particles? Usually modeling fluids with particles is the opposite of modeling fluids with finite elements

[wooden kiln]

I’m just going to say this right now, as someone that spends 8 hrs a day on fluid flow simulations, that if you try adding particles to fluids, your game will not run

Unless you completely cheat and fake it

[silk drum]

Has Dwarf Fortress been looked into? I think they started with something along the lines of finite volume methods to get water and lava flowing between cells.

[lavish stratus]

I'm doing an image compression project using matlab I've done the code how I'm pretty sure he wants it but 10% relative error is looking horrible like here's my plot but here's the image at k=11 and the original

[native atlas]

@neon herald

<@&268886789983436800>

[raven dew]

Don't ping moderators for maths help

And this is not numerical analysis

[native atlas]

Okay, sorry.

Have a good one

[brave saddle]

is the order always contained in the set of all natural numbers?

[knotty cloud]

the order α does not need to be contained in the set of natural numbers, as i see, well, i have seen this problem tho, it must here be a non-negative real number

hope that helps

[primal matrix]

how can i do this process of checking stability with a system of ODEs?

i can translate it if needed

[quiet crystal]

good reccs for books delving into numerical methods/numerical analysis for PDES?

[fathom rain]

depends on what pdes and what type of method you want. i dont think there are any good general books.

[fathom rain]

look at the eigenvalues https://eng.libretexts.org/Bookshelves/Industrial_and_Systems_Engineering/Chemical_Process_Dynamics_and_Controls_(Woolf)/10%3A_Dynamical_Systems_Analysis/10.04%3A_Using_eigenvalues_and_eigenvectors_to_find_stability_and_solve_ODEs

[brave crypt]

When do you ever prefer polynomial interpolation over spline interpolation?

And when Lagrange over newton?

Why do we learn all if apparently some are worse than others

[primal matrix]

ahh so its essentially the same thing as when you’re not solving numerically?

[primal matrix]

usually the first methods you see even tho you might not use them later for not being as good as some other methods, they are usually good understand what those methods are all about, it’s always better to start off with something simple

their construction is usually something very simple compared to some other methods, thats why its good to learn them first

[brave crypt]

so spline is just better than newton and lagrange?

or does spline have any negative sides too

[plucky kayak]

i think you derive finite difference approximations to derivatives from newton's interpolation?

which you can then use to solve (partial) differential equations

I think splines have the issue that you need to choose somewhat arbitrarily boundary conditions on the polynomials or their derivatives at the edges of your entire interpolation interval because the linear system ends up underdetermined

not sure if it's addressed in the literature

oh and also I'm pretty sure you can compute newton's interpolation polynomial in O(N^2) operations if you're clever, but blindly solving the linear system for splines has O(N^3) complexity

[brave crypt]

I see

thanks guys

[primal matrix]

from what i know, i could be wrong, splines were created to try to avoid runge’s phenomenon, but its a more complex algorithm, so its better in a way but not always

[next garden]

Does anyone knoe of a simple to follow proof, either in writing or video, that explains how and why the fourier basis spans L2 ?

[grave spoke]

There is also something called barycentric Lagrange interpolation, which is more efficient than Newton interpolation if the distribution of nodes is well-known, e.g. equidistant nodes or Chebyshev nodes.

[cobalt sky]

https://math.stackexchange.com/questions/316235/proving-that-the-fourier-basis-is-complete-for-cr-2-pi-c-with-l2-norm

[wild lily]

Hello, is this the right place to ask about BFGS?

[errant wagon]

yep

[wild lily]

Could I ask about the textbook then? Jorge and Nocedal numerical optimization. Or is it too simple for this channel?

[wild lily]

I'm having trouble figuring out why the first tangent is (0,-1)

Why the image doesn't seem to match the sequence of vectors

Or rather, I can't see the sequence

And why the second tangent is (0,\alpha)

I'm assuming the second tangent is actually (0,1) and he's replacing the sequence t_k with t_k/\alpha to get that.

It is not obvious to me at this point.

[errant wagon]

Where is problem 12.9

[wild lily]

@team132 here it is.

[eternal drift]

You are approaching the point (-sqrt(2), 0) on the circle from below

[eternal drift]

These are all tangents in that direction (it's a tangent cone so positive rescalings are also tangents)

So the tangents from below are actually (0, -alpha) for alpha in R+

The tangent cone is union of these two sets leading to the (0, d) with d in R

[wild lily]

That doesn't make sense to me since if I'm approaching from below, wouldn't my tangent be (0,1)^T?

[eternal drift]

Look at the definition of the tangent, it's z - x divided by some positive value

So suppose z is just below x, for example (-sqrt(2), -eps), then the tangent value is (-sqrt(2), -eps) - (-sqrt(2), 0) = (0, -eps)

[wary oriole]

hey guys

i have a question regarding finite difference for laplace

for four interior points it's easy to drive the system

but for 50*50 points, what's the system looks like? I have trouble to get A and b

please help. i got all iterative algos ready to go but i cannot drive the system. The due date is coming :opencry: Please Please

[fathom rain]

kron(toeplitz(n1,[2,-1]),I(n2))+kron(I(n1),toeplitz(n2,[2,-1]),)

the look of it depends on what n1 and n2 are

(and you can use kron for b too)

[wary oriole]

can you explain a bit how you drive it?

what's n1, n2?

[wary oriole]

hey i tried your code. N=50 for both dimension. It's not correct

array([[100., -1., 0., ..., -0., 0., -0.],
[ 0., 50., 50., ..., -0., 0., -0.],
[ 0., 0., 50., ..., -0., 0., -0.],
...,
[ 0., 0., 0., ..., -1., 0., -0.],
[ 0., 0., 0., ..., -1., -1., -0.],
[ 0., 0., 0., ..., -0., 50., -2.]])

A = np.kron(toeplitz(N, [2,-1]), np.identity(N)) + np.kron(np.identity(N),toeplitz(N, [2,-1]))
in python

this has dimension 50*100

[woeful jetty]

so someone told me to ask this here. is Lax Wendroff supposed to look this weird? the problem is to solve burger's eq with u(x,0)=exp(-16x^2) and periodic boundary conditions

[fathom rain]

n1 and n2 are the sizes in the two dimensions dimensions

I don't use python, the above was just pseudo code. what i wrote is correct except i left a comma in the end, soo here is correct version in julia with n1=n2=3: kron(toeplitz(n1,[2,-1],[2,-1]),I(n2))+kron(I(n1),toeplitz(n2,[2,-1],[2,-1])) 9×9 Matrix{Int64}: 4 -1 0 -1 0 0 0 0 0 -1 4 -1 0 -1 0 0 0 0 0 -1 4 0 0 -1 0 0 0 -1 0 0 4 -1 0 -1 0 0 0 -1 0 -1 4 -1 0 -1 0 0 0 -1 0 -1 4 0 0 -1 0 0 0 -1 0 0 4 -1 0 0 0 0 0 -1 0 -1 4 -1 0 0 0 0 0 -1 0 -1 4

(I have my own toeplitz generator)

and this would ne for n1=3 n2=4

12×12 Matrix{Int64}:
  4  -1   0   0  -1   0   0   0   0   0   0   0
 -1   4  -1   0   0  -1   0   0   0   0   0   0
  0  -1   4  -1   0   0  -1   0   0   0   0   0
  0   0  -1   4   0   0   0  -1   0   0   0   0
 -1   0   0   0   4  -1   0   0  -1   0   0   0
  0  -1   0   0  -1   4  -1   0   0  -1   0   0
  0   0  -1   0   0  -1   4  -1   0   0  -1   0
  0   0   0  -1   0   0  -1   4   0   0   0  -1
  0   0   0   0  -1   0   0   0   4  -1   0   0
  0   0   0   0   0  -1   0   0  -1   4  -1   0
  0   0   0   0   0   0  -1   0   0  -1   4  -1
  0   0   0   0   0   0   0  -1   0   0  -1   4```

i guess you might see the pattern

(block sizes are determined by n2 in this setting and n1 determines number of blocks)

[fathom rain]

https://en.wikipedia.org/wiki/Gibbs_phenomenon

[fathom rain]

this is using scipy toeplitz with N=3 >>> np.kron(sp.linalg.toeplitz([2,-1,0],[2,-1,0]), np.identity(N)) + np.kron(np.identity(N),sp.linalg.toeplitz([2,-1,0],[2,-1,0])) array([[ 4., -1., 0., -1., -0., 0., 0., 0., 0.], [-1., 4., -1., -0., -1., -0., 0., 0., 0.], [ 0., -1., 4., 0., -0., -1., 0., 0., 0.], [-1., -0., 0., 4., -1., 0., -1., -0., 0.], [-0., -1., -0., -1., 4., -1., -0., -1., -0.], [ 0., -0., -1., 0., -1., 4., 0., -0., -1.], [ 0., 0., 0., -1., -0., 0., 4., -1., 0.], [ 0., 0., 0., -0., -1., -0., -1., 4., -1.], [ 0., 0., 0., 0., -0., -1., 0., -1., 4.]])

[wary oriole]

hi if N=2 i get
array([[ 4., -1., -1., 0., 0., 0.],
[-1., 4., -1., -1., 0., 0.],
[-1., -1., 4., 0., -1., 0.],
[ 0., -1., 0., 4., -1., -1.],
[ 0., 0., -1., -1., 4., -1.],
[ 0., 0., 0., -1., -1., 4.]])

this doesn't look like matrix in the book tho. should it be producing
array([[ 4., -1., -1., 0.],
[-1., 4., 0., -1.],
[-1., 0., 4., -1.],
[ 0., -1., -1., 4.]
])

like the book

also N=50 you way will product 150150. not 5050

[wary oriole]

the dimension doesn't match. so N=50 in my problem means each axes has 50 point including boundary

so i guess 48*48 equations?

"Next the spatial and temporal domain are discretized. The x- and y-axes are divided
in intervals with length δx. Assuming that we have N interval in the x- and y-
directions, we have δx = 1/N , and x = iδx, y = jδx where i, j ∈ (0, 1, 2, . . .N ).
Furthermore we assume a small time interval δt, so that t = kδt where k ∈ N"

[fathom rain]

no: >>> np.kron(sp.linalg.toeplitz([2,-1],[2,-1]), np.identity(N)) + np.kron(np.identity(N),sp.linalg.toeplitz([2,-1],[2,-1])) array([[ 4., -1., -1., -0.], [-1., 4., -0., -1.], [-1., -0., 4., -1.], [-0., -1., -1., 4.]])

that is for N=2

you have to manually modify the arguments of toeplitz function to match your N add as many 0s as you need

there is probably some convenience function in numpy to do that

if you want a 48x48 system you can do n1=8 and n2=6 for example

[fathom rain]

You can for example do

n1=2
n2=2
v1=np.zeros(n1)
v2=np.zeros(n2)
v1[0:2]=[2,-1]
v2[0:2]=[2,-1]
np.kron(sp.linalg.toeplitz(v1,v1), np.identity(n2)) + np.kron(np.identity(n1),sp.linalg.toeplitz(v2,v2))```

[wary oriole]

oh finially i think i got it. Thank you for the help. I still would like to drive this matrix by hand but im running out of time so i'll just take it for granted

👍

[fathom rain]

this is called the kronecker sum https://mathworld.wolfram.com/KroneckerSum.html

[wary oriole]

i meant drive enough of these 48*48 equations in order to see the pattern. like the book did with 4 equations

[fathom rain]

should say n1=3 and n2=4 on top right 😛

[fathom rain]

I would not recommend doing it like that book, very inconvenient. simpler to look at it lika a combination of 2 1D laplace matrices

[sudden adder]

anyone using wavelets methods to solve ODEs, PDEs, Integral equations, Integro-Differential Equations numerically. I would really like to talk. Please response.

[wary oriole]

hi @fathom rain thank you for your answer

[ruby pike]

Numerical LA question: Is there an algorithm that takes an upper triangular matrix over the complex numbers as input, and spits out a bidiagonal matrix which it's similar to?
Don't misunderstand my question: The usual "bidiagonalisation" algorithms (which occur in the context of computing SVD for instance) don't output a similar matrix. They only preserve the singular values, which is not what I want.

[fathom rain]

assume you have distinct eigenvalues, then you can just take any upper diagonal you want since it will be similar

and if you want the specific tranformation you can compute it using the eigenvector matrices for the matrix you desire

is there a reason to have it bidiagonal?

[rich flax]

Given SVD of A
Find SVD of matrix
( I A)
(A^T 0)
where 0 and I are all zero and identity matrix

Can anyone give me hint on how to find singular values of above matrix?

[ruby pike]

Substitute A = USV^T. Then multiply by diag(U^-1, I) on the left. Then multiply by diag(U, I) on the right. Then multiply by diag(V,I) on the left. Then multiply by diag(V^-1,I) on the right. The result gives [I, S; S, 0]. The singular values have been preserved, as each manipulation was an orthogonal equivalence.

Permute the rows and columns so that you have a direct sum of 2x2 matrices of the form:

⎡1   sᵢ⎤
⎢      ⎥
⎣sᵢ  0 ⎦

The singular value set that is sought is precisely the union (with possible repetition) of the singular values of the above 2x2 matrices.

The singular values of those 2x2 matrices happen to be:

⎡      __________________________        __________________________⎤
⎢     ╱          ___________            ╱          ___________     ⎥
⎢    ╱          ╱     2                ╱          ╱     2          ⎥
⎢   ╱     2   ╲╱  4⋅sᵢ  + 1    1      ╱     2   ╲╱  4⋅sᵢ  + 1    1 ⎥
⎢  ╱    sᵢ  + ────────────── + ─ ,   ╱    sᵢ  - ────────────── + ─ ⎥
⎣╲╱                 2          2   ╲╱                 2          2 ⎦

I think that's the answer to your question/

That's more complicated than you'd think given the form of the matrix.

[fathom rain]

more readable form (4*s^2 + 1)^(1/2)/2 \pm 1/2

6-element Vector{Float64}:
 2.3355261064228325
 1.3355261064228316
 1.1429783479872964
 1.0825415690164035
 0.14297834798729703
 0.08254156901640394

julia> sqrt.(4*svd(A).S.^2 .+1)/2 .-1/2
3-element Vector{Float64}:
 1.335526106422832
 0.14297834798729703
 0.08254156901640397

julia> sqrt.(4*svd(A).S.^2 .+1)/2 .+1/2
3-element Vector{Float64}:
 2.335526106422832
 1.142978347987297
 1.082541569016404```

[ruby pike]

An even slicker solution. The reduction to the matrix:

⎡I  S⎤
⎢    ⎥
⎣S  0⎦

uses only orthogonal similarities, so preserves both eigenvalues and singular values simultaneously.
It's also worth noting that that matrix is equal to $I \otimes diag(1,0) + S \otimes \begin{pmatrix} 0 & 1 \ 1 & 0\end{pmatrix}.$
The tensor product operation is commutative up to matrix similarity. Thus my reduction to a direct sum of 2x2 matrices follows. Take their eigenvalues, and (as @fathom rain pointed out) the absolute values of these are the singular values.

The "tensor product" here is really the Kronecker product: https://en.wikipedia.org/wiki/Kronecker_product

[ruby pike]

Bump

[hoary light]

A bidiagonalization (via matrix similarity) is equivalent to computing jordan canonical form. In general this is numerically unstable, but in principle you could still try to do this, especially since you already know the eigenvalues through the triangular form that you are given. Are you familiar with computing jordan canonical forms?

[ruby pike]

A bidiagonalization (via matrix similarity) is equivalent to computing jordan canonical form
No it isn't. Why?

I don't care about "almost" getting the Jordan form. Give me the Jordan form if you think it's possible.

[ruby pike]

I'm aware of links to the Jordan Normal Form. That's why I think it's an interesting question.

I don't see how to find the JNF of a bidiagonal matrix in a computable/well-conditioned manner.

[ruby pike]

I think that the set-valued function
Upper triangular -> Bidiagonal matrix
is lower hemicontinuous. This is not true for the Jordan Form. If that's true, they are not at all equivalent problems.

Not upper hemicontinuous though.

[lilac crane]

for 2x2 it seems true

since in that case a bidiagonal is similar either to the jordan block (a 1; 0 a) or to a diagonal

anyway

[pine jetty]

memorylessfunctor

[lilac crane]

the jordan condition number is max(delta, 1/delta)

well ok like this is already bidiagonal

but i think you can modify this

[pine jetty]

memorylessfunctor

[lilac crane]

@ruby pike

[lilac crane]

the bidiagonalization condition number of this blows up

if im not mistaken

intuitively it should

[ruby pike]

No it doesn't. Like any non-scalar matrix, it has infinitely many bidiagonal forms

[ruby pike]

I don't know what that means

I think this is wrong

What's a "bidiagonalisation condition number"?

[hoary light]

I guess you are right. I was thinking that given an (upper) bidiagonal matrix one can compute the vectors for the Jordan cycles quickly by backward substitution. However, due to needing to pivot (since the only interesting case is when there are multiple eigenvalues, so multiple diagonal entries of A - lambda * I will be zero) it is not actually any simpler.

[still forum]

Suppose I'm solving a Laplace/poisson equation on a square domain. When considering Neumann boundary conditions (like on the left edge of a square, $\frac{\partial u}{\partial x} = 0$), should I use $O(h^2)$ order accurate approximation on the boundary if everywhere on the interior I have $O(h^2)$ order accurate approximation of $u_{xx} + u_{yy}$?

[pine jetty]

jamiecjx

[still forum]

or is it ok to use first order forward difference

how does it affect the overall accuracy of the solution?

[pine jetty]

nicholas zero nine eight

nicholas zero nine eight

[still forum]

makes sense

thanks a lot

[lilac crane]

|B||B^-1| for B which conjugates the matrix to a particular bidiagonal

i guess you are worried this wouldnt be well defined?

[next garden]

Why does this imply non continuity?

[plucky kayak]

the vector approaches zero, but it's image under L approaches some non-zero vector with norm 1

[wary oriole]

guys

for iterative methods of anything really. For example jacobi method for solving linear system

how do you log your iterations? in python i print them like this
print('iteration: ' + str(k) + '\t' + 'error norm: ' + str(error))

is there a better way? i mean not just log in a file, but use logging in python?

in a production environment

[fathom rain]

just save in csv or some other format without the "iteration" etc, just keep a readme file what different columns represent, (or use dataframes or something similar)

[fathom rain]

you can also use something like this, I assume similar systems exist for python https://github.com/JuliaDynamics/DrWatson.jl

[wary oriole]

👍

[wary oriole]

btw, since last regular boundary condition for my finite difference laplace equation. Our teacher let us put sinks (0 temperature) boundary inside the square. I kinda drived the matrix

three sinks solution, the matrix is annoying to drive. i put specific zeroes in regular matrix you gave me

[next garden]

Why does this imply discontinuity?

I am not following

[lilac crane]

L(0) = 0

but L(vn) does not converge to 0

[next garden]

Ahhhh that makes sense

Ty

[wooden kiln]

In my graduate CFD course, we are talking about Von Neumman Stability Analysis. From my professor's notes, he has the following:

My question is, where did the i (index, not imaginary unit) terms go on the left hand side, and where did the k (wavenumbers) terms go on the right hand side?

[fathom rain]

$g^ne^{ik\theta}$ is not zero so you have $((g-1)/dt-ia\sin(\theta)/x)=0$ then just rearrange the terms

[pine jetty]

sverek

[fathom rain]

your k=i

quite sloppy notation overall 😛

so just replace i with k in the Stability of:... equation (and add subscript k on the time derivative part )

[inland owl]

Hi I study English language at uni but I’m doing truth tables and really struggling can someone PLEASE help me with them ?

[gloomy gorge]

#1021175428326633542

this is the wrong channel

[still forum]

I've been using SOR on a large linear system (from an elliptic pde) and I'm interested in finding the optimal parameter. However afaik the matrix, while sparse, doesn't have nice structure (not symmetric, tridiagonal, block tridiagonal) as the domain I took it from is a weird polygon shape and has a mixture of boundary conditions. How does one then minimise the spectral radius of the iteration matrix?

I believe the matrix is positive definite (it's pretty diagonally dominant) so my thought was to do a bunch of power iterations or something

but theres presumably too many clustered eigenvalues for that to work well

[hoary light]

If you have an idea of the (numerical) multiplicity of the max absolute value eigenvalues, you can try using a block power method

There are many randomized methods, often variations on "multiply the matrix with a random set of vectors to learn properties of the matrix" that can work

Do you need to use SOR? It is not normally used for very large systems

[still forum]

it was about$ 10^4$ unknowns from like a 100x100 grid roughly, and it was a coursework exercise so yeah I had to lol

I did some asymptotics and assumed the spectral radius of the matrix was $\rho \approx 1 - \frac{\pi^2}{N^2}$ and got something that sped up convergence by a rough factor of $N$ which was good enough for me

[pine jetty]

jamiecjx

[still forum]

back to this problem

I want to see how big I need to make my grid until the solution I obtain is grid independent (up to some tolerance), and I don't have access to the analytic solution.

My thought was to take a sequence of increasing grids which have some shared coarse grid $C$ and interpolate the values onto $C$, whilst measuring the difference between successive grid sizes (in L2 norm). But I'm having trouble finding a sort of standard test for this (Leveque doesn't seem to have anything on this)

[pine jetty]

jamiecjx

[pure lake]

hello

i need some views on this series

[1, 3, 7, 11, 15, 23, 27, 31, 39, 47, 59, 63, 71, 79, 91, 95, 103, 111, 123, 127, 155, 159, 167, 175, 191, 199, 207, 219, 223, 231, 239, 251, 255, 283, 287, 303, 319, 327, 347, 359, 367, 383, 411, 415, 423, 447, 463, 479, 487, 495, 507, 511, 539, 543, 559, 575, 583, 603, 615, 623, 639, 667, 671, 679, 703, 719, 735, 743, 751, 763, 767, 795, 799, 815, 831, 839, 859, 871, 879, 895, 923, 927, 935, 959, 975, 991, 999, 1007, 1019, 1023, 1051, 1055, 1071, 1087, 1095, 1115, 1127, 1135, 1151, 1179, 1183, 1191, 1215, 1231, 1247, 1255, 1263, 1275, 1279, 1307, 1327, 1343, 1351, 1383, 1407, 1435, 1439, 1471, 1487, 1503, 1511, 1519, 1535, 1563, 1567, 1583, 1627, 1639, 1647, 1663, 1691, 1695, 1703, 1727, 1743, 1767, 1775, 1787, 1791, 1819, 1823, 1855, 1863, 1883, 1895, 1903, 1919, 1951, 1959, 1983, 2015, 2031, 2043, 2047]

i dont want to induce any confirmation bias

i need blind views

[fathom rain]

depends on grid numbering but assuming you have $n_k=(n_0-1)2^k+1$ internal nodes then you hav the indices $j_k=(k+1)*(1:n_0)$ in each level, this is 1D but you can use kron again for this. then compare 2 solution levels for the difference is an approximation of the error. once size is enough you have your solution

[pine jetty]

sverek

[pearl canopy]

not the channel for this

[slate sinew]

does anyone have notes or resources that can help me understand the finite precision system? i'm kind of getting lost as to when i need truncate/round when solving lineal systems with the Gauss elimination method.

Example of an exercise:

Consider the finite precision system F (10, 4, −∞, ∞).

Solve the linear system by the Gaussian elimination method using rounding as the method
approximation in the finite precision system.

[minor tree]

Hi Hi
does anybody of yall know about good sources for the Fourier spectral method.
Most of the ones i find only consider certain PDEs as examples. For context i am trying to implement that method in Python for general PDEs with periodic BC and am unsure as to what is a good approach for implementation.

[fathom rain]

https://epubs.siam.org/doi/book/10.1137/1.9780898719598 maybe

[minor tree]

thank you that looks promising :)

[heady mesa]

Boyd's book

Spectral Methods by Boyd is the Bible of Spectral Methods

[wild lily]

Does anyone know how I should interpret Jorge and Nocedal's algorithm 16.5? I've been trying to implement this but I can't seem to understand it. I'm completely unsure as to what I'm getting wrong

[wild lily]

I'm trying to use algo 16.5 for section 18.6, for algo 17.4.

i wanted to do with just the cauchy point first. it's much slower, but i want to get the inner loop right first before i try something else

but it seems like the penalty instantly forces my iterate to 0

[dull wraith]

Hey i am new to numerical analysis.can anyone suggest from which book should i start

[silk drum]

Numerical analysis is a pretty wide field and you'll want a text which covers the particular subject matter of interest to avoid wasting motivation on interesting but less useful areas. In particular, you may have applications in mind and want a text targeted towards getting you up and running as soon as possible. Other people will want a theoretical and foundational text if applications are not immediately relevant.

[silk drum]

https://dl.acm.org/doi/10.1145/103162.103163

Although the main aim is to address floating precision, from memory, it helped me to grasp fixed precision by way of comparison and analogy. Off the top of my head, it would be truncate/round after each discrete/atomic instruction or operation, so usually a floating point operation such as +, -, ×, ÷. Actual computers have FMA and GPU DMMA so multiple operands and less rounding. An important detail is that libraries can carry out algorithms to minimise loss of precision and accuracy, a simple example is to accumulate from lowest magnitude first to avoid overflow or high magnitude values impacting the contributions from low magnitude values.

[untold crag]

Hi guys!

Let’s assume that I have a given time series with equally spaced observations. I want to fit an Ornstein Uhlenbeck process on it (i.e. estimate the three parameters). How to proceed? I know Maximum Likelihood Estimation is a good option but if any of you has some ressources for me to learn about how this procedure works, I’d be more than happy. My second question is, let’s assume now that we have a more “ugly” time series, non equally spaced observations, how does it change the procedure?
In the end, I would implement everything using Python which seems ok for this.
Let me know guys!
Thank you

[pulsar cipher]

Hi, I'm reading into numerics and I dont understand this part:

[pine jetty]

Sciencenjoyer

[plucky kayak]

which part

[pulsar cipher]

First I'm not sure which $x_ {n+1}$ is used in the lower condition and why did they define $\tilde{f}(x,x_{n+1}) := x_{n+1}$ ?

[pine jetty]

Sciencenjoyer

[pulsar cipher]

Basically I don't see how those two problems are equivalent

[plucky kayak]

"lower condition"?

[pulsar cipher]

what comes after "subject to"

[plucky kayak]

there's only one x_{n+1}

you had a problem in n variables and you added another one

[pulsar cipher]

but here x_{n+1} is used to define itself

[plucky kayak]

the defined a function of variables x and x_{n+1} which is equal to x_{n+1}

this isn't the definition of x_{n+1}

x_{n+1} is just a new arbitrary scalar variable

[pulsar cipher]

oh ok, I get that now, thanks

But I don't see why both problems are equivalent

[plucky kayak]

well suppose you want to minimize $f(x)$ with respect to $x$ and $x^\star \in \mathbb{R}^n$ is the minimizer and also consider another problem $t \to \underset{x\in \mathbb{R}^n, t \in \mathbb{R}}{\min}$ subject to $f(x) \leq t$ and $t^\star$ is the optimal value of that problem, then I claim that the pair $(x,t) = (x^\star,f(x^\star))$ is a solution of the second problem

[pine jetty]

Transparent Elemental

[plucky kayak]

if that wasn't the case, then I could decrease value of t, but that's not possible because we assumed that f(x^star) is the minimum possible value

on the other hand if you consider second problem and it's solution (x^star, t^star), then x^star is a solution of the problem of minimizing f(x) because we found such (x,t) that ensure that f(x) is as small as possible from the constraint that we had

[pulsar cipher]

I see that t* is an upper bound for f(x)

Sorry but what is "the minimizer"?

[plucky kayak]

the point at which the function achieves it's lowest value

[pulsar cipher]

ohhh. This should have been super obvious

At least I can see that now. Thank you very much

[brave crypt]

I need some help regarding (PINN-Physics Informed Neural Network) to make an AI that solves a non linear differential equation. Someone with enough knowledge about this please DM.

[atomic meadow]

im in grade 9

how can i help?

[brave crypt]

do u know what PINN is

have u studied enough ML and its application in Physics

@atomic meadow

[atomic meadow]

im not good at it

take help of chatgpt

[brave crypt]

thanks for nothing,. Chat gpt is lame, good for kids only

[atomic meadow]

sorry if i couldnt help

[brave crypt]

it's okay. have a nice day! Sorry it that was rude, I am just stressed

[atomic meadow]

you are going to release the AI or you are making for someone else?

[brave crypt]

it's a college porject

not commercial or anything

also, I will make it for myself ofc

maybe use it as a general tool

[atomic meadow]

which programming language your using?

ig, Python

@brave crypt

[brave crypt]

yes

python

[pale gazelle]

The probability density function of the chi-square distribution involves the gamma function (Γ). From what I've researched so far, it seems that in order to understand what the gamma function is and how it works, I need to have a thorough understanding of how complex numbers work (which unfortunately, I don't). Unlike functions defined for real numbers where you simply input a value x and get a corresponding output that depends on the function definition, it's rather unclear what actually happens when you input an x to the gamma function. I have a good understanding of derivatives, partial derivatives, integration, and some intermediate-level calculus topics. If someone were to ask me, "What does the function f(x) = 2x do?", I would answer, "It takes in a real number x and returns a value that's twice the value of x". Similarly, could anyone explain what the gamma function does? For context, I only want to understand the gamma function to get an intuitive understanding of the PDF of the chi-square distribution.

[plucky kayak]

in order to understand what the gamma function is and how it works, I need to have a thorough understanding of how complex numbers work
no you don't

[pale gazelle]

Then, could you please explain how that function works?

[plucky kayak]

this is what the gamma function does to a number z, whether it's real or complex

[pale gazelle]

What's the purpose of it?

[plucky kayak]

well it's the same purpose as with any other function

[pale gazelle]

Also, what are t and z?

[plucky kayak]

this integral comes up a lot so they gave it a name

[plucky kayak]

the integration variable and function input

[pale gazelle]

Since I'm asking this in #numerical-analysis, I suppose the person answering here does not know the nooks and crannies of the PDFs used in hypothesis testing in statistics. But, in case you happen to have a nuanced understanding of the PDF of the chi-square distribution, could you please explain the purpose of the gamma function in it?

[plucky kayak]

this is just what the PDF works out to be, if you want a proof for why the PDF is this and not something else then you can search for the proof

there isn't always a specific purpose like you might think

it's not like somebody thought about it for a long time and suddenly he "understood" that there must be a gamma function involved, no

chi squared is not unique in using gamma function either

[pale gazelle]

Oh, alrighty. I shall look into its derivation and proof to get a more nuanced understanding of why it's defined the way it is defined. Thanks!

[hollow perch]

Does anyone have an intuitive understanding why many problems on 2d arrays can be solved in polynomial time, but generalizing to nd arrays is np conplete? E.g. the linear sum assignment problem, or decomposition into a sum of rank 1 matrices/tensors

[hoary light]

To my understanding the main issue is that several concepts concerning matrices do not have a nice generalization for multilinear maps. In particular, in your aforementioned example of decomposition into lower-order tensors, since a matrix can be viewed as both a representation of a bilinear form and a linear map on a vector space, all the notions of rank coincide, namely with the SVD. In higher dimensions, more than one concept of rank exists and these extra degrees of freedom makes computing any such notion of rank very difficult.

Perhaps one fuzzy intuition to use is that in the case of a matrix (bilinear form), computing something like rank amounts to some kind of linear or quadratic problem, which is convex. However, at higher orders, problems such as approximating one tensor by another is no longer convex because it's like a higher order polynomial.

[hollow perch]

The convexity argument sounds like it should be correct because matrices are n=2 dimensional, but I don't completely buy it as the rank function is not convex

[hoary light]

The rank function is not convex, but instead think of the problem of approximating a matrix by a sum of rank 1 operators

In particular, think of how maximizing and minimizing the Rayleigh quotient gives you certain eigenvalues (which are of course nearly like the singular values in the hermitian case)

This doesn't really explain why the multilinear case is hard, but it does explain why the matrix case is easy

[hollow perch]

I'm still not convinced. Sure the maximum and minimum singular values are a convex function of the matrix A, but the singular value map that maps A to the vector of singular values sigma is not convex, and the SVD map that maps A to (U, sigma, V) is not even continuous, where U and V are the matrices of singular vectors

[hoary light]

The SVD map is not really the correct way to think of it, but rather the problem of minimizing |A - A_k | where A is the matrix in question and A_k is a rank-k matrix.

When treated as a rank-revealing problem then each step is quadratic via the Rayleigh quotient.

In contrast, for higher order tensors, we don't quite have the same formulation. Tensor rank doesn't quite line with the generalization of SVD that we might want. At least that is my impression.

[plucky kayak]

Hillar, Lim, "Most Tensor Problems are NP-Hard" from 2013 might be of interest

the more I learn about tensors the more I'm disappointed

[hoary light]

from seminars I've been to, they still play a useful role in applications. It's just that generically you can't solve a lot of tensor-related problems easily. But for tensors that arise in applications there may still be additional structure, empirically

[knotty shadow]

Ja, check, I'm $\textit{not}$ trying to make shit here, but you were not the greatest person in that. You can give me trouble, I don't care, but I don't think that you were a nice person there. Anyway, what you're asking belongs almost entirely (I always take liberties with my ability to claim almost-$something$), within computer science. I'm aware of the connections between neural networks, genetic algorithms, etc. and numerical analysis, but I don't think that you're question belongs in numerical analysis.

Look, I'm not trying to be an arsehole to you. I just had some words to type and now that they're typed we move on with life, in my opinion.

I cannot guarantee that anything will come of this, but I'll ask some colleagues on Monday about PINN. Otherwise, there are good libraries for Python which can deal with DEs and, at the worst, you can use numerical techniques (there are good libraries for that, too).

[pine jetty]

Saurus

[knotty shadow]

The Gamma Function is the function which results as the analytic continuation of the factorial. For natural numbers (I don't include 0, so $\mathbb{N} = \left{ 1, 2, 3, ... \right}$), the Gamma Function (here, denoted $\Gamma \left( \cdot \right)$), satisfies $\Gamma \left( n \right) = \left( n - 1 \right) !$. This clearly gives that the Gamma Function applied to a natural number is, approximately, the factorial of the number (there's that $-1$ which appears, but it's close enough). This then allows for a meaningful definition (and understanding), of the factorial of a real number and, by the analytic continuation, a "factorial" of (almost) any complex number.

Are you good to go with that explanation?

[pine jetty]

Saurus

[knotty shadow]

The story with tensors is, currently (by currently, I mean for the last two years), the bane of my existence. There are many definitions of the tensor product. Some are constructive (as in, it gives a way to do a tensor product), but many are simply existence results. Now, the tensor product of two vector spaces over the same field exists and it is unique (this is from the universal property, but there is nothing which shows how to construct the tensor product). The tensor product which features in, let me claim it to be, useful places, is, in fact, the Kronecker product.

Importantly, matrices describe linear operators of one finite dimensional vector space into another finite dimensional vector space (same fields). In fact, each linear operator (between finite dimensional vector spaces over the same field), has a matrix representation and the converse is true. This, very importantly, shows no relation between matrices and multi-linear maps (multi-linear of order strictly more than one), unless you do those funny things that the physicists do (for example, the vec function - then matrices have meaning as multi-linear maps).

The best seen use of tensors is in general relativity (actually, the way things are done there is the use of the Kronecker product, but the Kronecker product is a fantastic representation of the tensor product for finite dimensional spaces), but this, the way that it is treated (and, hence, the fantastic use of the Kronecker product), boils down to a bunch of index notation and multiplication of members of specific multi-dimensional arrays.

I'm pretty drunk and I forgot what the purpose was of why I started typing a response, so I'll just stop here before I move to the abstract (algebraic) tensor product and then the Fremlin tensor product (that guy is what I used in a paper a year or two ago). There's also the story about the tensor product being used in pseudo-Riemannian geometry not actually being a tensor product (welome the Kronecker product), because the product used isn't a true tensor product (it isn't complete, etc.).

[stoic sun]

This question has to do with pade approximations, and using them in software. I'm am unsure if its better to ask my question in here, #geometry-and-trigonometry or elsewhere, so if there is a better channel for this question, please let me know.

Before asking, let me just say: Im really only going to be using the Pade approximations for approximating the function arctan. I might talk about the process for making a pade approximation from a general taylor series, but in practice I only really need to apply it to the taylor series for arctan. I mention this because maybe this simplifies things, or is possible where the general case is not.

In essence, my question is this:

Can the conversion of a taylor series to a pade approximation be made iterative. By iterative, I mean taking it one term at a time.
Looking at a few videos that shows the process, it would appear that: No, you'll have to know how many terms you are going to use up front before beginning to work out the pade equation.

Some background: I have written and maintain an arbitrary precision, base-10 floating point number type library in C#. A while back, I added trigonometric functions to the library, which sets it apart from other arbitrary floating point libraries in the same language. I implemented this by using the taylor series. It takes as a parameter, the desired level of precision specified as the n in the expression 10^-n, so basically n is the number of zeros past the decimal point you want the result to be good to. It takes the Taylor series term by term, calculates the numeric result, subtracts that result from the numeric result from the prior iteration, and if the difference is less than 10^-n, it returns the result.
This is an iterative process, each iteration further refines the running total, and calculating the next term in the series does not require recalculating the terms up to that point.
I implement many trig functions through this function. Yes, i take care to restrict the inputs, wrapping the inputs at pi or pi/2 or whatever the range of the function in question happens to be.

Well as some of you already know, for certain trig functions at certain inputs, the taylor series is VERY slow to converge. In particular, I am getting very long evaluation times for arctan (and arcsin and arccos, which use arctan under the hood). As you know, the domain of arctan is All real numbers.

What id like to do it replace using the taylor series for arctan with something else with better convergence, and is more consistant with how long it will take to converge for any input.
Looking into this problem, I find suggestions of the following possible solutions:

• chebyshev polynomials
• pade approximations
• continued fractions

chebyshev polynomials use integrals, and i looked into them a little bit, but the algorithm wasn't immediately understood by me, but i didnt put a lot of effort into understanding them. if you feel like chebyshev polynomials are absolutely what I need, then Im sure i can figure it out, i just need to apply more effort. considering what i found when i looked into pade approximations, maybe this warrants another look.

...

Looking into pade approximations, first of all, the only reason I am able to consider doing this algorighm is due to the fact that I have written a symbolic multivariate polynomial library which can be leveraged to do the multplying of the arbitrary polynomials with different coefficients. BUT, it appears i will not be able to work out each A/B term from the associated taylor polynomial term, without first deciding how many terms I am going to need. In particular, one of the steps involve taking the A/B = (taylor series) equations, setting B₀ = 1, multiplying each side by the denominator and arriving at an equation like:
A₀ + A₁x + A₂x² + ... = (1 + B₁X + B₂X² + ...)(taylor series terms)
and from this figuring out what all the A and B coefficients are.
Formulated this way, although I could take the taylor series one by one, id still need to know how many B terms to multiply it by. Perhaps this is not strictly true, but at each iteration i would be essentially adding another A term, B term and taylor term and then perform all the 'missing' calculations that are entailed by the presence of the new terms. I guess this is possible, since, for the right hand side, we are essentially taking a cartesian product of the B terms and taylor terms and adding their products, because that is the case, adding another B term and T(aylor) term doesnt change any of the products we already calculated, we just need to add some more: the new B term times all of the prior T terms, all of the prior B terms times the new T terms and then the new T term times the new B term. So I guess this is possible. I dont know how it will affect the linear algebra problems that follow, but i do not foresee any fundamental problems. Is this right? Is there any barriers to this that im not seeing?

The continued fraction approach. I've written a library that does continued fractions, so this is no challenge. I think this would be my preferred approach if I could wing it. It certainly seems the most straightforward.
...
Actually, I think my main issue here was I was searching for the continued fraction of arcsin or arccos, which I couldnt find.
But thats stupid. Since both arcsin and arccos use arctan under the hood, a continued fraction for arctan would serve me just fine, and I think that was the only example I could find when searching for arcsin and arccos. But duh, I should just use that, the arctan one.

Yeah, this might just work for me instead.

Rubber duck moment.

[cosmic karma]

Hello folks! I’m trying to formulate a discontinuous Galerkin FEM for a linear coupled system with one dependent variable having a non-homogeneous Dirichlet BC. I read that the treatment of this boundary comes as a numerical flux. Is that correct and if so, how do I do it? Do I introduce a new parameter, like the flux would be an artificial viscosity or something to the system? Thank you!
What about the variable with a homogeneous dirichlet BC? Does this one need a numerical flux?

[fathom rain]

just impose the BC weakly at the boundary the same way as you compute fluxes everywhere else?

[cosmic karma]

great thanks! I am using averages for now to make it very simple. Thanks!

[pale gazelle]

Yes, it absolutely makes sense! Thanks for explaining!

@knotty shadow

[dark pumice]

Hi

[upbeat eagle]

help

[pine jetty]

jacobjivanov
Compile Error! Click the reaction for more information.
(You may edit your message to recompile.)

[still forum]

Suppose I'm solving a wave equation (with constant c^2 = 1 for simplicity) but the initial condition u(x,0) = 4x^2 (-0.5<x<0.5), and u(x,0) = 0 everywhere else. (ut(x,0) = 0 as well).

I used 2nd order finite differences everywhere (leapfrog) but oscillations get introduced which is bad. If my initial condition is continuous, the code I wrote works fine. How does one deal with discontinuous wave conditions in finite difference methods?

Or even, if my initial condition is u(x,0) = 1 for (-0.5<x<0.5), i still get random oscillations

despite the wave moving roughly in the right direction

cranking up the resolution of the grid seems to make the oscillations smaller, but is not ideal because I need quite a fine grid

Ok, I thought a bit about the physics

the grid ratio really matters. Since the wave speed is 1, we kinda need the ratio of h and k to be 1

otherwise, the solution doesn't match the grid

so I think I'm good for now

[hoary light]

This relates to the CFL (Courant Friedrichs Levy) condition, so it is a necessary condition for the numerical solution to be sensible

[still forum]

I was already considering the cfl condition before, I was using something like k/h = 0.25 in the above gif

so i don't think the oscillations are because of that, though if I set k/h > 1 i get blowup

[fathom rain]

https://en.wikipedia.org/wiki/Gibbs_phenomenon

[still forum]

I guess that matches what's going on pretty well

because im trying to approximate a discontinuous initial condition with something continuous?

I'm wondering if theres any trick to reduce/avoid the oscillations near the discontinuities without going out of finite differences

[fathom rain]

artificial viscocity, shock detection, hp-methods

refine the grid close to discontnuity for example

here is an example pic using dg https://www.2pi.se/publications/pdf/p1.pdf fig 5 page 466

[still forum]

im interested in refining the grid, but the discontinuity is travelling around, so does this require refining the grid in different regions as time evolves?

[fathom rain]

yes, it has to be adaptive

but that can be automated

[still forum]

Do you have a good place to read about adaptive mesh refinement?

[fathom rain]

you can probably look at the residual and refine where it is large, here in paper 2 you add artificial viscosity based on that https://uu.diva-portal.org/smash/record.jsf?pid=diva2%3A1805903&dswid=9850

so this paper https://www.sciencedirect.com/science/article/pii/S0021999120308743?via%3Dihub

(just the one i know on the top of my mind since i know all the authors 😛 i can probably find others if you want)

[still forum]

What is a generally good way to plot solutions of a wave equation that are static (so no gifs) in the 1D case?

[hard escarp]

Hi. I'm a phd student in numerical analysis, just getting into the field. Are there good resources on current open problems and research directions within numerical analysis? Specifically numerical analysis for PDE's.

[fathom rain]

i think you need to be a bit more specific than that

[cosmic karma]

any PDEs or specific applications you’re interested in?

[hoary light]

This is a broad question. What aspects of numerical analysis for PDEs are you interested in? For example there are things like accelerating solvers via GPU or parallel computing, focusing on structured linear algebra, linear or nonlinear optimization for implicit solvers, determining good meshing schemes, machine-learning and other data-driven methods for accelerating solutions, regularization of ill-posed problems, developing numerical methods for new kinds of PDE (this is more about modeling than numerical analysis typically), and so much more

[heady mesa]

If you're looking at the algorithms then Spectral Methods is a nice area with lots of new research coming out all of the time.

[shadow spoke]

anyone has this book in pdf form? im looking for 1st edition but all i find is 3rd edition and i dun need that .-.

[fathom rain]

https://www.abebooks.com/servlet/BookDetailsPL?bi=31807991985&searchurl=kn%3Dgilat%2Bnumerical%26sortby%3D17&cm_sp=snippet-_-srp1-_-title1

[shadow spoke]

i want free pdf TwT

[raven dew]

This server isn't an appropriate place to share pirated material or anything

[shadow spoke]

ok

[hardy oar]

hello anyone has pdf for this book: Epperson: “An introduction to numerical methods and analysis”. Wiley. Third edition, 2021

[digital cape]

Read up two messages: this is not the place to ask for pirated books

[hardy oar]

more like ultimate beta

[brave crypt]

Are you part of the SBP community?

[fathom rain]

I know a lot of them (and Jan Nordström was my advisor for a while)

[brave crypt]

Oh I see

Well if i say more i prob gonna doxx myself XD

But i know Nordström too

[fathom rain]

i sit a few offices from Ken Matsson

😄

[brave crypt]

I visited uppsala at some point so we prob alr met during fika

😉

[fathom rain]

i am the wheelchair guy 🙂

[brave crypt]

Will say hi if i have the chance to visit uppsala again

[still forum]

when solving wave equations numerically, I find often that behind/in front of the actual waves, theres some small oscillations. What I want to know is if this is a physical thing that happens or a numerical artifact

and if it's a numerical thing, what are the most common ways to reduce the spurious waves?

[brave crypt]

Are you solving the second order wave equation?

[still forum]

2d spatial wave equation, which was transformed to polar coordinates

[brave crypt]

Then its a numerical artifact

[still forum]

I suspected that was the case

I had some waves that were travelling faster than the allowed speed

[brave crypt]

What is the numerical scheme you used

[still forum]

2nd order central differences for everything

[brave crypt]

I see

Try using upwind schemes

[hard escarp]

@fathom rain @cosmic karma @hoary light Thank you for your responses. Sorry for the delay. I was looking for some less specific strategies for finding those things, actually. For example, one thing I do is I follow arxiv rss feeds. Another thing is I get emails about news SIAM papers.

But, @hoary light according to your list, I'd be more interesting in the last item -> proposing numerical methods for new kinds of PDEs (and studying convergence, stability, etc... )

[hoary light]

I think for looking at new PDEs one is generally approaching the problem as less mathematical and more about engineering and science. In particular, first one needs to ask questions like 'why is this PDE important?' and 'what kind of dynamics are we trying to model?' Once you know what PDE you want to solve then you can think about how to solve it numerically. Otherwise you might get into the situation where you spend a lot of time developing analysis for some PDE that is ultimately not of interest to anybody.

[hard escarp]

So maybe pay attention to engineering journals, maybe? See what kind of equations they're dealing with...

[fathom rain]

are you a phd student, then it is your advisors job to give you appropriate problems until you are capable to form your own interesting questions

[hard escarp]

Right. And that is happening. I just wanted to see alternatives that could be interesting.

[fathom rain]

then probably pick somehing adjacent to that project

[hard escarp]

That is actually a good point.

[fathom rain]

your thesis needs to convey a story and should not be about very differemnt things

[hard escarp]

I actually hadn't thought of that.

And it seems obvious once said explicitly hahaha

Thank you.

Thank you very much. I think I'll start doing study sessions focused on that, actually. It's a really good point.

I don't think I'd be able to think that unless someone told me, to be honest

[cosmic karma]

I second @fathom rain ‘s comment. I started a personal project on Discontinuous Galerkin method for the same dynamical system I’m working on for my thesis. Working on a tangential problem is a great entry point for a new field/method since you already know a lot on the one you’re currently working on

[brave crypt]

Probably an unpopular opinion: it might be useful to look at PDEs theory papers too. Usually in good journals PDEs theorists analyze PDEs that are important.

There are various approaches to develop good numerics for PDEs, and I am on the side that believes good numerics require a solid understanding of the theory, in particular well-posedness and energy estimates.

[gaunt crag]

does anyone know what book to read to become an expert on this flavor of numerical pdes
https://www.sciencedirect.com/science/article/abs/pii/S1290072910002450

or pdf:

[brave crypt]

It seems that the numerical part is only solving the ODEs obtained from the linear analysis

So maybe what you need to find out more is linear stability analysis for PDEs

[gaunt crag]

Do you know of a good reference

[ruby pike]

What's the easiest way to prove that every complex-symmetric is similar to a complex-symmetric tridiagonal matrix?
(Note: NOT Hermitian matrix!!!)

[hoary light]

What's the not-easy proof?

[gaunt crag]

is there a python or matlab alternative of mathematicas ndsolve

[hollow perch]

solve_ivp in scipy?

https://pythonnumericalmethods.berkeley.edu/notebooks/chapter22.06-Python-ODE-Solvers.html

[gaunt crag]

i see

im supposed to write code to get Table 1 on the right but dont know where to start

is the first step to get R? I dont get how they did this

[cosmic karma]

have you looked at what solve_ivp can do to help you?

[brave crypt]

What are the applications of numerical analysis to pure math?

[tiny lake]

Do people do numerical analysis with finite fields (among others, the field F_2^n) or R-modules? Or is it pretty much always R^n or C^n?

[cosmic karma]

Maybe computational group theory? I have a friend who does algebraic coding theory which he uses Magma to either confirm results or even generate them. Basically solving algebraic problems using a computer.

[tiny lake]

Looks interesting, thanks! I'll look into it more :)

[hollow perch]

The complexity of matrix multiplication is highly related to algebraic geometry

[brave crypt]

Numerical analysis lies at the interface between pure and applied maths

But since mostly ppl design and analyze numerical schemes for applied problems, it is perceived as applied maths

[uncut shuttle]

Holomorphic (complex-valued) dynamics, including the Julia and Fatou sets can also be seen as a subset of numerical analysis

So that's where pure maths comes in

[full eagle]

anyone familiar with the Fixed-point iteration method?

[heady mesa]

Yes

[full eagle]

can you explain to me what are the conditions to choose the right g(x)=x?

[heady mesa]

You need a contraction mapping.

[warped fiber]

Could I have help with this please?

[brave crypt]

Hi, its been a while sense i did much linear algebra outside of this class and im having trouble proving or disproving if this second matrix converges or not for Jacobi / gaus seidel methods.
im beginning with jacobi and tried looking at the spectral radius but unsure how to actually compute it

can i get some advice on relevant theorems or ways i should approach this question?

[mighty lotus]

numeric,. idioms.. not cority

[brave crypt]

thanks I think I got it now

[wary oriole]

guys. i have an exam coming up but instructor does not give solution to the practice problem. Here is the problem

Consider the two-dimensional Laplace equation where u is a function u(x,y)
a) Explain how the Jacobi update procedure can be derived from this equation (explain
how the equation is discretized, how are the derivatives discretized?
b) How can you improve convergence of this scheme (describe update scheme, stencil,
memory requirements) ?

i know part (a) is

but how can you answer part (b)? I would like a precise answer to get full points if it appear on the actually exam

oh i think i get it. "improve convergence of this scheme" is wrong wording i think/ just stencil and how to storage solution u

[vapid plume]

@wary oriole Do you know about the gauss seidel method?

[wary oriole]

yeah what's up

[still forum]

since they asked about memory requirements, think about which one is better of gauss seidel and jacobi

[brave crypt]

any textbook recommendations on the topic For beginners ?

[vapid plume]

Numerical methods and analysis? James F. Epperson

[gaunt crag]

I don't know how to solve for eta when I dont know RF

44-47, 48 is my problem

first step is solve for eta and RF with shooting method and then runge kutta the rest

[dusk meadow]

If I have a pair of intersecting k-dim planes in n-dim space, represented as y = Ai@x = bi, where each Ai is a full-rank n-by-k matrix, what is a way I can find a representation for their intersection? That is, y = C@x + d.

I came across the Zassenhaus algorithm, however this is for subspaces.

[plucky kayak]

your set of equations is equivalent to one large system where matrix was obtained by stacking A_i on top of each other

so you just need to pick any particular solution of that system and add any element of null space of the stacked matrix to that, that describes the entire intersection

[opal mason]

can you elaborate on the step of adding an element of the null space?

[plucky kayak]

the solution set to Ax=b is x_0 + z, where x_0 has the property that Ax_0 = b and z is in the null space, i. e. Az=0

[opal mason]

ohoh ok I see

[opal ocean]

sup, is anyone well versed in the finite element method? I'm working on an implementation for a project but so far it hasn't been going so well, so it would be nice to find someone who knows the stuff

[grave spoke]

what kind of pde?

[opal ocean]

schrodinger equation in a two dimensional well

so rn just finding eigenfunctions for the Laplacian

[brave crypt]

discontinuous Galerkin or standard continuous Galerkin?

[grave spoke]

using linear triangular finite elements? what have you done already?

[opal ocean]

standard continuous

[opal ocean]

and yep, triangular elements

well, so far I have code to generate meshes of the shapes I want, and code to compute the stiffness and mass matrix

the data structure for the meshes is kinda dumb, but the code has to be read by people who aren't into programming

[grave spoke]

so, what's left to do? and what hasn't been going well?

[opal ocean]

if you have any questions regarding the code, please ask!

[opal ocean]

well I have my eigenvalue problem

but all my eigenvectors give me nonsensical solutions

[opal ocean]

this one

so

to be more precise

I'm unsure of how I should deal with the edge condition (i.e., making the edge test functions zero)

what ive done is make the corresponding rows zero and the diagonal values on these rows 1

and I get solutions that are zero on every node, but have a huge singularity on one node

kinda like a dirac

like that

keep in mind that the solution I'm looking for is some kind of trigonometric polynomial

[opal ocean]

oh also im sorry, this is old code, there are a couple of coefficient mistakes there that i've fixed, however, I still get the same spiky function problem

[grave spoke]

eigenfunctions for the laplace operator work fine, but i think there is a mistake with the mass matrix as it is not always symmetric

[opal ocean]

wait what?

did you run my code?

[grave spoke]

yes

well

modified

[opal ocean]

okay yeah that makes sense

also its very strange that the mass matrix isnt symmetrical

since I only fill the upper triangular half, then make it symmetrical

[grave spoke]

maybe you can send me the fixed version, does it change something about the mass matrix?

[opal ocean]

only the diagonal coefficients

M = M + np.transpose(M)```

line 153

also I can't thank you enough for giving me some of your time

[grave spoke]

for zero boundary conditions, you can just remove the rows and columns corresponding to the nodes at the boundary btw

[opal ocean]

it's incredibly kind of you

[opal ocean]

I thought so, but doing so gives me a non invertible matrix

since it has null rows (if null is the right word)

[grave spoke]

and use np.linalg.eigh (for Hermitian matrices) to get eigenvectors ordered by the magnitude of eigenvalues (np.linalg.eig will not necessarily order the eigenvalues and eigenvectors as they may be complex-valued)

[opal ocean]

oh you mean remove them remove them

[grave spoke]

remove them as in also reducing the size of the system

[opal ocean]

but if I remove them, isnt it the same as considering the problem on a smaller inner square?

[opal ocean]

so did you mean

something like that?

#   find edge nodes
    to_remove = []
    for i, node in enumerate(list(test.keys())):
        if node[0] == 0 or node[1] == 0: # check edge nodes
            to_remove.append(i)
    #remove corresponding rows and columns
    K = np.delete(K, to_remove, 0)
    K = np.delete(K, to_remove, 1)
    M = np.delete(M, to_remove, 0)
    M = np.delete(M, to_remove, 1)
    P = np.linalg.inv(M)*K
    eigenvalues = np.linalg.eigh(P)
    # fill back the vectors by adding zeroes to edge values
    eigenfunctions = []
    for i in range(len(eigenvalues[1])):
        vect = list(eigenvalues[1][i])
        for j in to_remove:
            vect.insert(j, 0)
        eigenfunctions.append(np.array(vect))```

[grave spoke]

sorry, i had other things to tend to

[grave spoke]

not quite, consider the global stiffness matrix resulting from three one-dimensional elements of the same size, for example

\[K = \arr{1 & -1 & 0 & 0 \\ -1 & 2 & -1 & 0 \\ 0 & -1 & 2 & -1 \\ 0 & 0 & -1 & 1}\]

[pine jetty]

戯言 MSC2020 49M41

[grave spoke]

the above matrix is singular

but if we remove the rows and columns corresponding to the boundary nodes, we get an invertible matrix

\[\tilde{K} = \arr{2 & -1 \\ -1 & 2}\]

[pine jetty]

戯言 MSC2020 49M41

[grave spoke]

note that you want P to be a symmetric matrix

what you can do to reduce the generalized symmetric eigenvalue problem to an ordinary symmetric one is to consider the Cholesky decomposition of M

\begin{gather*}
M = L L^* \\
K \psi = E M \psi \Leftrightarrow L^{-1} K (L^*)^{-1} \psi = E \psi
\end{gather*}

[pine jetty]

戯言 MSC2020 49M41

[grave spoke]

so, you define (P = L^{-1} K (L^*)^{-1}) instead of (P = M^{-1} K) to retain the symmetry

[pine jetty]

戯言 MSC2020 49M41

[grave spoke]

the psi on the right side of the equivalence should be L* psi (so you need to transform the eigenvectors back)

[fathom rain]

eigenvalues of K is just $f(\theta_j)=2-2\cos(\theta_j), \theta_j=(j-1)\pi/n$

[pine jetty]

sverek

[fathom rain]

eigenvectors also known in closed form

[grave spoke]

we're talking about a general finite element mesh

[grave spoke]

in 2D actually, this is only an example to illustrate the boundary conditions

[grave spoke]

    test = meshes.general_mesh(meshes.cube_nodes(1, 1, N=16))
    matrices = browse_mesh(test)
    K = matrices[0]
    M = matrices[1]
    n = np.size(K, 0)

    bnd = []
    for i, node in enumerate(list(test.keys())):
        if node[0] == 0 or node[1] == 0 or node[0] == 1 or node[1] == 1: # check edge nodes
            bnd.append(i)
    K = np.delete(K, bnd, 0)
    K = np.delete(K, bnd, 1)
    M = np.delete(M, bnd, 0)
    M = np.delete(M, bnd, 1)
    L = np.linalg.cholesky(M)
    Linv = np.linalg.inv(L)
    P = Linv @ K @ np.transpose(Linv)
    eigenvalues, eigenvectors = np.linalg.eigh(P)

    dof = [i for i in range(n) if i not in bnd]
    for i in range(4):
        sol = np.transpose(Linv) @ eigenvectors[:, i]
        vec = np.zeros(n)
        vec[dof] = sol / np.linalg.norm(sol)
        display_vector(test, vec)

this seems to work well using the fixed version @opal ocean

for the mixed zero Dirichlet and Neumann boundary conditions, you get the following discrete eigenfunctions

[opal ocean]

oh that's definitely alright, again i'm very grateful

thank you so much, that's very kind of you, ive been stuck on this for a while

Im unsure of why P should be symmetrical, is that a thing about generalised eigenvalue problems?

I havent come across that idea anywhere, so if you had resources you could point me to that would be great

[grave spoke]

If you want to use eigh, then P should be Hermitian. You can still use M^{-1} K (the matrix will be diagonalizable and will have the same real eigenvalues), but you won't be able to use the algorithm specialized for Hermitian matrices (which also orders the eigenvalues and eigenvectors in numpy).

[grave spoke]

https://www.netlib.org/lapack/lug_old/node56.html

[opal ocean]

wait so they were the right eigenvalues, just not in the right order?

[opal ocean]

thank you

okay so what i did was right, I just had the vectors in a random order, and probably fucked up the reconsitution of the vector with value 0 on the boundary

I really can't thank you enough, you've been very patient and nice to me

[grave spoke]

I have checked; apart from the boundary conditions, you should write P = np.linalg.inv(M) @ K instead of P = np.linalg.inv(M) * K.

In numpy, * is elementwise multiplication, while @ is matrix multiplication

[grave spoke]

It's even better to use P = np.linalg.solve(M, K) for efficiency

[valid iron]

Yo, my teacher gave us this homework to prepare our exam but i don t understand what he is asking. He gave us a code we have to complete the iterative refinement but i don t understand how, and what we need to use. why factorisation LU ? how to prevent cancellation error for the line r = b - Ax ?

Here’s the task
The project file Lefichierprojet2024INFOF205.m contains a simulation of a linear system A⃗x=⃗b, where the matrix A is ill-conditioned. The simulation first generates a diagonal matrix Σ with a specified condition number and two arbitrary orthogonal matrices U and V, to construct A = UΣVT. Then, the theoretical solution of the system is generated as a vector of integers. The values of ⃗b are found by executing the operation ⃗b=A⃗x. You are asked to complete the routine by recalculating ⃗x from A and ⃗b, and adding a number of iterations nit = 20; of iterative refinement (using a for loop i=1:nit instead of while norm(e)/norm(x) > tol). Use explicit LU factorization in your implementation (command [L,U] = lu(A);). Avoid the implementation from slides that suffers from cancellation errors in the line r = b - A*x. After each refinement, store the values of the residual and the error ehat (i.e., the solution of the system, not the exact error). Default values for n and condA can be used (n=10 and condA=1.e14). Include the completed routine in your folder.

[wooden kiln]

Is it possible to prescribe a Neumman Boundary Condition for a variable, if you have their discrete values in Spectral Space?

[dawn wadi]

Hey, I’m having trouble finding resources on the appropriate definition for a stable algorithm, the course book “scientific computing” by Michael T. Heath gives vague supposedly equivalent definitions using casual wording without definition, like “nearby problem” etc.
Without any examples too, all of this is just a small paragraph btw.

This is troublesome to me because one of the exercises asks to give a proof for why addition is a stable algorithm under the IEEE standard’s machine epsilon.

Thanks!

[next garden]

I think the FEM finally clicked for me:

[worn pivot]

oh

[brave crypt]

Maybe what you want to look for is backward stability

[dawn wadi]

Looking at the wiki definition, can I rephrase my exercise formally as

[pine jetty]

Aslan

[dawn wadi]

Here F is meant to represent the set of floating point numbers

[dawn wadi]

So since fl(x+y) = (x+y)(1+eps) = x(1+eps) + y(1+eps).

Can I just pick x^ = x(1+eps), y^ = y(1+eps) and we’re done?

[next garden]

There's something I don;t understand with the FEM.

say our equation is u'' = f

We do the weak formulation and get

$\int u'' w dx = \int f w dx$ And from here use integration by parts to get something in terms of just u'w'

Then we replace u by a piecewise linear approximation $u \approx \sum u_i v_i$

However, the second you make that approximation, u'' would yield 0 as a function. So the equation:

$\int u'' w dx = -\int u' w' dx$

That you get from itnegration by parts should be false. So how come the equality holds? The left hand side and the RHS should evaluate to different things in my head.

[pine jetty]

Makogan

[grave spoke]

What do you mean by u'' would yield zero as a function?

Oh, you mean piecewise in the interior of each subinterval, yeah

[next garden]

yes, your tent function should be 0 when taking two derivatives

so my head doesn;t understand how the equality holds under the approximation

[grave spoke]

It doesn't. The piecewise linear approximation is not twice differentiable at the breakpoints in the classical sense (the integral of u'' over the whole domain is not defined).

[next garden]

That doesn;t really remove my confusion though. Before the approximation we have an integral equation.

Then we make the assumption that $u$ admits a given approximation, but the approximation we got breaks the analytic equation. And yet we proceed to integrate the RHS as if all was fine.

[pine jetty]

Makogan

[next garden]

Why does this work in practice?

My head keeps telling me it shouldn;t because the fundamental assumption you made through the integration by aprts doesn;t hold

[grave spoke]

Well, that's the idea, the analytic equation only holds for functions that are twice differentiable in the classical sense, so you weaken it by considering the weak formulation instead, which is still defined for H^1 functions.

[next garden]

But how does the solution you find manage to approximate the problem

When it won;t have a "correct" second derivative?

basically I don;t understand how breaking the fundamental constraint from your equation allows you to nonetheless find something useful

[brave crypt]

Depends on how the exercise is written. But usually x and y are not floating point numbers. So what we want to investigate is fl(fl(x) + fl(y))

[silk drum]

It's a valid conceptual question. $u''v = (u'v)'-u' v'$ so, by the equality, we can claim a lot of the meaning and information on the LHS is completely (in the sense of equality) captured by the RHS. But obviously we see here the LHS has a second derivative whereas the RHS does not.

FEM doesn't know or care that you started with an equation with a second derivative so it would be a dig into the convergence theorems and error estimates of how well FEM would approximate the true solution with different error estimates. I would guess any error estimate including second derivatives would be tragically bad, as pointed out.

Ultimately, the equality holds by the product rule for the derivative as used in integration by parts and zero is a perfectly valid second derivative in the interior. It's actually on the element boundaries that you have a problem of non-existence of some derivatives, if you only enforce continuity.

To work through the conceptual hurdles, one could apply FEM first to $u=f$ (you just get a mass matrix) and use higher order elements with non-zero interior second derivatives to get a feel for things.

[pine jetty]

mikeliuk

[next garden]

Thank you for the help. I am rumminating over this.

[dawn wadi]

Yes x, y are not floating point numbers as ive written in my case.

Why would we want to investigate your expression?

According to def. of backward stability (the thing you mentioned) we went to investigate fl(x+y) and if it’s equal to an “exact” perturbed solution of the corresponding operation (addition), that is something like x^ + y^, as I’ve written above.

[pine jetty]

Aslan

[dawn wadi]

Or am I missing something?

[dawn wadi]

This is more thorough I think?

[brave crypt]

My interpretation of fl is a function that maps a real number to a floating point number

Since x and y are real numbers, to do a floating point addition we need to convert it to floating point numbers first

So, when you say x^, thats already fl(x) since x^ = x(1 + eps) = fl(x)

And then one needs to take into account the fact that adding two floating point numbers may not result in a floating point number (since the number of bits required may be increased)

So its fl(fl(x) + fl(y))

[dawn wadi]

Yeah, where am I not saying this?That’s why I asked if picking them like so is enough to show backward stability in the first reply above.

[dawn wadi]

I see, but isn’t this a consequence of the algorithm. That is the addition between floating point numbers is not the same type of addition between real numbers(?) like isn’t this already implied?

https://en.wikipedia.org/wiki/Numerical_stability#Stability_in_numerical_linear_algebra

Here’s the definition I’m using for backwards stability btw

[dawn wadi]

But say we use this, what would I exactly be investigating as you suggested?

That is, what is enough to show that the addition algorithm is stable using that expression?

[brave crypt]

$ fl(fl(x) + fl(y)) = (x(1 + \varepsilon_1) + y(1 + \varepsilon_2))(1 + \varepsilon_3) $

,tex so $ (1 + \delta)(x + y) = (x(1 + \varepsilon_1) + y(1 + \varepsilon_2))(1 + \varepsilon_3) $

[pine jetty]

Kamiya

[brave crypt]

you are missing \varepsilon_3, which is the rounding error after the addition

[dawn wadi]

I'm not really following, what does this exactly say?

Like so what? Why does this imply backward stability?

[dawn wadi]

Assume for simplicity that x and y are exactly representable in my eq. then that epsilon is the error after addition.

[brave crypt]

,tex the $ |\delta| $ is what you want to bound to get backward stability

[pine jetty]

Kamiya

[dawn wadi]

So long as its bounded the algo. is stable?

[brave crypt]

If they are exactly representable, then they are already floating point numbers

[brave crypt]

,tex it has to be $|\delta| = O(\varepsilon)$

[pine jetty]

Kamiya

[dawn wadi]

I mean I'm just saying that bcuz, i thought it was implied that fl(x+y) = x+y(1 + eps) was enough, that's atleast the standard model they showed us; i guess its just in a rewritten form.

[brave crypt]

Yea fl(x + y) is enough if x and y are already floating point numbers

[pine jetty]

Aslan

Aslan

[dawn wadi]

My point is, even if x, y are not already that, can't we just rewrite it as such?

So it seems reasonable why we dont bother with that (in this case)

[dawn wadi]

I see

[dawn wadi]

bounded up to norm*

[brave crypt]

,tex Well lets say for simplicity we assume x and y are already exactly representable (floating point), then $(1 + \delta)(x + y) = (1 + \varepsilon)(x + y)$, so $\delta = \varepsilon$, and then we are done.

[pine jetty]

Kamiya

[dawn wadi]

Yeah i was sorta half joking but realized now that it doesnt even matter it seems, since we can just rewrite (up to error below machine eps.) in the same form anyways

[dawn wadi]

See above

Like in our case since we're still below any machine eps. error, isnt it fine to write fl(x+y) = (x+y)(1 + eps) even if x,y are real?

[dawn wadi]

Since we can always just rewrite that error in terms of some other error say eps, and have it in a form like this ?

Im sorry if im being annoying lmao

[brave crypt]

,tex so we have $\delta = \varepsilon_3 + \frac{\varepsilon_1 x + \varepsilon_2 y}{x + y} + O(\varepsilon^2) $, yes in this case it will lead to the same conclusion

[pine jetty]

Kamiya

[brave crypt]

but if the question is x - y, you will find that the longer method leads to a different conclusion

[dawn wadi]

Thats true

[brave crypt]

,tex since its not stable when $ x \approx y $

[pine jetty]

Kamiya

[dawn wadi]

Thats intersting, for some reason our textbook says it fine to write it for minus the same way

Im only in undegrad so maybe this is something subtle they dont think was important in general?

[brave crypt]

its actually a well-known phenomenon

https://en.wikipedia.org/wiki/Catastrophic_cancellation

[dawn wadi]

Yeah was about to say that they instead talk about this separately but without any "rigor", but before this they just say yeah you can write fl(x-y) = (x-y)(1 + eps)

Let me double check

[brave crypt]

yea if they are exactly representable, nothing bad will happen for sure

[dawn wadi]

Okay i was just dumb, it seems its implied that x and y are floating point numbers already

[brave crypt]

frankly, we dont need to worry about this since most scientific prog languages are quite robust alr

but we do need to care about it when it comes to linear systems (especially the ill-conditioned ones)

[dawn wadi]

I see

I think i also realized that in the exercise x and y are not real numbers but also just floating point numbers 💀

I should of just written it down, never ever making this mistake again

[brave crypt]

awesome. then this suffices.

[dawn wadi]

thank you for wasting your time. I learnt alot from this

[brave crypt]

No worries

[wicked pecan]

Are there any coefficients for 8th (or higher) order RK methods (for stiff or non-stiff) that are generally considered good?
I want to try out some higher order methods for a project of mine and whenever you look for them online you find research papers with large butcher tableaus that require that you choose certain parameters yourself, and I neither know why I'd use certian coefficients over other coefficients (assuming both of them were published in research papers), or how I'd choose the variable parameters for those butcher tableaus where you need to choose certain paramaters yourself...

[grave spoke]

My advisor once told me that back in the days, determining the coefficients of some higher-order RK method might be worth a PhD

[opal ocean]

thank you! (that's what I ended up going with)

not the solve one

just multiplication

now i'm trying to implement a 3d version but jesus christ the coefficients are a nightmare

[hybrid mulch]

Hello, I have an assignment that I need help with , I posted it in help 5. This is really important for me passing the class. If there is someone who is good with analysis and has some extra time i would pay you to help me out.

[fathom rain]

paying somebody to do your homework is cheating

[hybrid mulch]

Not to do my homework

To tutor me

and teach me how to do it

[wicked pecan]

Well, did anyone claim that phd and publish some nice coefficients for me to use? :)

[grave spoke]

Sure, here's the historical overview:

The thesis mentioned in the overview for an 8th order RK method:

And the coefficients:

[wicked pecan]

Oooh

Thank you so much

That’s exactly what I needed

[gaunt crag]

F(y,h), h, and a are specified. R is unknown eigenvalue. how do i find eta via shooting method without knowing R?

[golden wolf]

Hi! Is anyone familiar with 'prolate spheroidal wave functions' or 'Slepian functions'? I am mainly looking for a use case for them and to understand their real applications. It seems like a niche topic in that I can't find any applications for it, but it also seems highly theoretically useful.

[heady mesa]

Spectral Methods or Finite Difference methods? Which are more powerful?

So far I've only studied Spectral Methods and I mean they are really OP but was wondering if it's worth looking into Finite Difference methods as well.

[brave crypt]

It depends on what you are trying to solve

[heady mesa]

Non linear ODEs

[brave crypt]

Spectral methods are not that great for stiff ODEs.

For PDEs, spectral methods are really bad for nonlinear conservation laws.

This is due to the Gibbs phenomenon.

[heady mesa]

Interesting

I have been using them to solve Non-linear ODEs such as the Lane-Emden equation

I'll look into this stuff more, cheers

[brave crypt]

The Lane-Emden equation is essentially a Poisson equation with a power law term.

My gut feeling says that the weak solution is in H^2 (or maybe H^1).

[brave crypt]

By the Sobolev embedding theorem, in 1D, this implies that the weak solution is continuous.

So its not surprising if you find spectral methods powerful for this problem.

[heady mesa]

Nice, yeah I mean it's defined the semi infinite domain and hence I've been truncating the domain when solving using a computer

Combined it with the quaslinearization method

[heady mesa]

Oh cool, so finding the natural function space of the solution as a PDE equivalent allows you to extrapolate information about the ODEs solution?

[brave crypt]

Well, a 1D Poisson equation (PDE) is an ODE.

[heady mesa]

Yeah exactly

That's cool, really helpful, I don't think I'll be solving any stiff ODEs as part of this project.

[brave crypt]

My research is computational but we develop numerical methods by starting from the mathematical theory of PDEs

So I always try to understand the problem im trying to solve at the continuous level first.

[brave crypt]

This should lead to good results. Since the ODE only admits continuous solutions, linear stability suffices.

[heady mesa]

That's cool and makes sense.

At my university we haven't done analysis on ODEs, only Greens function and some basics on distribution theory. But this year I've done a two term module on PDEs and it feels like the analysis in this area is more rich (surprisingly).

[heady mesa]

Have you used the quaslinearization method (QLM)?

My supervisor had never heard of it...

[molten knot]

anyone have any ideas for a good 45 minute numerical linear algebra/optimization presentation topic? we’ve gone over all the basic factorizations and descent methods in the course but can pretty much choose any related topic for our project

[odd minnow]

Does anyone have learning resources in the form of videos that use the Advanced Engineering Mathematics book as a reference? Can you share the channel?

or Do you have a answer key with the the solution of Advance engineering mathematic? i want to do some exercise

[heady mesa]

The QR algorithm for finding the roots of a polynomial is cool

Basically you turn a polynomial into a matrix (called a companion matrix)
where the eigenvalues of this matrix are the roots of the polynomial. Then you use QR factorisation to approximate the eigenvalues (which are the roots of the polynomial)

[molten knot]

i think i saw something like this for the dynamic mode decomposition if you’re familiar with that, but i think they use QR algorithm to find the least squares solution to coefficients in the last column of a companion matrix

this topic sounds cool though i’ll check it out thanks!

[brave crypt]

I haven't used it computationally.

Because its too expensive.

[brave crypt]

If you prefer applied projects, maybe image denoising is also a good area to look at.

Something like total variation denoising, it involves both linear algebra and optimization.

[errant wagon]

So there are multiple textbooks named Advanced Engineering Mathematics. There are solutions manual available for Dennis Zill 7th edition online (however cost money). There is also quizlet.com which has fully worked solutions however again you will need to pay for access.

[heady mesa]

Interesting, I mean my programmes for solving to super high accuracy with this method takes just under 1 min, is that too slow?

What method is faster?

[brave crypt]

Oh because I mainly solve PDEs.

Its very time consuming for PDEs.

[heady mesa]

Ah okay, makes sense 😂

[molten knot]

this actually seems really cool thanks

[shut vine]

Hi, does anyone have any idea of deriving the explicit formula of X{k+1} for the classiscal Kaczmarz alg? Thanks!

[rancid light]

yeah, we can start from the iterative update equation

x^(k+1) = x^k + (b_i - ⟨a_i, x^k⟩) / ∥a_i∥^2 * a_i

where i = (k mod m) + 1, a_i is the i-th row of matrix A, b_i is the i-th element of vector b, and ∥a_i∥ denotes the Euclidean norm of a_i

you can expand that

[boreal glade]

Now that is weird

[heady mesa]

What is weird about it?

[tiny lake]

What does the current research in numerical analysis focus on?

Fine-tuning Krylov subspace methods? Less precise, small-width floating-point algorithms for machine learning? Something else?

[heady mesa]

Lots of areas in numerical analysis

Better algorithms for solving ODEs and PDEs. Better algorithms for approximating eigenvalues of certain matrices ext.

[fathom rain]

randomized methods, matrix functions, mixed precision, quantum algorithms, new methods for fde/pde/ode, parallel-in-time methods, preconditioners/multigrid, inverse problems, new eigenvalue-eigenvector methods and analysis tools (just a few examples, just look at topics of random conferences)

[random plume]

parallel, streaming and distributed methods, adapting some methods to real-time applications, numerical approximations of various infinite-dimensional operators (using pseudospectrum to verify calculations)

anything that has to do with tensor decomposition

And I second the part about better eigenvalue decomposition. The singular value decomposition is heavily used, and I don't know (yet) of a good numerical algorithm to obtain the jacobian normal form decomposition in a way that's not known to be numerically very unstable

even though there's a lot of cases where applying the singular value decomposition will make you loose information about generalized eigenvalues/vectors (there's cases where the singular value decomposition is not the decomposition you'd want to have)

[cursive yew]

This is from "What every computer scientist should know about floating-point arithmetic" by Goldberg. Can someone explain how the β' is the digit β/2

how did he get the error in the next line

[north marlin]

need help with my problem for crank nicolson implementation of black-scholes transformed heat equation https://colab.research.google.com/drive/15vN9621QVy36_prErLeX40tvOWRXdybs?authuser=2#scrollTo=BfFRdeysDdL- i cant seem to properly transform my variables (commented out code at the bottom under function crank_nicolson_heat_matrix() and I could use some help

[spark falcon]

What would you consider the best algorithm to calculate SVD? The best in terms of how good it is in image compression compared to standard library implementation or something simple like power method

[plucky kayak]

analytical continuation is a complex analysis topic

then you should fully explain your problem and what you want to do

because under certain conditions there's only one analytic continuation so it doesn't make sense to talk about having multiple of them

ah I guess you're talking about analytic continuation of f_e, which is determined for each set of data

well the answers to first and second question is "it depends on the problem"

choice of norm would be different and sometimes greatly affects the solution, which may or may not be beneficial in a particular problem

you may also find f_e(x) that equals f_t(x) for every x in the data set, yet it's way off from true function at other x's

https://en.wikipedia.org/wiki/Runge's_phenomenon

I think uniform convergence should take care of this problem, but it's probably very hard to justify it outside of special cases

I guess gradient boosting could be the closest thing to what you're looking for that actually works

it also produces a sequence of approximating functions that are chosen in some optimal sense

[coarse field]

Hello guys!!! I have a question , how can i aproximate the solution of an equation like h(x) = 0

and lets say that h(x) has one zero in a closed interval [a,b] , h(c) = 0

and h'(c) also equals zero

[cosmic karma]

if you don’t need a really fast one, maybe bisection method? haha

or secant method is probably so much better

[coarse field]

i think i cant use it here

lets suposse the h(x) is something like this

[heady mesa]

Yea secant and bisection wouldn't work, Newton Raphson would have floating point problems but other than that is should work?

[heady mesa]

On the example you showed in the graph

[heady mesa]

You could find the root of $h'(x)$ instead

[pine jetty]

Max

[coarse field]

idk

[coarse field]

but suppose that h'(x) is also something like this

ahhh no it cant

yeah!!!! good idea

tysm!!!

[coarse field]

how can i implement sympson rule ?

where the function recieves the error the user wants , and the function ?

[fathom rain]

assume the function is nice enough, estimate the error using Richardson, and then find what approximate n is needed for the error tolerance?

[coarse field]

I didnt learn "Richardson" can you try to explain me what is that ?

My issue with this is that to find the number of sub-intervals, it's necessary to bound the 4th derivative of the function. I don't think there's a general formula to do this for an arbitrary function, so I'm not sure how to proceed. Does this 'Richardson' you mentioned have anything to do with this?

[fathom rain]

Error approximation will be $(S(h)-S(2h))/15$ since it is a 4th order method ($(S(h)-S(2h))/(2^p-1)$). Then you have that that error is a[proximately $Ch^4$ compute $C$ and pick an $h$ so error is smaller than your tolerance

[pine jetty]

sverek

[fathom rain]

where $S(h)$ is simpson of ypur function for some $h$

[pine jetty]

sverek

[dapper crystal]

$\int_{0}^{1}e^{-x^2}dx$

[pine jetty]

Ångst

[dapper crystal]

i'm trying to calculate this integration with simpson rule and i got stuck

[fathom rain]

what was the problem

[cursive yew]

does anyone know any good resources to study about rounding error analysis ?

I have this book called Rounding error in algebraic process by Wilkinson. The book discusses the errors in arithmetic operations, computations involving polynomial and matrix operation in great detail. But it was published in the 60s. So i am not sure if it is updated enough for one to study. I write my codes in C

[dapper crystal]

thank you for the response, problem was related to my basic integration knowledge, after reviewing some notes I've succesfully calculated that integration

[heady mesa]

Might easier to search by looking at "floating point arithmetic"

[cursive yew]

Well, I have been reading this article by Goldberg: What every computer scientist must know about floating point arithmetic. It gives a detailed study of floating point arithmetic, "special" floats like NaN, infinity etc. I would like to see examples cases being examined

which i haven't found as now

good examples

[heady mesa]

The quadratic formula is the classic example

$\frac{-b\pm\sqrt{b^2-4ac}}{2a}$
\
\
If $a\approx0$ we have problems

[pine jetty]

Max

[heady mesa]

So times top and bottom by the conjugate of the top

And you get $=\frac{2c}{-b\pm\sqrt{b^2-4ac}}$ which no longer has that problem

[pine jetty]

Max

[cursive yew]

Yeah this one's mentioned in the article I mentioned above

[heady mesa]

I mean, personally I used this rounding stuff to improve my differential equation solving

essentially you can solve differential equations using polynomials, higher power the polynomials often the better the solutiopn

but higher the power polynomial, more opportunity for floating point errors (especially for large x)

So I would "patch" the domain and use lots of low order polynomials

[cursive yew]

I see

I am just trying to learn a more general "hand-written" approach for finding bounds for rounding error since I am planning to make available "unified" notes, codes and demonstrations. "Unified" as in that I have found that many resources in numerical analysis skip over things (like a proper rigorous error analysis, ofc they mention it once in the beginning of the course, but later seem to forget it) and i wish to make a website (completely free and codes under GNU general public license, meaning that they can run, study, share and modify the codes) where anyone can come and not only start from the very basics, but also learn the subject in more depth for a better understanding of ideas. Standard books or courses are nice but I think all of them have some gaps needed to be filled. I will try to fill them as best as I can

[heady mesa]

That's pretty cool

[heady mesa]

Yeah I think because in general you just use more machine precision, such as by using mpmath in Python

Like if I feel the result is wrong I just make the precision insanely big

[cursive yew]

Yeah but I'm more into using C for my notes

[heady mesa]

There will be an equivalent library for C

[cursive yew]

Well, not going to use any library other than the standard ones

Going to write everything from scratch

[fathom rain]

what do you mean "standard"?

mpfr and arb are very standard

[cursive yew]

Whatever comes pre-installed

[fathom rain]

preinstalled in what?

most machines have no compilers installed...

[cursive yew]

linux usually has gcc and some libraries installed

Idk about windows

[fathom rain]

so "standard" for you is whatever linux distribution you chose for the moment. ok. that is quite non-standard definition...good luck

[cursive yew]

Well, wouldn't be using anything other than stdio, stdlib, math i think

[autumn sphinx]

chants cuda cuda cuda cuda

[tiny lake]

Probably being part of the standard library specified by the language standard

[heady mesa]

@fathom rain do you use C alot?

I use mainly python, but are there advantages to the other languages like C in numerical analysis?

[cursive yew]

Well, first and foremost, C is really frickin fast.

[fathom rain]

No, I never use c unless I have to. For research I use like 95% julia

[fathom rain]

Use the tool that is best for what you and your collaborators do

[cursive yew]

C can take some time for people to get comfortable with but I think once you get used to it, your code can do stuff really faster. Take this example (taken from mit 6.172 lecture 1): Matrix multiplication. Two NxN matrices, whose elements you generate randomly, and then multiply the two matrices using $c_{ij} = \sum_{k=1}^Na_{ik}b_{kj}$.
For N=4096:
Python takes about 6 hours.
Java takes 46 minutes
C takes 19 minutes

[cursive yew]

Yes, that is eventually what you should do. Keeping in mind what your collaborators use

C is more like a personal preference for me

C can be intimidating at first tho. It lacks several features that python does, like garbage collection. But again, once you get used to C, it should be easy

[pine jetty]

notssoshocking

[brave crypt]

Did that example use numba?

Otherwise it is an unfair comparison, since Java and C are compiled.

[cursive yew]

Nope

Only these 3

[heady mesa]

Interesting, I'll have to look at Julia, haven't heard of it before

[heady mesa]

Damn

Yeah I may have to learn a bit of C then

[cursive yew]

Not necessarily, you could go with Julia as well

C is more like a personal preference for me

People have been using C for decades but nowadays you have languages like Julia if you want code simplicity and speed both

If all you wish to do is numerical analysis without caring much about the code part, go with julia

[brave crypt]

Yea what I meant is they should have used numba in python

Since numba makes your code compiled

If they just used standard python, the result would be obvious without even testing it

There is no way an interpreted language can run faster than a compiled language

[cursive yew]

So, I found this. https://gist.github.com/jfpuget/6f8c82b729677b0173d0

This page uses LU decomposition as a test case. I didn't look at the codes in detail but skimmed through the page for now.

C seems to have performed well as compared to Numba for higher N. Numba performs really well as compared to C till around N=15. After that C catches up and becomes faster. Similar case for C vs cython, though cython couldn't outperform Numba. Somehow scipy and lapack did better than C for higher N.

Julia performed really well as compared to C for lower N. Even for higher N, they were comparable

I haven't yet looked into the validity of this page though

[tiny lake]

Honestly, I see no benefit in using C unless your target platform isn't supported by existing compilers

But I've heard Fortran can be pretty good even these days because of its very aggressive optimisations

[cursive yew]

Well, as I have said, it's more of a personal preference

[brave crypt]

Obvious. Julia is JIT compiled.

[wet magnet]

hi, I'm wondering if anyone can give me some suggestion on how to get started on modeling pdes. Especially I'm reading this paper https://arxiv.org/pdf/1802.06089 and it has a lot of plots, and I'm not sure how to get those. How do you get a numerical simulation of a equilibrium solution? They just generated some point particles randomly but the PDE only involves the population density, how does the population density affect a single particle? I guess I'm not sure how to model that

sorry if I phrased my question poorly, I'm just really confused about how to get started. If anyone could point me to useful resources so I could learn similar things I would be very grateful

[heady mesa]

Depends what method you want to use to solve PDEs, do you know what type of PDE you want to solve? @wet magnet

[wet magnet]

im just looking at the one in this paper

[wet magnet]

and the riemannian manifold we want to solve it on are the 2d sphere and the hyperbolic plane.

[heady mesa]

(I can't suggest a good method as I don't do PDEs)

Id suggest looking up this PDE and high precision solution or "solving ..."

Then read some papers around this area to find out the state of the art methods.

Then try to learn them methods first from a textbook or thesis

[wet magnet]

oki, i see, i dont really want the details tbh im kind of short on time, i guess maybe I should skip the numerical parts ... but anyway, ty for your suggestions

[heady mesa]

Some methods will be quick to learn, like basic finite difference schemes or basic spectral methods

But if you want high precision then that will take some time to learn about

[wet magnet]

i see, ill look it up

[grave spoke]

look up meshfree methods

[wet magnet]

oh, this looks cool, where is this from

[grave spoke]

and Runge-Kutta methods for the integration of the resulting system of ODEs

[grave spoke]

the paper you sent

[wet magnet]

lol i only saw familiar words and got excited

didnt really find anything on numerical results of this pde on general manifolds

[grave spoke]

numerical analysis on manifolds is still a relatively new field

[wet magnet]

like this?https://scholarworks.boisestate.edu/cgi/viewcontent.cgi?article=3156&context=td

[grave spoke]

yes, topic-wise

i skimmed the document and i don't think it's a good introduction to it

[wet magnet]

what would be a good intro

like a youtube video?

[grave spoke]

a course on it would be best

or pick up a book about meshfree and particle methods rather than a dissertation

[manic pollen]

hey there, does anyone know the difference between numerical analysis and numerical methods? are both of them the same? or if there is a slight difference, what would that be?

[heady mesa]

The key difference is the word "methods" and "analysis"

In numerical analysis we study the analysis behind the methods. Ie. their rate of convergence, the conditions for convergence, the theory behind optimally using these methods ext.

Numerical methods are just the methods we use to solve stuff numerically.

For instance the Newton Raphson method is a method for finding roots to equations. In numerical analysis we study when this method will find a root, show that it has quadratic convergence, find conditions for when we would expect the converge to be even faster than quadratic ext.

[indigo fox]

at the bottom of the image where it mentions theorem 18, I don't quite understand the wording. What does it mean for it to interpolate f(x) at the gaussian quadrature points?

[heady mesa]

It means the interpolating polynomial = the function you are interpolating at those points

[wary oriole]

I asked this question but would like to hear your opinion. Why always C instead of C++ for scientific computing? Appearently not performance. What if you want to develop large scale library? You would give up OOP for simplicity? I would stick to C++ for benefit over the long run. But yet all scientific community pick C for some reason

[cursive yew]

Well honestly I'm somewhat new to the field. Just to give you context: completed my bachelor's level of physics major a couple of days ago. Not a computer science guy either, only know some basic stuff. My professor advised me C or FORTRAN over c++ and Julia when I was in my 2nd year. Back then I chose Fortran. Recently i developed a taste for C. For some reason it feels a lot more intuitive to me than OOP. Maybe it is compatible with the way I think of algorithms in my head. Not entirely sure.

Now coming back to your question on why the scientific community prefers C. I'm not sure of this but back when I saw this online archive of open source scientific computing codes, I found that c, c++, Fortran were very common. So not just c, c++ is used a lot as well

One of my friends, who works in fluids, had to use some fluid dynamics code (I forgot the name, something along the lines of "basilisk" i think). As far as I remember, she had to learn c++ instead, since the code is in c++. My supervisor advised C for my many body simulations project

[wary oriole]

Most classes i saw use C or Julia. or python and matlab but that's scripting and prototyping side. I fancy late develop large project so want to stick to C++ and OOP

[cursive yew]

Yeah i guess. I cannot talk about the courses but at least my reason for C is just the general interest

[wary oriole]

oh some people mentioned C has much better API compatibility with other languages, so a library written in C is much easier to use in projects written in other languages.

[cursive yew]

Oh yeah that. The starting of this video: https://youtu.be/443UNeGrFoM?si=x7cbQdsN98f3hSvw shares why the person uploading the video prefers C

I somewhat agree with him

[cursive yew]

(though his ideology comes from experience and mine from the interest I have and less on experience)

[wary oriole]

i think it's a just math thing. in other fields like game programming and likes C++ is a lot more

[cursive yew]

Yes I think for games, oop is a better option, though functional programming can work as well (not sure how better).
Roller coaster tycoon is written in assembly language

[dense echo]

The API compatibility thing is real, I'd say. If you take a random programming language, it will fall into one of two classes:

1. It has reasonably good facilities for calling a C library.
2. It doesn't let you easily call code written in any other language than itself.
   C is generally the lowest common denominator that all language implementors make interoperability for first.

Doesn't mean it is any easier to write a library in the first place (C has a wide selection of footguns to pick from), but once you manage to write it, you will have maximized the possible audience for it.

[tropic garnet]

just as a meet-and-greet, who around here is the person to ask for numerical verification in PDEs? as in, using numerical methods to establish rigorous theoretical results about existence, regularity, blowups etc. in PDEs?

some also call it computer-assisted proofs

[plucky kayak]

watch this if you want to learn how differently C++ compiler handles function names compared to C https://youtu.be/NeOTr0u7ALk?si=AZeosQJFpQ3YMsk9&t=263

[fathom rain]

proof-assitants don't use any numerics

[tropic garnet]

I am not referring to Lean or Coq
but rather, for instance, Galerkin numerical calculations with interval arithmetic for an approximation (with bounded errors) and then using fixed point theorems to obtain true solutions etc.

[tropic garnet]

people are really jumping on this trend of using the computer to discover and study specific PDE scenarios rigorously

[fathom rain]

who is "people" and don't choose subject because of trends... 😛

[tropic garnet]

well I don't want to name people directly but arxiv should tell you

[fathom rain]

the only person i know is warwick tucker (took a course for him on interval arithmetics and validated numerics many years ago). never seen much about that on any conferences, but probably just a different crowd

[tropic garnet]

interesting. australia huh.
I'm thinking of some people from East Coast and Canada.

[fathom rain]

he was at cornell and uppsala before

[tropic garnet]

Mitsuhiro T. Nakao wrote a book about the subject in 2019 (Numerical verification)

[fathom rain]

good luck, seems nice

[indigo fox]

For this question, by definition the weights of the gauss quadrature is the sum of squares of the orthonormal polys at each point, however the solutions suggest that you can use lagrange's polys to find the weight. What does lagrange's poly have to do with the weights?

[wary oriole]

oh wow thanks a lot for explaining. Interesting people consider C first for interoperability but ignore C++. i thought they are so closely related

[brave crypt]

OOP is quite overrated. I worked in a $10bn tech startup (e-commerce) before I went to academia, and we actually had to give up OOP for more performance.

I am sure that the scientific computing community also prioritizes performance.

[brave crypt]

Iirc there is a group at waseda university that has been doing this for a long time

[autumn sphinx]

Maintainability is becoming more of a concern in science community

At least in my corner of structural biology

[dense echo]

"OOP or not" is a false dichotomy anyway. One often gets trade-offs between performance and maintainability, but nobody says you have to resolve them the same way everywhere in your project. Figure out where the performance-critical parts of the code are and make a deliberate decision to sacrifice readability there; that doesn't mean you can't still strive for a more straightforward style in the surrounding code.

[random plume]

I mean, to me there's a freaking ton of "one-commit" code. Put out for the paper in a github that will never be touched ever again but never iterated upon, even if it sometimes isn't fully working or contains egregious performance pitfalls (like calculating the n-1 smallest eigenvalues of a symmetric matrix when it's equivalent to first just the largest one)

[tropic garnet]

oh? do you have names or references

[rain furnace]

can anyone help with applying Gauss Legendre numerical integration to the elemental mass, advection and diffusion matrices?

[brave crypt]

I will just refer you to this paper
https://www.sciencedirect.com/science/article/abs/pii/S1007570421004883?via%3Dihub

Nakao is one of them, which you have mentioned before.

[worthy geyser]

When using 1d newton raphson, we terminate once you either reach N iterations or when something enters a threshold range. What is that something and what would it be for 2d?

[plucky kayak]

there are many criteria to choose from

the naive is $||x_{k+1}-x_k|| \leq \varepsilon$ for any dimension

[pine jetty]

Transparent Elemental

[wicked pecan]

And what exactly is |||| doing

[plucky kayak]

for 1 dimensional variant you can use the geometric sequence-based approach $$\left|\frac{x_{k+1}-x_k}{1 - \frac{x_{k+1}-x_k}{x_k - x_{k-1}}}\right| \leq \varepsilon$$

[pine jetty]

Transparent Elemental

[plucky kayak]

norm

[worthy geyser]

thanks 👍

[heady mesa]

Usually the infinity norm, the converge rate proofs for the Newton Raphson method in Rn uses the infinity norm (ie. the maximum error of the components)

[wicked pecan]

Are most norms whose names also exist for metrics just equal to the distance to 0 (or the 0 vector)?

[plucky kayak]

every normed space is a metric space, so yes

[indigo fox]

bump

[grave spoke]

What's a sufficient and necessary condition for an n-point quadrature rule to integrate polynomial functions of maximum degree n - 1 exactly?

[raven sequoia]

I was trying to find the intersection of these curves

this is what wolfram showed me, are they real ?

it should be real as the graph clearly touches, but not sure why it shows up like this

[plucky kayak]

cardano's formula gives real roots even if you have to work with complex numbers along the way

[raven sequoia]

i mean, whether the solutions given real or complex?

I feel like all the 3 solutions should be real

[plucky kayak]

they are real

as you can see this one simplifies to 0.854 * e^(i * pi)

[raven sequoia]

I was trying to solve this in sage math

It's showing it as a complex number, but i notice the complex part is tiny tiny

[plucky kayak]

yes that's just an artifact of floating point operations

[raven sequoia]

any ideas how to get a decimal approximation for this?

using sage?

i can take the real part i suppose, but do you know any other ways to get it done

[plucky kayak]

well maybe the solve function has some parameter to specify real solution, I've never used sage

[fathom rain]

julia> using Roots
julia> @. f(x)=x^3 - x^2 - x + 1/2
julia> fzeros(f,-1,2)
3-element Vector{Float64}:
 -0.8546376797184614
  0.4030317167626848
  1.4516059629557767

[wary oriole]

hey guys im wondering why all these linear algebra library (scipy for example) doesn't have iterative solver for linear system? like jacobi, gauss-seidel, SOR, or efficient ones like Conjugate Gradient, multigrid?

it seems well studied, parallelized, easy to implement. wondering why that's the case

it has advantages compare to direct solvers. so i think it's legit to include in a library?

[plucky kayak]

at least according to scipy docs, they call LAPACK and as far as I can see LAPACK doesn't provide iterative solvers

[fathom rain]

https://iterativesolvers.julialinearalgebra.org/stable/

[plucky kayak]

also https://eigen.tuxfamily.org/dox/group__IterativeLinearSolvers__Module.html

[wary oriole]

interesting, so gold standard LAPACK doesn't have any iterative solvers