numerical-analysis

[wary oriole]

Heath

[final egret]

basic analysis

then read a lot of papers on implementations of numerical methods today

[lime moat]

Kyle Novak's Numerical Methods for Scientific Computing is good for intuition and a broad survey of the field, but I find it to not develop the underlying concepts well enough, at least for the NLA portion.

I used William Ford's Numerical Linear Algebra with Applications: Using MATLAB for my class, and Golub & van Loan's Matrix Computations was recommended material that I didn't pick up. I found the Ford to be a good foundation, but certain things (Givens, Householder) weren't explained particularly well IME, and others (CG, GMRES, MINRES) were, at least in the textbook, more a barrage of dense algorithic slop with no intuition or connections than a good exposition. The coverage of everything else (floating-point theorems, LU with partial pivoting, QR, implicit QR, SVD methods, Krylov, Arnoldi, Lanczos, etc.) was fine.

[wide spear]

I like Demmel for NLA

[bitter musk]

just a notational question, why is it that the kth approximation in jacobi method is denoted as x^(k) and not x_(k)

[mystic hazel]

1. continuity: $ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 $ $\$ $\$
2. x momentum : $-\frac{\partial p}{\partial x} + \frac{\partial}{\partial x}\left(2\mu \frac{\partial u}{\partial x}\right) + \frac{\partial}{\partial y}\left(\mu \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right)\right) + \rho(g_x) = 0 $ $\$ $\$
3. y momentum : $-\frac{\partial p}{\partial y} + \frac{\partial}{\partial y}\left(2\mu \frac{\partial v}{\partial y}\right) + \frac{\partial}{\partial x}\left(\mu \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right)\right) + \rho(g_y) = 0 $ $\$ $\$
4. z momentum: $\frac{\partial p}{\partial z} = 0$ $\$ $\$
5. where $ \mu = \frac{1}{2}(A^{-n})(\dot{\epsilon}^{n-1})$ $\$ $\$
6. where $\dot{\epsilon} = \sqrt{\frac{1}{2}\left(\frac{\partial u}{\partial x}\right)^2 + \frac{1}{2}\left(\frac{\partial v}{\partial y}\right)^2 + \frac{1}{4}\left(\left(\frac{\partial u}{\partial y}\right)^2 + \left(\frac{\partial v}{\partial x}\right)^2 \right)}$

[pine jetty]

NOTHING

[mystic hazel]

well, I am working on fluid mechanics of glaciers, and then I went for 2D problem of NS equations,
Assumptions:

1. steady flow
2. Reynolds number is very small as of the magnitude of 10^-6, because of which the inertia term is very small and hence 0
3. no velocity and no gradient in z direction....
4. here epsilon dot is strain rate, and rho is density, u and v, x and y component of velocity, , mu is viscosity, p is pressure
5. A, n , rho, g_x and g_y are constants

now what do you think will be the best numerical method to solve this coupled PDE

[wide spear]

Superscripts for iterations is pretty common in various areas of numerical analysis

[wide spear]

Have you already tried any methods

[wide spear]

I've just had the realization that RK4 is quite similar to Simpson's rule

[crisp lake]

Why are two-point Padé affairs rigidly bolted to 0 and infinity?

[smoky rapids]

I'm not sure this fits here, lmk if it doesn't

[crisp lake]

I may have been convinced by Hoang: https://arxiv.org/pdf/1305.6104

[wide spear]

The marginal gain seems quite small

[crisp lake]

Most of it seems to be happening around the interval endpoints, where both CGL & scaled Chebyshev spike hard, but CGL not quite as hard as scaled Chebyshev.

[wide spear]

There is no spike for chebyshev and scaled chebyshev

Only equispaced points have error spiking at the interval endpoints

[crisp lake]

That’s one diagram. There are many.
Figs. 7 & 8:

[wide spear]

Those errors are quite mild

Ok like speaking practically

Both are reasonable sets of points to use

If a method works with one set of points and not the other, I would seriously question the validity of the method because the error behavior between the two is not that different

[crisp lake]

The claim is that despite scaled having a theoretical advantage, CGL is just as good in practice. Comparisons are also made to 9 & 10 point optimal Lebesgue nodes.

The theoretical advantage is also clarified earlier in the paper.

[spice valley]

Hey guys, I am kind of stuck on this task on numerical integration. For subtask i) I tried to just gather whatever I could find in the script and wrote it down. But I am not quite sure if that's what I was supposed to do and how to use it in ii).

[wide spear]

I don't quite understand what the equations you've written mean

Could you explain

[spice valley]

yes! so formula 4.11 is $I_n(f)=\sum_{i=0}^n \lambda_i f\left(x_i\right)$ for approximating $\int_a^b f(x) \mu(x) \mathrm{d} x$. And then through $\int_a^b p(x) \mu(x) \mathrm{d} x=\sum_{i=0}^n f\left(x_i\right) \int_a^b L_{i n}(x) \mu(x) \mathrm{d} x=\sum_{i=0}^n f\left(x_i\right) \lambda_i$ you get $\lambda_i=\int_a^b L_{i n}(x) \mu(x) \mathrm{d} x, \quad i=0, \ldots, n$ for the weights of the quadrature with $L_{i n}(x)$ being the Lagrange interpolationpolynomial.

[pine jetty]

aName

[spice valley]

and I don't really know how integrating on a 2 dimentional space works tbh. so far I've only ever differentiated in multiple dimensions so I just guessed and just used a double sum and both the weights and multiplied them because it I didn't know what else to do?

[wide spear]

You never did multivariable calculus?

[spice valley]

I definitely did for derivatives

but I not for integration

[crisp lake]

You don't need to care about the Lagrange interpolants because they haven't asked you to derive error bounds or anything. As far as I can tell they're just asking you to do the obvious application of a quadrature rule to multiple dimensions repeatedly, independently in each dimension. A double sum and multiplying the weights is basically what they're asking for.

Gauß-Legendre for n=1 is just the midpoint rule.

I don't know what the point of that is. Gauß-Чебышёв-Lobatto a.k.a. Clenshaw-Curtis quadrature all the way.

[spice valley]

haha thank you!!

[solemn cobalt]

@crisp lake U an Arab too?

[crisp lake]

No.

[solemn cobalt]

A Desi or a Slav?

[crisp lake]

Cubana

[tired depot]

I know Desi, but Slav is kinda new word for me

[crisp lake]

Hmm. I could use a closed form for iterated integrals of Chebyshev polynomials.

[crisp lake]

This paper makes it look a little like one might be able to use powers of a sequence of bordered tridiagonal matrices but I'd need to find their eigenvalues.
https://jnaiam.org/wp-content/uploads/2024/02/Modified-Finite-Integration-Method-Using-Chebyshev-Polynomial-for-Solving-Linear-Differential-Equations.pdf

[fathom rain]

what is a "bordered tridiagonal matrix"?

[crisp lake]

Bordering means that some of the edges of the matrix are nonzero.

The first and n-th rows and the first and n-th columns are also potentially nonzero beyond the tridiagonal.

[fathom rain]

ok. what matrix in that paper is of that form?

[crisp lake]

The indefinite integrals of the Chebyshev polynomials.

[fathom rain]

I don't see any tridiagonal matrices in there. But goodluck.

[crisp lake]

It’s not written as a matrix. cf. Lemma 2.2 (iv)

[wide spear]

The trig identities don't help with this?

[crisp lake]

I’m not sure what you’re saying. The recurrence for the indefinite integral of Chebyshev polynomials expresses the indefinite integral of T_n in terms of T_{n+1} and T_{n-1}.

[wide spear]

$T_n(\cos\theta)=\cos(n\theta)$

[pine jetty]

守沢 千秋

[crisp lake]

$\int T_n(x)dx = \frac{1}{2(n+1)}\cdot T_{n+1}(x) - \frac{1}{2(n-1)}T_{n-1}(x)$

[pine jetty]

Nadia/Надя/नादिया/娜迪亚/ نادية

[crisp lake]

The trig stuff isn’t the issue in play.

[wide spear]

I don't know what you mean

Ok you have this form

What form do you want it in

[crisp lake]

That’s one step. Iterating integration an arbitrary number of times is the issue.

[wide spear]

Well you can integrate this can't you

Integrate a few times, see if there is a pattern

[crisp lake]

In principle, major families of orthogonal polynomials should have extensive results about them.

[wide spear]

And if you can't find something online, it's up to you to derive the results that you need

[crisp lake]

Deriving results is not that casually able to be done.

[wide spear]

No but in this case I think it is feasible

[crisp lake]

The basic move is for the sake of integration by parts for Chebyshev expansions of the two-point equal-order Lagrange-Hermite basis polynomials.

[crisp lake]

See the below about integration by parts.

[crisp lake]

I'm not really coming up with anything convenient. Thus far what I have looks like this:

[mossy oyster]

(reposting from #advanced-lounge as I figure this channel is probably a better match)
Any recommendations for a good numerical analysis book? My college uses Numerical Analysis by Richard Burden, is it good or would you folks recommend something else?
Also interested in a book on numerical linear algebra. Looking at either Demmel’s Applied Numerical Linear Algebra or Trefethen and Bau’s NLA, can't decide.

[wide spear]

Both Demmel and T&B are good numerical linear algebra books

Is this a book for a class or for self study?

[mossy oyster]

I'm taking a class on numerical analysis soon. I'll have a textbook (Burden) but I normally like to have extra resources to do additional learning out of. I'm hoping to implement my own (very bare-bones) scientific computing library by the end of the semester, so I'm hoping to find books with a healthy balance of theory and implementation.

[mossy oyster]

Any distinction between the two or is it just a coinflip?

[wide spear]

I've only read Demmel so I can't comment on any differences

[mossy oyster]

Ah fair. I suppose I'll just get the cheaper one lmao.

[grave spoke]

<@&268886789983436800>

another one

snow mod

[crisp lake]

My plot to use integration by parts seems to have foundered.

[harsh wave]

What in the holy chat gpt is this.

[wide spear]

It is not chat gpt

[harsh wave]

It looks like something nasa would write

💔

[hollow lily]

I have a question. I have to perform the GN-method, how do i find the best statrting value? meaning B_0, on which the i iterate for B_1 and so on? I have read the wiki page but i did not find an answer

[wide spear]

Which GN

[cloud tusk]

LaTeX ❌
ChatGPT font ✅

[sonic oyster]

i am havig a hard time calculatig the median. can someone go over it with me

[wide spear]

#probability-statistics

[stable trellis]

does anyone know how people generate the really nice plots of pseudospectra?

if i wanted to write my own code to make pseudospectrum plots is there any advice anyone has

[wide spear]

Is there something you have in mind?

This can be done in matplotlib

[fathom rain]

look at eigtool or Pseudospctra.jl

[stable trellis]

well i think i just want to know how to plot pseudospectra since idk how to specifically generate the contours

[wide spear]

Matplotlib contour plot

[magic spruce]

i liked "elementary numerical analysis" by atkinson when i started off, but it reads more like a methods book tbf. it is still well written though, and there are nice exercises

[tall solar]

I don’t do numerical analysis but I TAd for it before and we had Burton. Seems like an ok book I guess. I’d probably consult it again if I had a reason

[wide spear]

Is this for $u_t+f(u)_x=0$

[pine jetty]

守沢 千秋

[wicked hamlet]

More like $u_t + (k(x)f(u))_x = 0$ where $k$ is a spatially dependent discontinuous fumction

[pine jetty]

Messiah

[cosmic karma]

Hello! Is there an inverse inequality which involves H^{1/2} and L^2 norms for finite elements? If so, can someone point me to a good reference?

[brave crypt]

Hello

Anyone familiar with dirichlet

Neumann or robin

[wide spear]

[brave crypt]

3d diffusion

Initial boundary conditions

I have 2 main cases

Cfl

And dirichlet homogeneous

Neumann and Robin mixed

The code is wrong

[cosmic karma]

I guess I am asking for a reference to compute the discrete H^{1/2} norm as was stated in this forum: https://scicomp.stackexchange.com/questions/42513/how-to-compute-numerically-the-h1-2-norm-of-a-function?noredirect=1&lq=1

[wide spear]

The one in the answer?

Or the one in the question

[cosmic karma]

The one in the answer 🙂

$h^{-1/2}| g |_{0,h}$

[pine jetty]

emphatic_wax

[wide spear]

Ok

[cosmic karma]

I have been trying to find one but I am quite new in the field so I am still not good at this

[wide spear]

So

Well ok

Can you compute the norm of g

[cosmic karma]

Yes

[wide spear]

Ok and h is a constant right

[cosmic karma]

Yes that is the mesh size...

[wide spear]

Ok so you can compute that

[wide spear]

this

[cosmic karma]

I really just want to see a reference for $|g|{1/2} \leq C h^{-1/2}| g |{0,h}$.

[pine jetty]

emphatic_wax

[cosmic karma]

Yes

Do you know a book or a paper which mentions this inverse inequality?

[wide spear]

I don't think so

[cosmic karma]

okay thank you

[wanton dawn]

I know nothing about mesh stuff, but in PDE there is the trace inequality $|g|{H^{s - 1/2}(\partial D)} \leq C |g|{H^s(D)}$

[pine jetty]

L

[crisp lake]

What is the question? Indefinite integrals of Chebyshev polynomials are a little unusual but not any sort of AI thing.

[wicked hamlet]

Nobody has an idea? 🥲

[wide spear]

When do you need an answer by

I can look in 2 weeks

[uneven fossil]

hi, who do you know about fractal interpolation

[lime moat]

@bitter musk Also, if you want to reference a particular element, this allows the typical subscript notation for matrices ($x^{(k)} = [ x^{(k)}_{ij} ]$)

[pine jetty]

Zach

[lime moat]

Never mind I'm stupid, x^(k) is a vector. Still, same idea, you may want to reference a particular entry x_i^(k) in the vector

[cosmic karma]

Hello! For which values of s does this hold?

[wanton dawn]

$s > 1/2$. See Taylor's PDE book volume 1 for a proof.

[pine jetty]

L

[steady lantern]

Hi, if you have to solve a PDE which has all boundaries as neumann, how do you solve such a PDE?

I know a solution to all neumann boundaries do not exist so something has to be done about it.

[uneven fossil]

I'm not sure what exactly you mean, can you give me example?

[steady lantern]

This is an example

https://math.stackexchange.com/questions/1690842/numerically-solving-a-poisson-equation-with-neumann-boundary-conditions
also attaching the problem for reference

[uneven fossil]

ok, meet DM

[wide spear]

Are you doing finite differences

[steady lantern]

Yes

[wide spear]

Wait I don’t understand where your confusion is

Is the MSE post unclear somewhere?

[uneven fossil]

I think it would be good to set up a hierarchy diagram.

[wide spear]

(Also I’m about to board a flight so responses may be delayed)

[steady lantern]

uh yeah so the problem I am trying to numerically simulate is a bit complicated, it is transient and it has to go to eventually go to a state X. Now the issue with MSE post is this:

1. I don't know how to solve this in a transient manner.
2. Making sigma = 0 in my case doesn't make sense as my state X (eventually what it should reach) is not zero, and is some finite quantity I know beforehand.
   I don't how to solve this numerically

[wide spear]

Transient = time derivative?

[steady lantern]

yes

[wide spear]

So you don’t have a poisson equation

You have something else

[steady lantern]

It's a bit complicated, not a pure poisson, but something similar.

you can consider something similar to navier stokes

[wide spear]

Ok so you have an initial condition

[fathom rain]

just impose the neumann bc with 1st order fdm, standard way (set top left and bottom right element to -1)

(and keep the other ones 1)

fully diagonalizable

[steady lantern]

set top left and bottom right element to -1
This was one of the cheap ways of doing it, however it pollutes the solution around the neighbourhood (not that it affects my solution, but I feel it's a really weird way of converging)

[fathom rain]

?

[uneven fossil]

You can change the second term to a second-order difference formula and divide the given section to write a positive or negative difference formula.

[sonic wolf]

anyone have advice for b?

[reef patio]

what did you do alr at b?

Did you show the abs(en) smaller than ...

[reef patio]

?

[crisp lake]

The two-point equal-order Traub-Lagrange-Hermite basis looks tremendously better in the Bernstein basis than the Chebyshev basis:

load("bernstein");
load("functs");
load("orthopoly");
f(n) := block([],
    W:wronskian(makelist(bernstein_poly(k-1,2*n-1,x),k,2*n),x),
    W1:subst(x=0,W),
    W2:subst(x=1,W),
    Rs:makelist(list_matrix_entries(if k <= n then row(W1,k) else row(W2, k-n)),k,1,2*n),
    W3:lreduce(addrow,map(matrix,Rs)),
    return(invert(W3))
);
g(n) := block([],
    W:wronskian(makelist(chebyshev_t(k-1,2*x-1),k,2*n),x),
    W1:subst(x=0,W),
    W2:subst(x=1,W),
    Rs:makelist(list_matrix_entries(if k <= n then row(W1,k) else row(W2, k-n)),k,1,2*n),
    W3:lreduce(addrow,map(matrix,Rs)),
    return(invert(W3))
);
f(4);
g(4);

[crisp lake]

The Bernstein expansion of the two-point equal-order Traub-Lagrange-Hermite basis elements looks like:

(%i7) f(4)
                     [ 1  0  0    0   0   0   0     0   ]
                     [                                  ]
                     [    1                             ]
                     [ 1  -  0    0   0   0   0     0   ]
                     [    7                             ]
                     [                                  ]
                     [    2  1                          ]
                     [ 1  -  --   0   0   0   0     0   ]
                     [    7  42                         ]
                     [                                  ]
                     [    3  1    1                     ]
                     [ 1  -  --  ---  0   0   0     0   ]
                     [    7  14  210                    ]
(%o7)                [                                  ]
                     [                     3  1      1  ]
                     [ 0  0  0    0   1  - -  --  - --- ]
                     [                     7  14    210 ]
                     [                                  ]
                     [                     2  1         ]
                     [ 0  0  0    0   1  - -  --    0   ]
                     [                     7  42        ]
                     [                                  ]
                     [                     1            ]
                     [ 0  0  0    0   1  - -  0     0   ]
                     [                     7            ]
                     [                                  ]
                     [ 0  0  0    0   1   0   0     0   ]

[crisp lake]

The Chebyshev expansion of the two-point equal-order Traub-Lagrange-Hermite basis elements looks like:

(%i8) g(4)
      [   1       81      17       5       1        81     17        5   ]
      [   -      ----    ----    -----     -     - ----   ----   - ----- ]
      [   2      1024    2048    12288     2       1024   2048     12288 ]
      [                                                                  ]
      [   1225    201      13       5     1225     201     13        5   ]
      [ - ----  - ----  - ----  - -----   ----   - ----   ----   - ----- ]
      [   2048    4096    4096    49152   2048     4096   4096     49152 ]
      [                                                                  ]
      [           175      47      5              175       47     5     ]
      [   0     - ----  - ----  - ----     0      ----   - ----   ----   ]
      [           2048    4096    8192            2048     4096   8192   ]
      [                                                                  ]
      [  245     245      21       3       245    245       21      3    ]
      [  ----    ----    ----    -----   - ----   ----   - ----   -----  ]
      [  2048    4096    4096    16384     2048   4096     4096   16384  ]
(%o8) [                                                                  ]
      [           7       7       1                 7      7        1    ]
      [   0      ----    ----    ----      0     - ----   ----   - ----  ]
      [          1024    2048    4096              1024   2048     4096  ]
      [                                                                  ]
      [    49      49      9        5      49       49     9         5   ]
      [ - ----  - ----  - ----  - -----   ----   - ----   ----   - ----- ]
      [   2048    4096    4096    49152   2048     4096   4096     49152 ]
      [                                                                  ]
      [            1       1        1              1        1       1    ]
      [   0     - ----  - ----  - -----    0      ----   - ----   -----  ]
      [           2048    4096    24576           2048     4096   24576  ]
      [                                                                  ]
      [   5       5       1        1        5      5        1       1    ]
      [  ----    ----    ----    -----   - ----   ----   - ----   -----  ]
      [  2048    4096    4096    49152     2048   4096     4096   49152  ]

A 4-block block-diagonal matrix with triangular blocks is a great deal more tractable than just a completely dense matrix, though maybe there's a pattern to it I've not recognised.

[crisp lake]

Proving that matrix structure looks more relevant than bringing Chebyshev polynomials into the mix.

[tall solar]

What are you interested in? Idk if I just don’t know enough NA or whether there is more context in this thread.

[crisp lake]

It's pursuant to efficient implementation of a scheme for C^d piecewise rational approximants from a 1975 paper by Jan Gelfgren.

[crisp lake]

I think it’s this one.

[crisp lake]

Sorry it’s this: they all got named the same

[stuck rock]

hey so has anyone here heard of polynomial approximation method by using medians? im really struggling trying to find any relevant information about it

[stuck rock]

guess not

[lime moat]

Taylor I mentioned

[vivid hemlock]

I'm not sure if this is the right channel to ask; but why is it advise to use windowing with STFT?

[brave crypt]

what does absolute tolerance mean here?

[vagrant jackal]

as in, the maximum absolute error (actual - approximate) rather than relative error (absolute error / value)

[brave crypt]

thankyou

[mighty turret]

Can anyone please recommend any material, textbooks, YouTube lectures or universities to learn numerical analysis from. I may have to drop the class and come back next semester. However I truly do want to understand how to use numerical analysis. Right now I’m using teatime numerical analysis as my textbook. Unfortunately it doesn’t derive or explain when to use one method over another, only here’s a bunch of tools. My professor also seems to have just used AI to generate a bunch of questions for the homework, and I can’t even answer 10 out of the 152 part assignment for the first homework. Yes I counted. There are 4 more assignments like this in addition to coding assignments in Matlab which essentially take a week to finish. I’ve spent every free moment I’ve had for the last 4 weeks reading and trying to understand what to do, but I learn the best by doing and seeing be done. Both things that unfortunately are in sort supply.

I will give anything a try, I truly just want to understand what I’m doing, and how to use numerical analysis.

Thank you for your time and have a great day.

[wide spear]

Iserles has a book

What topics does your class cover

[mighty turret]

I don’t know were my lapse in understanding is but we are right now we just got done with root finding methods. Newton method, fixed point and secant. I get bisection though. My biggest hangup is choosing one method over the other.

I unfortunately struggle to identify big Oh, the pace of convergence, and apparently newtons difference is not meshing.

What was the name of that book by iserles by the way

[wide spear]

What is a problem that you have struggled with

[mighty turret]

I have no material that I can use to test if I actually understand what I am doing. The question I have from my professor include knowing the difference between big-oh for computer science and numerical analysis. It’s not crazy to figure out, just not something we talked about. This is basically the same pattern for all questions he asks us to do. I think I have the formulas I just don’t have an understanding of where I stand because I either don’t know how to answer a particular question. Or at the very least do not understand how to visualize what I need to do.

I understand it’s difficult to help me solve this issue when I don’t know if or where I am lost.

Also that book you recommended, is it only on differential equations or is there a more general book that I am missing.

[wide spear]

Yes it's focused on differential equations

I was unsure of what sorts of numerical analysis you are covering

[mighty turret]

I forget can we post pictures or no

Well never mind, I wrote down the 25 topics that we have covered so far.

Rate of convergence (equation), behavior of convergence (equations), error types, machine arithmetic, Taylor polynomial approximations (error term and boundary), newtons divided difference, root finding (biscection, fixed point, newtons method, and secant method). And we are about to cover ronga-kuta here soon, I think

[prisma sandal]

anyone did french canadian studies for numerical analysis?

[naive sonnet]

hey guys i need a little help with the following problem, is there some clever way to solve it? I can solve it but it just require a whole lot of algebra and i wonder if it's what they expect from us.

Find $c_{1,2,3}\in\mathbb{R}$ s.t. $$ \intop_{a}^{b}f\left(x\right)dx\approx h\left(c_{1}f\left(a\right)+c_{2}f\left(a+\frac{h}{4}\right)+c_{3}f\left(a+\frac{h}{2}\right)+c_{2}f\left(a+\frac{3h}{4}\right)+c_{1}f\left(b\right)\right)$$ for $h=b-a$ will be exact for degree 5 polynomials.

[pine jetty]

Henry_quite_hungry

[wide spear]

You can make this easier by setting a=-1 and b=1

But no, no clever ways to solve this

Plug in polynomials and solve the resulting linear system

If you're observant you can reduce the amount of work a bit

[naive sonnet]

wouldn't this reduce this proof to work only with this interval?

[wide spear]

Would it?

Try deriving the coefficients in this interval and then applying them to a different interval

[lime moat]

Kyle Novak, Numerical Methods for Scientific Computing.

Text is good and modern but highly approachable. Bit of a "topics" textbook than crazy depth but good regardless. PDF is free, paper book is about $20

https://www.equalsharepress.com/media/NMFSC.pdf

Mix of NLA, NDEs, and general NA

[median flicker]

I am working on an optimization problem to find the best-fit parameters for my model. I'm using the NSGA-II alogrithm, since have 3 objectives to be satisfied. I have experiment data and a corresponding physics model to predict and simulate the experiments. As objective functions, I'm using the RMSE of the experimental and model values at discrete data points. I'm stuck at a local minimum, as none of the patero front solution sets are fitting the model's predictions. So, can anyone help me understand whether I can take the RSME objective functions or if it's not the best way?

[wide spear]

This is not numerical analysis

Maybe ask in #discrete-math

[quiet wren]

[teal eagle]

Hi great math fellows.
I'm preparing a presentation regarding Richardson's Extrapolation (RE), and I like to prepare one or some slides regarding the geometric aspect of RE if possible (eitehr generic formula of RE, os using an example from differentiation for example). I need some help regarding this. Is it actually possible to give a clear geometric/visual view of RE? Any thought is appreciated.
Thanks.

[tiny lake]

Looks fantastic

[lime moat]

Also I will mention I used William Ford, Numerical Linear Algebra with Applications: Using MATLAB in my numerical linear algebra class FWIW

The van Golub Matrix Computations book is supposed to be very good but I haven't used it personally

[frigid barn]

is numerical analysis related to linear algebra?

[wide spear]

Yes

[crisp lake]

With the help of mechanical codemonkeys, I have results comparing the combined use of Gelfgren (1975), Traub (1964) and Farouki & Rajan (1987) to ordinary cubic splines:

[crisp lake]

It tells you how to implement things like the SVD. The concepts of what the e.g. SVD et al are at a theoretical level probably are better covered elsewhere (Hoffman & Kunze? not sure it goes very deep into various things important in numerics) & the details of coding the SVD et al are basically already done in existing libraries.

[lime moat]

Then yeah, the William Ford book is good, then

[crisp lake]

Oh wow, the vibe coding thing to see what it did without much handholding did indeed go astray.

It's not actually using the Traub Bell polynomial -based formula.

[wide spear]

Are you interested in getting feedback

[crisp lake]

Sure, go for it.

[wide spear]

Are the blue and orange overlapping here

There's a few plots like this

[crisp lake]

I suspect so. I'll have to point that out as a place where it went astray beyond just not actually using the algorithm I told it to.

[wide spear]

Where the blue line appears to be missing

Also chap 2 and chap 3 are both called mathematical background

[crisp lake]

It's all maths. I'm not sure either word really makes sense as a heading.

[wide spear]

Well they appear to be duplicates

[crisp lake]

It's an example of LLM output to critique.

At some point I'll have to complain/remark/etc. that when the particular method is consequential, it's not good at following the directions to produce the end-to-end output in that manner.

[crisp lake]

The AI tells me that what I just cooked up for univariate boundary value problems may be publicable. I have my doubts.

[wide spear]

It is not

[crisp lake]

I haven't shown it to you. But if it's any vindication, the above was meant to say that I don't believe the AI when it says that.

[wide spear]

Ah ok

[crisp lake]

It hadn't even run the benchmarks yet, done any testing, or even put the constraint setup into code yet when it said that. I'll ping once it's got far enough to push to GitHub & regenerate the benchmark report.

[crisp lake]

The exercise is mostly meant to give an idea of what the AI things can & can't do more than to produce anything original.

[crisp lake]

Here we go. It finally finished. There are issues with the degree types chosen that might be giving false overall results but here:
The source is at https://github.com/NadiaYvette/gelfgren/

[wide spear]

Some thoughts:

1. For the discontinuous RHS that is tested first, it would be nice to have some plots of the actual computed solution
2. Also plots of the solution for the highly oscillatory RHS

[crisp lake]

The write-up is written horribly. Once I get through telling it how to do more basics I'll have it rewrite the whole thing. It might be worth having it zoom in on non-pathological cases because a lot won't be visually distinguishable from a distance when everything achieves high accuracy. I think there's still tuning that would help with discontinuous forcing functions.

[crisp lake]

If it really is only worthwhile for smooth problems, nothing has issues dealing with those to begin with.

[wide spear]

Oh it would also be nice to see comparisons to a simple finite difference calculation

[crisp lake]

I'll have it add that, though I think finite difference corresponds to spline model, no?

[crisp lake]

Yeah the polynomial finite differences are splines.

[crisp lake]

I had it do a full re-organisation & tweaked some of the explanations of the implementation I was trying to give. I'll send it spinning on comparing to finite differences & maybe also spectral methods which might take a bit to generate. I might also have it work on more than 1D Poisson, which is just double integration.

[crisp lake]

When I told it to attach it seems to have forced the text I'd already lined up as the message to be sent prematurely. Here is the re-organised write-up.

[crisp lake]

I think it's spun off in a wrong direction & missed the dimensionality vs. degree exchange I wanted to have it do in order to keep the constraints quadratic. Basically

• P = Qu
• P' = Q'u+Qu'
• P'' = Q''u+2Q'u'-QL[u] = Q''u+2Q'u'-QF
  for a 2nd-order linear differential operator L[u] where u=P/Q and the BVP is L[u]=F plus boundary conditions.

[crisp lake]

It sort of doesn't too matter too much if the proposed method etc. actually works for getting a sort of case of development to illustrate the contours of where the AI needed handholding & didn't about the method. Either way, it also kind of did the grinding work to get the thing I wanted to see more about off the ground.

[crisp lake]

Woops that assumes L[u]=u''; maybe the revision should look like:

• L[u] = a(x)u''+b(x)u'+c(x)u
• M[u] = b(x)u'+c(x)u
• L[u]=F(x)
• P = Qu
• P' = Q'u+Qu'
• a(x) P'' = a(x) (Q''u + 2Q'u') + F(x) - M[u]

[crisp lake]

@wide spear Okay, here is a fresh report with comparisons vs. 4th-order finite differences & Chebyshev spectral methods on Helmholtz, advection & Sturm-Liouville problems:
I'm still not sure of everything that's going on b/c cutting the degrees of the constraints to all be quadratic and the degrees of all inequality constraints to be linear should make for a very tractable instance of a KKT problem.

[wide spear]

Lol so I should just stick with finite differences huh

[crisp lake]

I'm not sure how that's a conclusion to draw from the table. Even if you rejected the rational methods in all cases, there would at minimum be a distinction between problem types suited for spectral methods vs. finite differences.

[crisp lake]

I think the extreme advantages it claims in HelmHoltz & ReacDiff would probably wash out in a different form of spectral method from Chebyshev.

There are a lot of ways the thing kind of drives off into a ditch without guidance.

[crisp lake]

@wide spear I updated the special function evaluation things a little. Maybe not completely. It's back into a revision mess but things can get lost in those reorgs so here's a snapshot of that update:

[crisp lake]

Do you get the idea of what I'm doing in all that stuff?

[wide spear]

Roughly yes

Oh it is probably also worth trying a nonlinear problem

[crisp lake]

It keeps it down to quadratic so that when they're translated to KKT constraints, they come out linear. Though there's a bit more to it. There's a moderately complex formula from Traub for equal order Lagrange-Hermite interpolation involving Bell polynomials. The indeterminate always appears differenced with a knot. In the two knot case, the Bell polynomials at the arguments given simplify radically. Then the expressions involving products of powers of the indeterminate differenced with the knots are expressible as Bernstein polynomials of varying rank. Farouki & Rajan then give formulas for rank elevation. With all that in hand, special function approximation then involves interpolation up to some order, additional continuity constraints and consuming the last few degrees of freedom with optimisation or higher derivative evaluation at the endpoints of the grand interval. The boundary value problem case is a little more interesting because it goes something like:

• L[u] = a(x)u''+b(x)u'+c(x)u = F(x)
• M[u] = b(x)u'+c(x)u
• P = Qu
• P' = Qu'+Qu
• P'' = Qu''+2Q'u'+Q''u = (F-M)Q/a+2Q'u'+Q''u
  and then clearing a(x) from the denominator cuts all constraint degrees down to quadratic so the gradient is linear. From there BVPs have a much easier time. And this kind of pattern continues with keeping everything at worst quadratic so optimisations are well-behaved.

[lethal wave]

When solving questions i sometimes see “within the accuracy of 10^-3”, what does that mean exactly?

[wide spear]

Relative error

[silk drum]

Probably $|x_{n}-x_{n-1}|<0.001$

[pine jetty]

mikeliuk

[vapid aurora]

i don't know where to start...

[wide spear]

Do you know what Richardson extrapolation is

[vapid aurora]

barely

[silk drum]

In order to get a good understanding, you probably want to compare two appropriate Taylor series expansions to see how equation (1) was derived and then repeat the process while keeping higher-order terms.

[magic cape]

How are we supposed to know u_(n-1) for the first step of this algorithm?

I presume we're given u_0; then what the heck is u_(-1)?

Will that data be provided too?

[grave spoke]

Do a warm start: Compute u_1 from u_0 using a one-step method of your choice (e.g. forward Euler) before applying the central difference method.

[summer kelp]

What even is numerical analysis about?

I don't understand

I'm seeing like floating point representations, LU-decomposition, QR Householder composition, etc., but what is this entire branch of maths even about?

[wide spear]

It's about numerically solving problems

[summer kelp]

Then how does something like Householder QR decomposition help us with that?

[summer kelp]

And it's not "solving" problems, it's approaching solutions where we want our error to be very small

That's not the same

[wide spear]

QR is useful for linear least squares and eigenvalue problems

Householder is a method for computing QR

[wide spear]

Yes you are correct but if you go around telling people (especially people in other fields) that they aren't solving problems, just approximating the solution, they will roll their eyes at you

[summer kelp]

Sorry I don't want to be mean. I just have a hard time understanding why I'm getting taught this subject

[wide spear]

Ok

Is it a numerical differential equations course or numerical linear algebra

Or more general

[summer kelp]

It's just called "numerical analysis"

The first chapter happened to be about linear algebra

Next one is about "iterative methods"

Like Newton Raphson and stuff

"Iterative methods for equations and systems of equations"

There's a later chapter about differential equations

[wide spear]

In most applied mathematics and related fields like engineering, physics, etc..., exact solutions cannot be found to whatever problems you are interested in

If you need a solution, then numerical solutions are as good as you can get

[silk drum]

From a basic/simple perspective, analysis is about looking at things carefully and numerical is regarding making approximations with rational and real numbers.

So in total, you are looking carefully at different ways of making approximations with rational and real numbers.

People can have differing opinions on which topics are worthy of study in numerical analysis and the prevalence and usefulness of computers makes finite precision arithmetic of interest.

[silk drum]

Historically, it was found in many applications of simulation and modelling in engineering and physics that you ultimately use systems of linear equations to find the approximation to the solution of interest. Direct methods of approximation/solution were used, and later iterative methods of approximation/solution for large and difficult problems became important.

Sadly, it seems your textbooks are doing a poor job of conveying what the subject is and motivating the various topics being taught within the subject matter.

When one looks back on the history of the field of numeral analysis, numerical linear algebra and numerical "solution" of linear systems of equations (and matrices) are very much central and front and center.

[summer kelp]

There is an intro paragraph tho that kinda says why it's important, but I still don't get it tho

[summer kelp]

I just don't like the "it's good enough" mentality when it comes to maths

Like, 10^-24 for example is not 0, tho in numerical analysis it is

[wide spear]

10^-24 is not 0 I agree

However when you are working with people in other fields they do not care about 10^-24

[summer kelp]

(Like physics)

I guess that's just the pure mathematician in me

[wide spear]

When doing anything empirical, you are contending against the messiness of reality

[summer kelp]

I don't think I'm going to like this subject very much :/

[wide spear]

What sort of math do you like

[summer kelp]

Analysis

[grave spoke]

🙄

[summer kelp]

Like

Real and complex analysis

And ig abstract algebra

[wide spear]

The epsilon factors floating around are the same as the epsilons floating around in analysis

[summer kelp]

Yeah I know

But my professor gave this example once

0.000000001 + 0.000001 + 0.00001 + ... + 1000000000000000

Apparently this is really hard to do for computers

So they usually just take the big numbers and say that it's good enough

I don't like that :/

[grave spoke]

Good thing that's not what numerical analysis is about

[summer kelp]

My first chapter was about computers and floating point representations and such

[grave spoke]

Sure, you have limited precision when working with machine numbers

[summer kelp]

My professor also makes references to computers regularly

[wide spear]

Lowkey your professor doesn’t sound very good

[summer kelp]

I think he's fine

[wide spear]

Are computers objectionable

[summer kelp]

Idk what that word means

[wide spear]

Do you dislike the references to computers

[summer kelp]

Not really

They're just whatever

I don't do anything with them

[wide spear]

Numerical analysts do care about floating point error in various ways and there is a lot of work that has gone into mitigating various sources of it

You probably learned that floating point addition is nonassociative

This is a problem

There are various ways of actually implementing floating point addition in hardware, and you can make it associative

[summer kelp]

Yeah

So basically this subject is used more for applied maths?

[wide spear]

Additionally, you can prove that if you are adding a bunch of numbers and you sort them from smallest to largest and then add in this order, you minimize the amount of floating point error among all orderings of the addition

[summer kelp]

Yeah, that's the solution that he gave as well

[wide spear]

Pure mathematicians certainly care about computations as well

https://arxiv.org/abs/1902.05518

[summer kelp]

I'm not saying that I don't care about computations

[wide spear]

I'm going to be honest I have no clue what you're saying

[summer kelp]

Why do we need methods of approximating a solution for a problem, if we can just solve the entire problem precisely?

[wide spear]

Frequently you can't solve the problem exactly

In most situations

[summer kelp]

Okay, but a 'solution' given by a method over 1000 iterations is not a true solution

So why work with that?

[wide spear]

If you need a solution then an approximate solution is better than nothing

You can also quantify how much error working with an approximate solution will introduce

[summer kelp]

I don't like this subject :c

I hope after this semester, I never have to see it again

No hate to you guys, ofc

[cursive heron]

is this chatgpt?

[wide spear]

I think so yeah

[late kettle]

I passed the original text that I initially wrote to an AI grammatical checker so it certainly sounds AI-ish. It didn't add more content; nevertheless, I'm okay with it since the content isn't changed but more organized.

I usually avoid doing that but I had careful reasons at the moment of posting it.

[cursive heron]

I honestly have no idea why someone would want to sound more like a robot and less like a human. You'd think it makes the text more organized, but it actually just hides the original point behind overly flowery adjectives and extremely overused chatgpt cliches

[late kettle]

I agree with you and I think all of my discussions in various channels seemed more clear than this

I care so much about the quality of our math channels here. So, I think I'd rather delete it and as I'm having the time now to discuss all the ideas I'll just prefer to do it on my own pattern. Thank you!

[cursive heron]

It's good practice to write everything yourself, and I think people shouldn't be so worried about having poor grammar or sounding unsophisticated

[wide spear]

[late kettle]

Sure! I've been always siding with that

I also generally like to avoid being sophisticated and put things simply as they could. We're hopefully spending some time here to have a fruitful discussion within the boundaries of our experience and not to sound more intellectual than others. This is good to emphasize on!

I'll be keeping on this manner while I'm returning home and I'd be happy to discuss the same ideas with them more obviously and provide some recommended books or resources that I had came across regarding it further.

[magic cape]

So we're learning about the stability of the central difference method applied to the system [M][\ddot(u)] + [K][u] = [f^(external)]. To do this, we assume f^(external) = 0 and study the homogenous system instead

I want to check if my understanding of why we do this is correct

Because the solution can be written as [u] = [u]_homogenous + [u]_particular, then we know for the algorithm to be stable it must be true that the homogenous part is stable

This is why we can safely set f^(external) = 0 to study the stability

Is this correct?

[wide spear]

Yes

[magic cape]

Awesome, thank you!

[warm lake]

any hint? -_-

i know how to show this :

but i don't know how to show that the norm is less than 1.

[next garden]

I am trying to find a way to approximate highly discontinuous/fractal like functions on R^3. This set of SDFs has some good examples of the kind of complex function I want to try to approximate:
https://tovacinni.github.io/sdf-explorer/

I have recently learnt about Mollifiers, which seem like a nice approach. However, the documents I have on my atm only really look at them in a very abstract way as a means to prove properties of non smooth functions.

What I would like to learn more about is algorithms to define/construct a sequence of mollifiers that progressively approximates the non smooth function more and more.

In R my easiest example would be the weirestrass function. I'd like something like the taylor approximation. i.e. a sequence of functions such that their linear combination becomes a progressively better approximation of the Weirestrass function, as measured by some metric, e.g. the L^2 norm.

A very practical setting I am looking at. If you see the tree ont he set I shared, imagine using an easy to compute "blob" when it is far away from a viewer and adding increasingly more terms as the object gets closer, so that details start popping up more and more.

[wide spear]

Look at e^-(x^2/eps)

For the weierstrass function it’s probably to just truncate the Fourier series though

I’m boarding a flight now but I can discuss alternative approaches later as well

[late kettle]

Wouldn't the discrete wavelet transform work too I think

[next garden]

I am looking for a absis in R^3 though, the Weirestrass was an illustrative example to express what I am looking for

[next garden]

Inw hat sense? What am I doing with the DWT?

i.e. how am I deciding which wavelets to project onto on R^3

[wide spear]

From what I know of SDF work, a common approach is to divide your domain up into small cubes and use some sort of interpolation in each

So polynomial interpolation of various sorts in each box

Or fourier interpolation

[late kettle]

DWT is a way to transform the thing to a different basis for analysis (usually thought about in the sense of Fourier-like Transform).

I think what I was trying to mean is using "wavelet basis". And you can still use ""wavelet basis"" to approximate non-smooth functions using orthogonal-like basis polynomials (and it is often done in various domains for simulations).

The wavelet basis or piecewise linear basis idea is also similar-ish to the interpolation approach. I think they can all be used to define a "signed distance function" as mathematically speaking - you just end up changing the associated "associated metric/distance function" to describe your geometry

Or I might just be too stupid about this approach overall

[next garden]

I agree with the idea, but since function bases are infinite I need a way to pick a subset onto which to project

[wide spear]

Polynomials up to a certain degree, certain number of Fourier modes, etc…

[next garden]

For the wavelets I mean :p. I have some code that does what you say with legendre polynomials. But that kind of code struggles with sharp angles

[wide spear]

Everything will struggle with sharp angles

[pulsar spoke]

Hey, I never got the opportunity to take numerical analysis in my undergrad. Is there a practical book you guys would recommend?

[wary oriole]

Heath is a great one

[keen kraken]

I have PDFs of my math books if you want them.

[primal scarab]

Yeah send.

[craggy nexus]

atkinson we have followed

[mint siren]

i need a book rec for learning to solve pdes numerically, I'm an undergrad

[grave spoke]

what kind of pdes?

[wide spear]

Try Iserles, come back if it doesn’t cover what you need

[mint siren]

any kind. elliptical ones preferably.

[mint siren]

thanks for the rec

[austere glade]

any1 know a good python library for applied numerical analysis with lots of tools for building algorithms

[wide spear]

Numpy and scipy?

[austere glade]

thanks ive only used numpy but I’ll have to try scipy

[wide spear]

Python in practice is not great for numerical analysis because it's slow

There are ways of getting around the slowness by using python libraries which interface with C/C++/Fortran but at some point it's more effective to code in C/C++/Fortran directly

[stable trellis]

do u happen to know any intro references that talk about numerically doing analytic continuation or contour deformation

i'm trying to implement things related to the GW approximation (some comp-chem method) but i'm not sure how those numerical methods work

[wide spear]

You're likely looking at reading papers

Numerical analysis is not a super old field so I don't think the book literature is that well developed

[austere glade]

might try that since it would be advantageous to work in c/c++ anyway for the project 🫡

[wide spear]

[quartz jay]

I thought julia was the current meta for that stuff

[wide spear]

Julia is a complicated topic

Some people really like it

But there are concerns around the stability of the language as a whole

Also if you need to do parallel computing, the lower level you are, the easier it is to do

[tiny lake]

How common is writing stuff in Fortran these days? Is it still used because of its aggressive optimisations? Or is it niche at this point?

[wide spear]

That is also a complicated question

And varies quite a bit by subfield

If you don't need to run stuff on GPU then Fortran is still pretty good

I once rewrote my C++ code in Fortran and got a 4x speedup or something

But

Fortran compiler coverage for various modern hardware architectures is a bit spotty

And the fortran compilers are a lower priority than C/C++

[late kettle]

It is certainly quite common. Some would pretty much prefer the array syntax in Fortran. And also some (including me) would find Fortran a decent nightmare for certain stuff beyond numerics. A lot of CFD codes are written and still being shaped in Fortran. But for my own totality in my personal needs I prefer C/C++

[tiny lake]

I see. I've never written any Fortran, but honestly I'm surprised it's still a thing

I understand that forcing a C/C++ compiler to do some things is difficult at times, but it's still hard for me to imagine

Perhaps I should try to learn it, write some code, and compare the asm output

[late kettle]

Well, Fortran did develop over time didn't it? Also certain people usage boils down to familiarity when the gap isn't huge. It's not bad to try some Fortran too

[tiny lake]

Yep, it's probably something worth experiencing even if only for the sake of seeing it

[wide spear]

It's really not that hard to write

[tiny lake]

Oh, I didn’t mean it’s hard. I don’t know the language, but I imagine it’s quite limited in its functionalities

[wide spear]

It's pretty reasonable

The only thing which is a bit hacky is if you need a C++ vector style resizable array

One really nice thing about Fortran is that there's a compiler directive to check for array bounds when doing array accesses so no tricky seg faults to debug

[austere glade]

i wonder if there’s been any research to find the efficacy of different languages like Fortran for different numerical methods, to determine how significantly better they might perform

it’s probably already a thing actually 😭

[spice kestrel]

If I have a singular 2D integral on a closed contour (say finding u(x,t) with kernel ln|x-x'| or 1/|x-x'| for x' in the contour in the case of Euler equations and SQG equations if that helps), what are some methods to evaluate these numerically? I am thinking assigning nodes around the contour and deal with the contribution of each segment independently, but what should one do for the segment that contains the point of singularity?

[wide spear]

This is complicated because for any given numerical method there are many ways of implementing it

[wide spear]

The key term to search for is kernel regularization

There are many ways of doing it and lots of papers discussing various techniques

Common techniques include skipping the point of the singularity or regularizing the kernel by adding a small epsilon

[spice kestrel]

ty!

[serene finch]

u can do singularity subtraction/product integration on the seltf panel. split
K=S+R where S is the known singular part, integrate S on the self panel analytically & just quadrature the smooth remainder R

[wide spear]

Integrating the singular part analytically is not always so easy

[quartz jay]

Thats funny given that it was supposed to replace common lisp for scientific computing if i recall properly

[scenic wren]

hello.

i wanted to see if anyone would be willing to review my code base for a project I built for extneding my knowledge of python and distribution workloads.

https://github.com/huggablehacker/Collatz-Frontier

let me know if this is allowed or not.

[late kettle]

#computing-software is for programming related questions. but it is heavily not for debugging purposes

[silk drum]

Have you taken into account https://en.wikipedia.org/wiki/Fundamental_theorem_of_software_engineering and Fortran C bindings?

[tiny lake]

Yes. If there’s little code to be written, an inconvenient language might not be a problem, but then people will often be very reluctant to integrate multiple languages into their projects (in general). It there’s a lot of code to be written, then the inconvenience of the language becomes a problem on its own, so unless there’s something to be gained, it might be a bad idea.
I’m talking in general here: I don’t know Fortran. Maybe it’s not that bad of a choice anyway. C is used for its excellent portability. Even if you have a super niche architecture, you can always write your own C compiler for it. Maybe it’s the same with Fortran

[silk drum]

Fortran is good for structured domains and reasonably good for unstructured domains. C is more suitable than Fortran for programming operating systems, and where pointers are useful and required as part of the language design (less friction).

In terms of portability, historically Fortran has been a first-class citizen in HPC so equally as portable as C for applications, and initially probably more portable as a condition of the tender (pointless having a supercomputer that can't run the Fortran applications you bought it for).

[wide spear]

Historically yes but with GPUs everywhere Fortran is not so first class anymore

[spice kestrel]

Ty so much for this. This helps a lot! The simulation is working now!

[kindred drum]

Do yall memorize the error terms for interpolation, root estimation, numerical derivates and intergrals? I just had an exam where I had to recall the error terms for a newton iteration to finding roots, but for the life of me i could not remember. This then lead to me wondering do people really need to memorize these terms, because I dont think it would be hard to just look it every once in a while, or whenever you need it.

[wide spear]

I think most people have a general idea of the error bounds

Maybe not the exact coefficients, but error scaling yes

Trapezoid rule and midpoint rule both have error that scales like O(h^2) and depends on the magnitude of the second derivative

stuff like that

[kindred drum]

Oh I see, think I was too focused on being able to understand the methods and ignored the other major part of why we prefer certain method to others, being the error is less in less iterations and stuff like that

[wide spear]

I would say that understanding the error behavior is part of understanding the method

[wide spear]

I'm not sure that u is the best thing to use for unit roundoff given that u is also your solution variable

Usually it's epsilon

Approach 2 is the better one to use here

Keep in mind that some operations are floating point exact

For questions 3/4/5, one usually ignores the impacts of floating point error on stability/CFL

[summer robin]

yas u its the epsilon machine here for the flaot 16

[wide spear]

I'm saying do not use u because you are already using u in your discretized heat equation

[summer robin]

ah sorry, i ll change it

The question I was asked is to find the condition that ensures the stability of the scheme in Float16

So, I did all the estimates and found this

here G will be the Von Neumann classical condition and \Delta G with errors,

[late kettle]

I think your bound is more of a roundoff growth estimate than a true modified stability criterion

[summer robin]

’Im not even sure what I wrote makes sense. I was asked to apply the von Neumann method to the scheme with errors, examine how those errors affect it stability, and try to find a numerical approximation of the CFL condition to determine

[vagrant jackal]

@dry atlas This channel is for on-topic conversation pertaining to numerical analysis. Please relegate any shitposting in #chill.

[dry atlas]

Ok

Uh

I’ll try to remember that but I kinda forgot what I said here

Pretty sure it was just a question on what that thing was

[granite geode]

I’ve done some numerical analysis before, but I’ll be starting a new course (which covers quite a lot of material, from the basics until solving pdes) from a quite theoretical perspective soon. Do any of you have book or more generally any kind of resource recommendations?

[radiant salmon]

Im trying to compute eigenvalues of $\nabla^2u = -\lambda u$ numerically using Chebyshev Differentiation matrices, however my numerics/code don’t seem to agree with what’s expected $\lambda_n = \frac{\pi^2}{L^2}(m^2-n^2)$ for $m,n\in \mathbb{N}$

[pine jetty]

KySquared

[radiant salmon]

Here’s an image of the numerically computed values & my code:

WAIT I THINK I GET IT
Im making these polynomials on the interval [-1,1] x [-1,1] but on paper im doing [0,1] x [0,1]

[radiant salmon]

Okay, I think I fixed but, but I’m still getting nonsense :(

[wide spear]

Sorry I would try to help but my mind is unable to perceive matlab

[pine jetty]

WildNapkins

[copper ore]

Hey guys I would appreciate some yt channel that explains numerical undergraduate topics with solved examples I am aware that is rare to find that sort of channel but anything helps! thank you

[primal rapids]

+1

[mossy tusk]

Hi guys Does anyone have a beginner-friendly resource for understanding what ‘moments’ mean in the context of Krylov complexity/Krylov methods? Im not sure if this is the best channel to ask the question here but i guess its quite releated to numerical linear algebra

[woven palm]

is this for numerical methods?

[mossy tusk]

idk but it has to do with numerical algebra i guess

[wide spear]

Yes

[mossy tusk]

do you know much about it ?

[wide spear]

https://www.karlin.mff.cuni.cz/~strakos/download/2008_Durham.pdf does this help

[woven palm]

Ty 🙂

[mossy tusk]

I will read it sure but my question is more about something else you know where they define moments which is defined as (O| L^n| O) what i do not get how is that suppose to describe the universal operator growth of a hamiltonian because are not those just the definiton of excpected value raised to some power E(X^n), and secondly i saw in a book that that each term in the series means that the operator growth more or something in the krylov space which i do not get.

where L is the super liouvillian operator

idk if you know much about it

[wide spear]

Oh ok this is not the same

[mossy tusk]

its kind of releated because i looked a bit to what you send. Its sort of releated though

[wide spear]

Oh is it

That's good

[mossy tusk]

yes

[slate jackal]

Hello everyone.
Does anyone know of a resource, whether videos, notes, or a book, that explains difference equations well?

[north wharf]

not exactly

[brazen monolith]

Search up professor Leonard on YouTube

He has a good introductory differential equation course

[grave spoke]

difference equations ≠ differential equations

[brazen monolith]

oh oops misread

sorry nvm then

[half heron]

What are the convergence properties of a discrete-time transition operator when it is constrained to be a symplectic map?

[tropic scaffold]

Hi all. I'm a contributor to a Factorio mod that runs a matrix solver to help users with in-game calculations. Would anyone here know the most stable way of converting a matrix to reduced row echelon form? We recently ran into an edge case with a matrix that had extremely large differences in values (ranging from around 0.1 to 1e9) and a hard-coded tolerance of 1e-10 x (max matrix entry) for zero-checking that wasn't low enough in this case. I could probably lower this further, but I'm wondering if there's a way I could implement the algorithm more safely.

The code is here if anyone's curious, though it's in lua and is kind of clunky: https://github.com/ClaudeMetz/FactoryPlanner/blob/master/modfiles/backend/calculation/matrix_engine.lua#L760

I mostly tried to copy how octave does it, which is a little more concise here: https://github.com/gnu-octave/octave/blob/default/scripts/linear-algebra/rref.m

If there's any good reference algorithms to look at that would be super helpful. Thanks.

[stable trellis]

when solving generalized eigenvalue problem with scipy.eig where the metric is not positive definite, is there a way to get eigenvectors normalized with respect to the metric (i.e. not normalized with respect to euclidean norm)?

[wide spear]

Why do you need RREF?

[tropic scaffold]

Are you asking why we need a matrix solver at all? Or why we're using rref in particular, as opposed to something else?

[wide spear]

Yes why rref as opposed to LU with partial pivoting

[tropic scaffold]

thanks, I'll look into this

[wary oriole]

LU decomposition is the most efficient way of solving linear system. RREF is more solve on paper thing

[tropic scaffold]

I'm not sure how much switching from rref to lu would solve this specific issue, I think it's more of a units scaling thing

[wide spear]

Partial pivoting is what fixes your issue

[tropic scaffold]

I'm already doing that

I can try to describe

[wide spear]

Are you using doubles?

[tropic scaffold]

yes

[wide spear]

I don’t see why there’s an issue then, you can have a range of 10^16 with doubles without floating point issues

[tropic scaffold]

I set a tolerance of 10^-10, I suppose I can make it lower. In this particular case it wasn't low enough.

Also the solver sets up "dummy" variables for things like net inputs or outputs, that are just a column with a single 1 in it and the rest zeros. The units in these could be set to match the max value of its type, which would also solve this issue.

[fluid cliff]

Anyone here who can recommend / list some of the best universities / research groups for numerical analysis? I am researching locations for doing my PhD

[wide spear]

This is kind of a challenging question

Do you have any geographic restrictions

Do you have any particular areas of numerical analysis you want to do

[fluid cliff]

I’m interested in computational fluid dynamics algorithms so relatively practical, I’m however not sure about geographical location. Im European so in Europe would be easiest and in terms of salary it would be a LOT better than whatever is going on in the US. It would also be more line with my masters degrees. Im however debating this topic and therefore looking into which universities have the best offer…

https://www.math.cit.tum.de/en/math/research/groups/numerical-analysis/

Something like Munich does ‘sound’ not bad but I ofcourse dont know anything about experiences by others or possible supervisors

[wide spear]

Fluid dynamics is still a pretty large area

I work in geophysical fluid dynamics and I know a lot of PhD positions are announced on ES_jobs_net which is an email list that you can subscribe to

But I don't know anything on the aerospace, plasma, etc... side

[fathom rain]

do it with somebody connected to dlr, probably around braunschweig

or vki

in the us, stanford, mit, caltech etc

[fathom rain]

if you are interested in staying in academia kth n stockholm is very good

[kind steppe]

https://github.com/0x32-l3git/accelerated-leibniz

How does this create a speedup? I know why it converges to 3-pi but not why it is faster at convergence.

[round aurora]

Guys is it possible to solve this Fredholm equation with Nystrom Method?

[wide spear]

I’d imagine yes

[round aurora]

None could do it lol
I guess imma send it wrong then

[wide spear]

Well what do you have

You might need some modifications because the resulting system is probably nonlinear

[round aurora]

The issue is in MATLAB code
It should compile to a comparing graph between numerical solution and Exact solution

[wide spear]

Good luck

My mind is incapable of perceiving matlab code

[subtle night]

|| ~~What are you supposed to do? An exact or approximate solving of this equation?

I'm asking because it seems you need to find the optimal contraction constant from the integral you're given, which seems to be 3√3 / 16, a rather small number (a tiny bit less than 1/3). Could be relevant if you're currently dealing with a chapter introducing fixed point methods or picard iterations. Idk.~~ ||
I realized this wasn't the answer you were looking for, my bad.

[rich galleon]

|| I don’t know what the nystrom method is, but I’m going to point out the following:

(u(t) + 0.25t - 1)/t = a constant. Sure, a constant equal to the integral involving u, but a constant. Call it K. Then, u(t) = (K-0.25)t + 1.

Thus, K = the integral from 0 to 1 of s/(1+((K-0.25)s+1)^2) = the integral from 1 to K + 0.75 of 1/((K-0.25))^2 * (u-1)/(1+u^2).

Thus, K*(K-0.25)^2 is equal to the integral from 1 to K + 0.75 of (u-1)/(1+u^2). From here, you can get an explicit integral in terms of logarithms and the arctangent function. You would not be able to solve for K analytically, but a numerical approach should get you there pretty quick if you massage the computation right||

||Regardless, u(t) must be a linear function||

[wide spear]

<@&268886789983436800>

[analog basin]

Von karman institute is great
Imperial college, eth and delft are all great too

I think the US is particularly stronger in this area ngl

[rotund yoke]

Hello everyone,
I'm working in research paper on mathematical epidemiology ( fractional calculus).
If you see that you have a great background to contribute, so contact me .

[pearl talon]

is this something that someone would study in numerical analysis, or did i miss it when studying linear algebra?

wish i could learn more about cool applications of linear algebra like this

or would i get more of this if i were to study more cs

[wide spear]

I have no clue what is going on here

[sinful orchid]

wtf am i reading

but this is probably more of a math thing. eg numerical linear algebra

i don’t think the average comp sci student or like software dev really give a shit about matrix factorizations

[wide spear]

Lol I don't think any numerical linear algebra knows what is going on here either

[pearl talon]

SadHampster

[sinful orchid]

TRUE !

i imagine Pr(this twitter post is a shitpost) to be high

[analog basin]

engineers heavily do !

[wide spear]

Average SWE doesn't but engineers definitely do yes

[primal rapids]

hey

im studying numerical and im trying to grasp the concept of it

basically the goal is to get approximated values to optimize the computation time because sometimes u cant get the exact answer or the exact answer take so much computation time

we are focusing on the non linear equations and we have two types of methods:
closed root method : bisection method
false position method
open root method:
newtons method
secants method

how do we know which method is the most efficient to use

and for linear equations , do we not need methods bc they are simpler to solve?

[wide spear]

You don't know which is most efficient

For linear equations you can solve them without any needed approximations

Via solving linear systems such as Gaussian elimination/LU factorization or related techniques

[thin prawn]

basically what they are doing is transforming a matrix A in its closest pure rotation, that is, a matrix in SO(3) that minimizes the distance (this is related to Procrustes orthogonal problem if you want to see), it can be shown that the matrix that satisfies that is the orthonormal matrix present in the polar decomposition A = QS.
in the tweet they use an iterative approach to compute this Q which is Q_0 = A, Q_n = (Q_(n-1) + Q_(n-1)^-T)/2.
you can see that at 3 iterations the rotation is almost identical to the Polar decomposition.
what i dont really know is the quaternion, gramschimdt and cayley thing? i suppose they are extracting the "rotation" components of these factorizations and showcasing it.
I may be wrong but the reason they do this i think is because after applying operations on an orthogonal matrix A that shouldnt alter its orthogonality property it will regardless stop being orthogonal due to accumulated floating point errors, so you have to re-orthonormalize it constantly, and it is faster to use the iterative approach every time you want to re-orthonormalize instead of computing the polar decomposition every time.

and as for your question if you would study this thing in numerical analysis, you may or not, i dont think it is a central part of numerical analysis nor numerical linear algebra, it´s a thing you can just browse if you know the basics of numerical linear algebra, and also the tweet is very poorly-explained so you can´t get much insight from it.

[radiant salmon]

Is there a well to “tell” when a diff. eq. is stiff? Preferably just by looking at it, but anyway to determine it would be nice

I ask since I recently learned that the diffeq im trying to solve numerically is stiff, so finding an appropriate timestep is a pain as I’ve just been doing trial-and-error

[grave spoke]

in practice, when an explicit adaptive numerical integration method performs poorly in the sense that it tries to shrink the time step over and over again

[wide spear]

You should change your time stepping method instead of just shrinking the time step

[wary oriole]

Isn't that well-studied? There are pros and cons for each, from bisection the most robust but slow to newton the fast but not robust. Newton's method has quartic convergence, so if the problem works you should use this. Secant and other's has superlinear convergence.
Most general algo you see in commerical software is hybird. Meaning it will try to run newton as much as possible but if it doesn't work then go back to more robust method, then try fast ones again later.

Safeguarded Methods it's called. like Brent.
If it's system of nonlinear then it's tough. but there are also robust newton-like methods

[sinful orchid]

yeah basically trial and error. your step size has to get smaller and smaller if you want more accuracy. also consider the size of the eigenvalues if it can be expressed as a linear system. large eigenvalues = stiff

[sinful orchid]

my numerical analysis professor mentioned that there are also some hybrid methods where you do bisections for a few iterations then switch to newtons or vice versa.

[kindred drum]

Hello, I'm working on Numerical Linear Algebra by Trefethen & Bau. I don't understand the statement of this problem on chapter 12, am I trying to map data vectors with n coordinates to the interpolating polynomial of p(y)? But the points are equispaced so why are the vectors x called data vectors if they are presumably in the form -1+2i/n i=1,...,n

[fathom rain]

”at”

[summer kelp]

Where can I find comprehensive and reliable information about splines?

[twin sail]

it is not possible to do LU pivoting both for numerical accuracy, and sparsity, yes?

since each of them require their own row ordering?

so short of additional, heuristic, logic (on already greedy pivoting), which would be entirely from scratch, ie there is no magical formula that balances the two, it's best to just pick one or the other?

like, technically if your sparsity algorithm gives you a set of rows to choose from next, you can pick the biggest pivot, but in practice you're not going to get a choice? (ie there will be a single row for best sparsity) so there's no point even considering it?

[wide spear]

If you have a sparse problem then LU is likely not the best method at any rate

[twin sail]

what other methods should I be considering then?

[wide spear]

Iterative solvers tend to be a lot faster for sparse problems

[twin sail]

note that the LU decomposition would be re-used many many times

or, well, hmm

[wide spear]

Maybe look at super LU then

https://portal.nersc.gov/project/sparse/superlu/

What sort of problem is your matrix coming from?

[twin sail]

circuit analysis

[wide spear]

Oh hmmmm

Is it symmetric?

[twin sail]

nope

[wide spear]

Positive definite?

Ok yeah use super LU or gmres then

Super LU is available in scipy

[twin sail]

hmm I'm trying to write all the algorithms myself though

it's a learning project, not for some super optimised purpose

and I'm not using python or C/C++ regardless

[wide spear]

https://dl.acm.org/doi/abs/10.1145/1089014.1089017

You can implement from this description then

GMRES is not too hard to implement either

[twin sail]

hmm yes I see that

in all my courses there has been a huge emphasis on LU, though we did cover iterative

particularly in circuit analysis, I've only heard reference to LU

I don't suppose you know of a special case of iterative methods being worse for circuit problems?

[wide spear]

Well it will work without too much thought

Easy to use is what ends up determining what gets used in practice

[twin sail]

iterative methods aren't guaranteed to converge, yes? or is GMRES different in some way I didn't see?

[wide spear]

They are

[twin sail]

since I assume it's a problem of generality, where circuit systems can end up in states that are extremely nonlinear and thus the corresponding jacobians have high/worse spectral radius?

[wide spear]

The convergence of iterative methods is a bit of a delicate matter and the more ill conditioned a matrix is, the slower an iterative method will converge

But ill conditioned matrices are not good for non-iterative solvers either

[twin sail]

ill conditioned matrices produce greater error in deterministic algorithms, and slower in convergence in iterative methods, yes?

[wide spear]

Yes

[twin sail]

also, gmres requires solving least squares? Is sparse qr/(householder?) many times really still more efficient than sparse LU?

[twin sail]

but the error will still be present in iterative methods, right? Or do iterative methods not suffer as directly from error from ill conditioning?

[wide spear]

GMRES doesn't require solving least squares

Or

The least squares problem is much easier

[twin sail]

oh yes I see, the qr decomposition is effectively an outer product update or whatever it's equivalent to, that kind of minimal update

[wide spear]

Because you don't have to do a new least square solve each time

You update the QR factorization at each step

[twin sail]

and gmres is guaranteed to converge, with infinite precision arithmetic? as long as A is invertible?

hmm ok I'll consider iterative then

regardless, I did want to know for LU, if sparsity pivoting intrinsically makes it nontrivial to pivot for big numbers/numerical stability?

[wide spear]

GMRES converges even if A isn't invertible

Min res

[twin sail]

oh whoops yes I see

[summer robin]

Hello, could you please tell me if it is correct to write this, knowing that i took all the rounding errors equals to \delta same everywhere? (the professor who told me to take the same delta)

I feel like it’s not mathematically rigorous, but I can’t quite explain why!!

[summer kelp]

Anyone have an easy proof for this theorem:
Let $B \in \mathbb{N}{\geq 2}$ and $x \in \mathbb{R} \backslash {0}$. Then $x$ can always be written as
[
x = \sigma B^N \sum{i=1}^{\infty} d_{-i}B^{-i} = \sigma B^N(0.d_{-1}d_{-2}d_{-3} \cdots)B
]
With $\sigma \in {-1,1}, N \in \mathbb{Z}$ and $d{-i} \in {0,1,2,\dots ,B-1}$

This representation is unique if the following conditions are met:
\begin{itemize}
\item $d_{-1} \neq 0$
\item $\forall n \in \mathbb{N}, \exists i \geq n$ s.t; $d_{-i} \neq B-1$
\end{itemize}

[pine jetty]

Nico

[summer kelp]

?

I have a proof for this in my syllabus, but it kinda sucks

[cunning raven]

of which part

existence is reasonably easy

uniqueness is harder, but can more or less be proved via algorithm

[summer kelp]

Uniqueness, yes

[cunning raven]

clearly we may restrict to positives. suppose we have x!=0. then it is easy to see that for exactly one integer N, x*B^(-N) is in [1/B,1), which is implied by d_{-1}!=0 (1 is excluded by the second rule).
thus we are left with showing uniqueness in this interval. suppose that 0.a1a2a3... and 0.b1b2b3... are equal as numbers. for the sake of contradiction, suppose they are different as sequences. then there is some minimal n such that bn!=an. assume without loss of generality that bn>an. we now reason that 0.b1b2b3...>0.a1a2a3..., because 0.b1b2...bn>0.a1a2...an..., because the sum of all digits strictly after an is at most B^-n, and it cannot be equal to B^-n because of the second rule.

as such we have a contradiction and uniqueness is proved

the important detail here is that we need to prove that the only expansion of 1 of the form 0.a1a2a3... is with B-1=a1=a2=a3=...

[summer kelp]

Thank you!

[still forum]

https://arxiv.org/abs/1903.12642

made my day

[wide spear]

😵‍💫

[still forum]

one day we will have spectral element methods using sphered tetrahedrons

[wide spear]

We do already

[quartz pumice]

I’m gay

Numerical analysis goes here?

What about optimization?

[tardy ember]

Yeah that stuff pretty much all goes here

[next pecan]

Oh yessssss

Whenever anyone discovers an interesting preconditioning method, please don't hesitate to ping me

[plucky junco]

How can I search for the most optimum way to cut this green rectangle to form the arch between the two red arched lines.

[prisma sluice]

3d printing?

[slender star]

What do those curves need to satisfy

[brave crypt]

Does anyone know how to go about bounding the global truncation error in a solution, to an initial value problem, produced by the RK4 algorithm? For my particular case I would be satisfied if I could show that the error is less than 10^(-4). For anyone interested, you can find a more detailed description of my problem here https://math.stackexchange.com/questions/3228133/bounding-the-error-in-a-solution-to-an-ivp-produced-by-rk4

[quartz pumice]

this algorithm has time complexity T(n) = c1 + n*c2 + n*n*c3 + 0.5*n*n*(n-1)*c4 right?

where c1 = time taken for things outside i-loop,
c2 = time taken for things outside j-loop
c3= time taken for things outside k-loop
c4 = time taken for things inside k-loop

A,B are arrays of n integers

Hence it is O(n^3) right?

[glass solstice]

Ye, assuming c_whatever don't mess up that complexity

[quartz pumice]

wdym

the c_i are just constants representing the time taken for the primitive ops, which is irrelevant to the big O

[glass solstice]

Ah, i assumed the c were constant, but just in case they werent

ya know

constants don't mess it up

c usually are constant, I know, but I've seen some very weird notations before

[summer egret]

aight, I’m stuck at an exercise. this is gonna take a moment to write everything down though

So, we’re looking at solvers for hamiltonian systems, and there’s this exercise:

the symplectic euler is this one

so I started to compute the truncation error, componentwise (unsure if that’s what I have to do, actually), and this is how far I got (for the first component):

\begin{align}
T(\Delta t) &= \frac{1}{\Delta t} \left( p(t_{k+1}) - p^{k+1} \right) \
&= \frac{1}{\Delta t} \left ( p(t_k) + \Delta t \dot p(t_k) + \mathcal{O}(\Delta t^2) - p(t_k) + \Delta t \frac{\partial H}{\partial q}(p^{k+1}, q^k)\right) \
&= - \frac{\partial H}{\partial q}(p^k, q^k) + \frac{\partial H}{\partial q} (p^{k+1}, q^k) + \mathcal{O}(\Delta t)
\end{align}

[pine jetty]

Sascha Baer:

[summer egret]

and I can’t simplify this any further, or at least I don’t see how

any ideas?

[summer egret]

access denied, I think the image is from a private discord server or sth

(not that I could help, the only search trees I'm familiar with are ab-trees)

[fathom latch]

nvm thank you anyways

i found an online visualizer

[jovial flame]

what's the fastest algorithm for computing a legendre symbol? use euler's criterion?

assume its modulo a large prime p

more specifically, for x in F_p for p a large prime, i want to determine if x has a square root in F_p, but don't actually need to compute the square root

[median plinth]

I think using the techniques described in the computational examples section of the wiki page on legendre symbols is the best way to go

Think there's no one best way and the best algorithm would differ depending on how large your number was, or if your number was prime or not, x that is

[elfin arch]

well modular exponentiation is pretty fast, no?

so the euler criterion isn't going to be that bad unless the prime is truly huge

[golden mango]

Depending on the implementation, yeah

For example, in Python you wouldn't want to use the default exponentiation operator, you'd use the pow() function

[median plinth]

I wonder if using Gauss's lemma would be any faster

There are optimization to be made rather than calculating each individual term in the sequence

[golden mango]

I don't know Gauss's lemma so I can't help you there

[viscid crag]

anyone good with minimax?

algorithm?

[rough perch]

@viscid crag sure

[distant pewter]

hi guys is there anyone who would help me trying to figure out how to calculate random numbers for a slot machine software that we are developing?

[slender star]

Depends on what kind of random numbers you're looking for

[distant pewter]

basically i have three major groups that i have to connect to each other in order to choose a combo

the problem is that some numbers that i could have every spin are way too big

im thinking in using hexadecimal numbers

[tame trellis]

what language are you using?

if you are using some language with native bigint like Python, I don't see an issue

if you are using a language which does not, just use a library for bigint

[distant pewter]

mh for the reduction of the random number that we will generate in order to match the combo, it will be good to use the modulus operator?

[tame trellis]

it depnds, sometimes doing a modulo introduces a bias towards certain numbers

you would want to make sure you are careful about that

[distant pewter]

correct we actually need a bias custom that we can regulate

i think

maybe random

[tame trellis]

what exactly are you trying to do

can you describe the original problem instead

[wet river]

advanced math

[distant pewter]

Ok basically I'm creating a slot machine, to create the results of the spins we made a file with the permutations of the numbers written in blocks. Now, the combination is retrieved from the file by a random value. What we need to accomplish in essence is a good luck rng algorithm. So I wanted to know if it was safe to use a simple mod operator to reduce the number generated and if it was correct to use this operator for this algorithm, for example I have a range of numbers: 0.50, but the random number is 26543245

[tame trellis]

what are you using for random numbers

why not just scale the range

your rng's range -> [0, 100), your required range -> [50, 100) => divide rng/2 and add to 50

[distant pewter]

we use cluster filled with entrophy

ye but if i have 3 milion

and it has to be inside 0,50

it cant work

[tame trellis]

divide by 6 million?

I don't understand what you mean

[distant pewter]

ok sorry im dumb as fuck, makes sense, i have to be sure that every number will be inside the range tho

[tame trellis]

yes, divide and handle the edges accordingly, if you have trouble, paste your code, I'll take a look

[distant pewter]

because it dipends on the linear scale otherwise every reduced number will be 50

well handling the edges is the issue

i should divide 3 milion by every zero till 50

[quartz pumice]

https://turingmachinesimulator.com/shared/ibwabwmmsc

@quiet sparrow could u explain how u colored the input even though it is not known beforehand how many symbols it has?

i know how to color input when input string is of length 3, 6, 9, etc. individually (taking 3 states for length 3, 6 states for length 6, ...) but i don't know how to color it in general

[quiet sparrow]

one by one

you bascially go
xxxxyyyyzzzz ->
xxxxxyyyyyzzzzz

if that makes sense

[quartz pumice]

ok got it i recoded the coloring part http://turingmachinesimulator.com/shared/jqyhrteflv

thanks

[fresh oxide]

this seems fun

wut u talking about??

[jovial flame]

Hello! i have two events A_1 and A_2 and each of them has an associated binary value v(A_1) and v(A_2) such that not both of them are 0. is there an easy way to choose an element uniformly at random from the set {A_i : v(A_i) != 0} in constant time? meaning that no matter what the values v(A_i) are the operation to make the decision takes the same amount of time, or is otherwise indistinguishable based on the input?

essentially i don't want any kind of realistic eavesdropper to know which event i've chosen

[fresh oxide]

w8 wut?

i dont even wut the question is

[jovial flame]

what part is confusing

there are 3 possibilities: v(A_1) = 0 and v(A_2) = 1, or v(A_1) = 1 and v(A_2) = 0, or both are 1. in the first case i want the chosen event to be A_2, in the second i want to choose A_1, and the third i want to choose randomly between the 2

but the choice needs to be made in way that doesn't give away information about the v(A_i) if someone can see, for example, the power consumption of the machine, or the time taken to make the decision

[median plinth]

why not just flip a coin and look at A_1 or A_2 with equal probability

and then if that one is 0, then do the other one

[jovial flame]

because that last branch is not constant time

[fresh oxide]

i dont get the part of power consumption or time

[jovial flame]

i'm by no means an expert on constant time coding but i think basic arithmetic can be made constant time

[fresh oxide]

but i guess u could make a seperate clock, decide, and then use the clock to update in timed interval

[jovial flame]

so i think something like this might work
(v(A_1) and not v(A_2)) + 1 - (not v(A_1) and v(A_2)) + (v(A_1) and v(A_2)) * random({0,1})

to select an index

[fresh oxide]

may i ask? why u want to set it for a constant power consumption or time?

[jovial flame]

it's for cryptographic purposes so its needed for security

[fresh oxide]

ah ok

makes sense then

yea so could u do that?

have it time like 1 second

and update after 1 second

but decide before hand?

cuz i dont think even random always operates on uniform time, though i may be wrong

and it depends wut u actually compiling ur code with

cuz some psuedorandom number generators use functions and have variable processing times

not only that, but the way u manipulate variables, even if u add some empty code to the last branch as a fluff, u cant be sure it will require same time cuz cpu also does predictions and caching based on those predictions

[jovial flame]

i guess i meant that the time is independent of the values of the v(A_i). so pretty sure if random({0,1}) has variable time, that would be ok

[fresh oxide]

ok i cant figure out wut u mean, so imma just give up, but best of luck, hopefully u find ur solution 😃

[tame trellis]

I don't see why that is not constant time

don't you get an infinite series that sums to a finite value

[fresh oxide]

@tame trellis wut?

[tame trellis]

rerolling on getting 0-0

leads to constant time overall

[brave crypt]

@jovial flame first, you shouldn't implement your own cryptography for anything other than just learning. Second, flipping the coin should do it.

If you're paranoid about memory access to different places having different timings / power, do a XOR swap with the second argument always multiplied by the value of the coin flip

[dawn viper]

xor swop takes up more time than using temporary variable

Writing to registers is a lot faster

+implementing own crypto not only insecure but also probably slower cuz itll prob not utilise special crypto instructions(aes-ni)
even using preexisting code may be insecure especially rsa(just dont use rsa its too hard to implement properly)

use urandom for secure random numbers

[tame trellis]

not really relevant, but https://github.com/lemire/fastrand I personally like this for random numbers

[golden mango]

That isn't secure though, so you can't use it for cryptography

[steady flint]

use openbsd arc4random for secure random numbers

[dawn viper]

Time tho

its risky to use a black box

[steady flint]

i don't actually know if openbsd's /dev/random code is any better than linux's

[dawn viper]

Tbh jus treat urandom as a random file

[median plinth]

Just ask here it's fine

[dawn viper]

its ok it takes time and some staring at it going like wth to understand some algos

[tulip meadow]

Does anyone know if it's possible to pass in more parameters into a function with only a few input arguments in matlab? I am trying to use ode45 to solve a second order ODE I have, but I want to pass in random values for the constants in my ODE. If anyone here can help, I'd highly appreciate it

For some further clarification, ode45 can accept 3 input agruments: a function, a range, and initial conditions. I want to pass in 5 other constants: m, k, c, F0, w. I tried to create a nested function that takes in the inputs, then the other function would only pass the needed stuff, but it hasn't been working. Thanks.

[blissful crater]

I'm not framiliar with the function, but if I am understanding you right. You want to look at the result for multiple input conditions?

One way to do this is to construct an array with each initial condition you want to look at, then run a for loop with the function you are working with inside of it.

On each iteration of the for loop you can plug in a new initial condition into that function based off of it's index in the array.

[tulip meadow]

I didn't see this before but I managed to get it working by creating a function that takes in all the inputs, then another function with a nested function that only inputs 2 inputs while the constants are defined.

[broken coyote]

How do you program an algorithm like the Riemann Zeta function where by you input a number and the algorithm can tell you what number it approaches ?

[tame trellis]

brute force

most languages natively support complex number operations

[broken coyote]

brute force seems a bit vague since it could take forever to find the number

[tame trellis]

what do you mean?

do you want the exact value?

what do you mean by exact value then?

zeta(3) for example

[broken coyote]

yeh then you want to know what number the function approaches

[tame trellis]

so you do not want exact value?

[broken coyote]

can it ever really get an exact value, you can always add another 1/(n+1)^s term to it

[tame trellis]

the question is too vague, so I'll see myself out

[broken coyote]

plus precision issues of computers on top of that

it really isn't

the function converges for any real number > 1 and you're looking to find what number it converges to - i don't see how that can be vague

[tame trellis]

zeta(3) for example

[broken coyote]

sure

that is a valid example to use for the algorithm i guess

nevermind i'll ask another server

[spring plinth]

Well, you need the right formla

[broken coyote]

well i know that much

the wrong one won't help 😛

[fresh oxide]

but how would i convert any array of numbers, given a maximum value within the array to base 64 string
?

[sage vapor]

convert any array to the empty string

[fresh oxide]

? @sage vapor

[spring plinth]

You can just pretend it is text and convert directly

[fresh oxide]

yea but then it actually makes the shit longer, i wanna shorten it

@spring plinth

[brave crypt]

You haven't even clearly stated what you mean by "longer," what properties you want your "conversion" to have, or what the array numbers is supposed to represent.

[dawn viper]

@broken coyote depends on the function, find a fast converging sequence, like a continued fraction, a sum over exponentially decaying terms

tbh i cant understand how my code for computing zeta works anymore so i cant quite explain lol but i computed eta then zeta

[fresh oxide]

@brave crypt yea nvm posted it on stack overflow and got some amazing answers that im pretty happy with 😃

i already stated array of numbers (with given max value m) to base 64 string, idk wut more properties i could give u

either ways doesnt matter

[earnest ice]

I'm studying this as a degree this September

Anything I should know?

[summer egret]

“this” is a bit broad

you’re starting an applied math degree? at what level?

[earnest ice]

Freshmanevel

The degree is literally called applied and computational maths

[indigo musk]

you should know that degree names aren't really important, better check the curriculum

[summer egret]

yea if you tell us what classes you'll have we might be able to give you more specific tios

honestly, all the advice I can give is pretty ETH-specific, I've kinda gotten a feeling that ETH is way more brutal than most unis

[earnest ice]

First year is Lin algebra, two intro programming modules, functions+numbers, diff eqs, calculus, two physics modules and some other stuff

[acoustic slate]

what are you interested in?

[brave crypt]

hi! how can I extract the forward vector from a transformation matrix?

I have tried obtaining it by multipling it a vector (0, 0, 1, 0)

but it was in vain

[earnest ice]

@acoustic slate mostly ML/AI and quantum computing but that's mainly because that's all I really know about through like pop science books

[supple sandal]

@VerdGehirn does it work if you try vector ( 0, 0, -1, 0 ) ... I'm thinking that if If x is to the right, and y is up, then forward is -z in a right handed coordinate system, and +z if in a left handed coordinate system. Maybe you want -z and you're in a right handed coordinate system?

Maybe you have to normalize the output vector, in case the upper left 3x3 matrix is not orthonormal (e g. Scale != 1 )

Idk what happens if the upper left 3x3 has shear

But if the columns are orthogonal and length = 1 then I think it's a rotation that reorients that vector without scaling it

[vagrant stream]

Can someone link me an example of how to approximate the solutions of a system of equations using scaled column pivoting?

My textbook has none and googling isn't helping me.

[fresh valve]

If I have factorization A=UU' where u is upper triangular, can anyone see how I can rewrite it as A=LL' where L is lower triangular?

[spring plinth]

hmm

U' is transpose of U, so that's lower triangular...

(UU')'=UU', so that means that A is symmetric...

[low sandal]

So this is a little out of the blue

I know Newton's method is looking for fixed points of an iterative method

Does that have anything major to do with Brouwer's fixed point theorem?

[tardy ember]

If anything it'd be closer to Banach fixed point theorem I think

I remember hearing that part of the idea of the inverse function theorem was very similar to Newton's method

[summer egret]

yea iirc you prove convergence of newton's method with banach

[summer egret]

yea, just worked it out again real quick. Here’s the basic case, multidimensional stuff gets uglier but should ultimately be analogous

Let $f \in C^2(\R)$ (weaker conditions might suffice?) and let $x^$ be a root of $f$. The Newton-Iteration is given by $x^{k+1} = x^k - \frac{f(x^k)}{f'(x^k)}$. In other words, the function $F(x) = x - \frac{f(x)}{f'(x)}$ maps each value to the next step in the Newton-Iteration, and we have $F(x^) = x^*$, so that convergence to a fixed point of $F$ means convergence to a root of $f$.

Further, we have $F'(x) = 1 - \frac{f'(x)^2 - f(x)f''(x)}{f'(x)^2}$, which at $x^$ becomes $0$. Thus, by continuitiy of $F'$, there exists an open interval $I \ni x^$ on which $|F'(x)| < 1$, and therefore $F$ is Lipschitz with Lipschitz constant $<1$ on $I$. Thus, by Banach’s Fixed Point theorem, the iteration $x^{k+1} = F(x^k)$ eventually converges to $x^*$.

[pine jetty]

Sascha Baer:

[summer egret]

@low sandal

[tame trellis]

why not?

I mean

can you show some code maybe or, what you have done

like if you wanted I's three times more often

I, I, I, J, L, O, S, T, Z

should do what you want, right?

just map numbers in a randrange, to the letter distribution?

I don't understand the question, so I'll see myself out, hopefully someone else will reply.

[brave crypt]

hello i am building an LSTM for stock prices (OHLC/volume data)

how would you normalize this data? i'm using to normalizing datasets linearly to [0., 1.]

[spring plinth]

subtract min

divide by max-min

@brave crypt

[brave crypt]

yes but test set is not going to lie within that range

not necessarily

[spring plinth]

ah

[brave crypt]

is that a problem?

[spring plinth]

then german tank problem your way out

[brave crypt]

O_O

explain.

[spring plinth]

okay similar but not exactly that

you can just subtract min - one standard deviation

and then divide by max-min+2 standard deviations

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

similar in nature, but here you are dealing with a continuous variable

@brave crypt worst case if it can be very unbounded, you can consider taking the log of that

[brave crypt]

i might take the log anyways

it makes more sense with stocks

thanks for showing me the german tank problem, i havent heard of that as i don't have a degree yet

this is actually really cool

im reading it with my grandma

[spring plinth]

neither do I

but I just read up on it before

and it might be useful here since you are predicting the maximum of a random variable

[brave crypt]

i'd like to not do so

but rather have another way to represent it to my model

will it just outright throw an error if the input is outside my range

or will it just give erratic responses?

i think i could actually just like

normalize to a reciprocal

or something to that effect?

eh theres no way to do anything like that

i guess i COULD do that with interlacing for representation of values inside/outside the range

no meaningful loss in fidelity

would just reduce precision from 32 bits to 31

or 30 actually due to signed-ness

🤔 that actually does make sense, using signing. i would imagine the model would have a hard time learning patterns from that though

i think it would be useful enough of a result to just clip at +/- some number of standard deviations

or even use a signed reciprocal of stdevs

no, that doesn't work either for sub-1stdev

so i need to use clipping or interlacing

i'll norm to 2 standard deviations and clip there, i think.

[brave crypt]

oh no no no

what i'll do is scale linearly such that mean = 0

and standard deviation = +/- 1

2 stdevs = +/-2

after applying natural logarithm

[brave crypt]

you don't have to normalize if you use batch normalization layers as a cheap trick 🙂

that also takes care of any lookahead bias you can possibly introduce by messing up your preprocessing

you should expect your model to have next to no predictive value

[brave crypt]

hello

I have questions my profs seems to be having trouble answering

[spring plinth]

okay ask

@brave crypt

[brave crypt]

sorry was battling with mathematica

one sec driving home

[neon snow]

Hey! I'm doing a simple program for a course in introduction to programming in Python where I'm trying to find (a) and (b) such that (a^3+b^3=n) where (n) is given. Of course I can write a program that does for example (a=\sqrt[3]{n-b^3}) and then iterates (b=1, 2, 3 ,4,...) until criteria is met when (a^3+b^3=n) but that isn't very satisfactory, or efficient. Any ideas if there's a smoother way to do it?

[pine jetty]

Not Euclid:

[hoary pebble]

I'm thinking this:
a³ + b³ = n
(a + b)(a² - ab + b²) = n

Taking mod n - 1:
(a + b)(a² - ab + b²) = 1 (mod n - 1)

With that, we can search for numbers coprime to n - 1. Once we find them, we figure out which numbers are multiplicative inverses of which (by bezout or by trying each of the numbers together).

So let's say we have those two numbers, x and y. All working in mod n - 1, we know:
a + b = x
a² + b² - ab = y

a + b = x
(a + b)² - 3ab = y

x² - 3a(x - a) = y

3a² - 3ax + x² - y = 0

I'm kinda stuck on solving that (mod n - 1) for a, we're reaching the end of my number theory knowledge

@neon snow

But that should at least hint towards solutions

[neon snow]

@hoary pebble Thank you.

[quaint parrot]

Im having a hard time programming finite difference of heat equation into matlab....idk why :/

[brazen flume]

@neon snow I'm late to the party but have an idea, unless you're already beyond that problem...

As someone who is better at programming than I am at math, here's what I'm thinking:

The largest a or b can be is the cubed root of n, correct?

[median plinth]

No, a or b could be negative

[brazen flume]

Ah shit, well never mind. I was going to say generate a set of all cubes between 0 and the max value of a/b and then iterate through the set, looking if n - (current item in set) is also in set.

Like cubes = {i**3 for i in range(0, max_value_of_a_or_b)} in Python. But never mind, should've thought about it more.

[topaz garden]

Why is this true? Particularly the equality between the union of two intervals and the interval from negative infinity to positive infinity?

from https://en.wikipedia.org/wiki/Interval_arithmetic#Simple_arithmetic

Is it a mistake?

[median plinth]

The first equality makes the sense

The second is uh

[topaz garden]

Yea...

[tardy current]

I think I understand why they wrote that

they are replacing the union with their previously defined addition of intervals

of course the equality isn't actually true, but I believe that's what they did

they probably only care about bounds in one direction and wanted an interval as the answer

but they should have written subset

[topaz garden]

Well actually...

If you remove the intermediate inequality and just focus on $\frac{1}{[y_1,y_1]}=(-\infty,\infty)$

[pine jetty]

Icohedron:

[topaz garden]

Because 0 is in the interval, and we're interested in the range of values that the fraction takes

We'll have 1/0 in there

Which can be negative or positive infinity

[tardy current]

yes there is no question that you can get large numbers of either sign

the problem with the equality they stated is only that if y1 and y2 are small in magnitude, you can not get numbers close to zero as the result of division

hence the two separate intervals