numerical-analysis

[violet sundial]

Yes to R^3 and convex polyhedra. By a lot of work, I'm mostly thinking that it means calculating each starting normal (though I suppose those are constant once the polyhedron is specified), making note of that, and then a matrix multiplication plus calculating the angle from the y-axis for each normal once the polyhedron is rotated.

[wide spear]

Calculating the normal that is the most vertical is simple, take the dot product with the upwards unit vector and keep track of which one gives the largest value

[fallen relic]

Function fzero on Matlab doesn't work for finding a root of $e^x-x-1$ what can be done?

[pine jetty]

Suraj

[wide spear]

#computing-software

[wide spear]

Computationally, I don't think this is actually too bad

If you are using blas for the mat mul

If you have already implemented it and it's too slow, then ok we can discuss further

[violet sundial]

That's fair. I was just working through how to implement it and it seemed like a fairly brute force method, so I thought it was worth checking into it.

[wide spear]

There are probably some techniques using the symmetry group of the polyhedron if your polyhedron has additional symmetry

Hmmmm

There are probably some clever things that you can do with matrix decompositions (https://en.wikipedia.org/wiki/Rotation_matrix#Decompositions)

But probably a bit overkill for a 3x3

[violet sundial]

Thanks for talking through it with me. Much appreciated.

[wide spear]

Let me know if the brute force technique is not satisfactory

[pine jetty]

Farm

[wide spear]

Schrodinger equation in spherical coordinates?

If you aren't doing this numerically, this isn't the correct channel

It looks like you're trying to do separation of variables, so #odes-and-pdes

[neat haven]

Okay

[molten knot]

when you make the ansatz $u(r, \phi, \theta)=R(r)P(\phi)T(\theta)$ and plug it in, distribute the derivatives, etc you’ll get a lot of cancellation

[pine jetty]

ansh

[molten knot]

bc like say d/dr P(\phi)=0 since P has no dependence on radius

[wide spear]

Bad

[neat haven]

I'll try it out

[wide spear]

:spray bottle emoji:

[wary oriole]

@Farm This is solved in details in all QM books. better read than ask here cause the whole derivation is a lot. As my prof. said, take the whole afternoon 😄

[pine jetty]

ArgR4N

[wide spear]

What have you tried

[meager dome]

We have seen a proof, the difficult part was the inequality of $1/Cond(A)_N \geq \frac{dist(B \textit{ singular})}{N(A)}$

but I am not very satisfied with some intuitively justified steps

[pine jetty]

ArgR4N

[meager dome]

For this we constructed a B singular such thath $N(A - B)/N(A) = 1/Cond(A)_N$. I don´t really thing this proof is at my level for me to think a proof 😅, but I wanted to see a more rigourus one. (The unjustified steps were related to proving the existence of some vector with some condition relative to the norm. This condition is easily seen by taking some known norm, but what is difficult is to see that such a vector exists for any norm.)

[pine jetty]

ArgR4N

[next garden]

This text is confusing me. In particualr when it talks about divided differences. The value $k_{ij}$ should be a scalar. But if I blindly do a difference of the normals divided by the edge length I get a vector. I think this is not properly explained.

[pine jetty]

Makogan

[wide spear]

Can you share the entire document

[outer jacinth]

does anyone know what happnes to L if you divede lets say A R1 by 6 to get the postiont 1,1 to be 1 so its eaiser to the rest of lu decomp

[wide spear]

What

[outer jacinth]

like say i divied row 1 of matrix A by 6 what would happen to L

[cyan wing]

I have a pretty silly question, but it just involves a standard numerical way to compute derivatives at some point x

So my professor wanted us to find the gradient at some given points and he told us to compute them numerically given a tolerance of eps = 10^-7

So i just wanted to know how this method

$f^\prime_{x_i}=\frac{f(x+\epsilon e_i)-f(x-\epsilon e_i)}{2\epsilon}$

[pine jetty]

till

[cyan wing]

Would be different from just repeatedly going through the estimates by having smaller and smaller values of a step size h and just stopping when the difference between two consecutive estimates would be less than epsilon

[wanton dawn]

what are " the estimates" you are going through here?

[cyan wing]

oh i mean just the approximations from the central difference formula D_k

my code for that in matlab would be just a while loop iterating through smaller and smaller values of h until $|D_n - D_{n-1}|<\epsilon$

[pine jetty]

till

[cyan wing]

hmm i dont know if my question's making sense sorry. purely just something i was curious abt

[vapid plume]

It does make sense

If you know what f is, you can use f'(x) and then use | d_n - f'(x) | < e for the tolerance

But if we don't know what f is or f' is difficult to compute then I believe we would use | d_n - d_n-1 |

[wide spear]

Another alternative is to instead compare derivative approximations of two different orders

[wide spear]

This is a second order approximation

So you can compare it with a 4th order approximation with the same epsilon

[cyan wing]

oh right that would make sense

[wide spear]

This is the idea behind adaptive runge kutta methods

[old compass]

How would I find the rate of convergence of f(x,n)?. I already was able to get rid of the absolute value signs because it alternates being above and below W0(x). It appears to be 1/W0(x).

[wide spear]

Have you tried a taylor series

[old compass]

nope, but how would that help me? and do you mean the taylor series for W_0 or f?

[wide spear]

For f

Taylor expand f around W_0 in terms of n

[old compass]

oh, that makes sense

wait how would I do it around W_0 if W_0 is not an integer?

and n is?

[wide spear]

Hmmmm

Actually I think it's probably easier to consider a continuous version

The continuous version of this is $y'(x,t)=\ln(x/y(x,t))$

[pine jetty]

Angetenar

[old compass]

how is that continuous and how do you interpret that?

[wide spear]

Continuous as in it's an ODE and not a recurrence relation

Hmmmmm

[old compass]

oh

that makes sense

[wide spear]

No I don't think this is the right way either

So the solution is $W_0=xe^{-W(x)}$ right

[pine jetty]

Angetenar

[wide spear]

Where W(x) is the product log function

Or the limit I should say

[old compass]

no the limit is 1/W(x)

i think

thats what the graph implies

[wide spear]

Well, mathematica is telling me that the solution to y=log(x/y) is y=xe^(-W(x))

Anyways

We're doing to do a perturbation analysis

[old compass]

oh

yah

I thought you mean to my original question

[wide spear]

So for the solution W_0, we are going to perturb this slightly

So we have W_0+eps

And then we see what happens to it

So let's say that f(x,n-1)=W_0+eps, we want to know what W_0 is

So we compute the taylor expansion of $\log(x/(W_0+\eps))$ with respect to $\eps$ around $\eps=0$

[pine jetty]

Angetenar

[wide spear]

To 3 terms, this is $\log(x/(W_0+\eps))\approx \log(x/W_0)-\eps/W_0+\eps^2/(2W_0^2)+O(\eps^3)$

[pine jetty]

Angetenar

[old compass]

that makes sense

[wide spear]

Then, when we plug this into the limit, we have [\frac{W_0(x)-f(x,n)}{W_0(x)-f(x,n-1)}=\frac{1}{\eps}\left(W_0-\log(x/W_0)-\frac{\eps}{W_0}+\frac{\eps^2}{2W_0^2}+O(\eps^3)\right)]

[pine jetty]

Angetenar

[wide spear]

Now, we know that W_0=log(x/W_0) so those two terms cancel

And we are left with $\frac{1}{\eps}\left(-\frac{\eps}{W_0}+\frac{\eps^2}{2W_0^2}+O(\eps^3)\right)=-\frac{1}{W_0}+\frac{\eps}{2W_0^2}+O(\eps^2)$

[pine jetty]

Angetenar

[wide spear]

So the limit is -1/W_0

[old compass]

ty

that makes sense

[wide spear]

Oh wait there is a sign error

It should be 1/W_0

Hopefully that made sense

[old compass]

yah it did

ty

[wide spear]

Perturbation analysis is a useful and powerful technique

[old compass]

I see that now

[sacred mural]

I'm a Maths Undergrad Student and intrested in the field of Numerical analysis, I select the topic for my undergrad thesis "Numerical solutions for unsteady laminar boundary layer flow and heat transfer over a horizontal sheet with radiation and nonuniform heat Source/Sink". When i read this Paper there are a lot of topics from Complex fluid dynamics and only numerical technique has been applied to solve governing equation. My question is that is this good topic for Numerical analysis field or should i change it, also consider future scope. Thanks

[wide spear]

Seems like a reasonable undergrad thesis topic

[crisp lake]

Hmm. I suspect there should be theorems about convergence of piecewise rational approximants but am hard-pressed to turn up the references or rates of convergence or bounds on the constants involved.

[wide spear]

Have you looked at padé approximants?

[crisp lake]

Each interior piece of a C^d piecewise rational clearly has to be a 2-point Padé approximant, but I don't know what that says esp. about rates of esp. uniform convergence.

[wide spear]

Convergence with respect to what? Shrinking interval size?

[crisp lake]

L^p norms I'm hoping, even if L^1 or L^\infty .

[wide spear]

That's the norm that you are interested in convergence in

What is changing to cause the convergence? Degree of approximant or interval size?

At any rate, for each piece of the interpolant, you should have an error bound that depends on the degree of the interpolant and the size of the interval

You just take all of these errors together

[crisp lake]

The degree types would stay constant. More knots & piecewise intervals would be what increases, like piecewise polynomials of having a fixed degree within each piece.

It appears that lower degree types make more sense to avoid accidentally introducing singularities on the real line & whatever the stiffness issues are with high-degree polynomials & rationals.

[wide spear]

Have a look through this

[crisp lake]

The missing piece may be the piecewise rational aspect of it.

I can't actually find the analogous convergence results about piecewise polynomials / splines.

[wide spear]

You can derive it then

[crisp lake]

There should be similar results with convergence as the piecewise intervals proliferate for plain old splines / piecewise polynomials, too, but I somehow don't know where to look for those, either.

[wide spear]

Google scholar?

[crisp lake]

I appear to be a particularly bad prompt engineer.

I've already been throwing queries at Google Scholar for a few days. Somehow ”piecewise rational» turns up very little.

Jan Gelfling or similar wrote something in 1975 & nothing else has touched the subject since, or so it appears.

[wide spear]

https://www.sciencedirect.com/science/article/pii/S0021904502000606 ?

Perhaps searching piecewise rational approximation would be better

[crisp lake]

n-term rational is a single grand rational that's got some sort of duality or homomorphic or diffeomorphic aspect between the two. So a bunch of things will go about chattering about piecewise polynomials & then will have a grand unified rational as the thing they're establishing some sort of correspondence with.

So that's a sort of search term overlap that's a total nightmare to prompt engineer my way past. And I've thus far been unable to prompt engineer my way past it.

[scarlet tide]

What’s a prompt engineer?

[crisp lake]

It's an inflated title for someone who issues commands to AI systems.

[crisp lake]

I was looking for the 1D case but the paper looks close to applicable.

[crisp lake]

How does one say that they're representing a polynomial space in a particular basis? e.g. Chebyshev, Bernstein etc.? What are degree d polynomials in one of those bases?

[wide spear]

What do you mean

Not quite sure what you mean for the first question

[crisp lake]

What are the names for those spaces?

[wide spear]

For the second, if you are asking about representations of monomials x^d, then just take inner products

[wide spear]

It's just polynomials of degree d

The basis doesn't change the vector space

[crisp lake]

That's why it's a distinct concept from a vector space, yes.

[wide spear]

I have no clue what spaces you are referring to then

[crisp lake]

Neither do I, hence my asking.

[wide spear]

What space do you want a name for

Polynomials with specific choices of bases do not have names

[crisp lake]

There's more & different structure, like Fourier-Plancherel bits where they're a sort of leading fragment of a larger space, different bases corresponding to different weights and inner products etc.

[next garden]

Is there a technique to sample a surface, similar to a poisson sampling, but instead of having a possoin disk aorund each sample, you have a poissin ellipsoid that follows the curvature of the surface?

I am trying to cover an implicit surface with samplse and trying to get it to match curvature directions.

[crisp lake]

Upon closer examination, it suffers from the same problem that makes it all ungooglable: it establishes a duality between piecewise polynomials and one big grand rational function instead of ever even mentioning piecewise rational functions.

[crisp lake]

Hmm. Maybe this terminology bit is an abstract algebra question.

[brittle breach]

State Lagrange's interpolation and show that Newton's forward formula is a special case of it

[vapid plume]

Hello.

I am interested in pursuing graduate studies in numerical analysis and scientific computing. Particularly, finite element methods, numerical linear algebra.

I have some background in numerical analysis, finite difference methods, and geometric multigrid.

As for numerical linear algebra, I have not had a proper course or self-studied it.

I was wondering if I could have a rough, or possibly accurate, guide on how to prepare for graduate studies as a start, as in which books would be suitable to get familiar with finite element methods, the theory with functional analysis, and an introduction to numerical linear algebra.

Thank you.

[wide spear]

For numerical linear algebra, Trefethen and Bau or Demmel are standard books

[vapid plume]

For finite elements and its theory?

I havent had a proper course in functional analysis as well, I was thinking of beginning with Kreyszig or is it too introductory?

[wide spear]

You don't need that much functional analysis to do numerics

There are many many books about finite elements

Iserles has a book on numerical methods for differential equations which covers finite/spectral element methods

[vapid plume]

Thank you.

[crisp lake]

What is this?

@wide spear Jan Gelfgren had two manuscripts in the 70s & it appears the one titled for multi-point Padé approximants has exactly what I'm looking for.

[wide spear]

Nice

[crisp lake]

This is remarkably ungooglable. I think it's mostly a coincidence that I ever turned up anything at all because all of the search engines get totally derailed by the numerous papers etc. mentioning the duality between splines on subdivisions of an interval & one grand rational spanning the whole unsubdivided interval.

[crisp lake]

There's some kind of arrangement with $L_{\overrightarrow{d}}(x)=\prod_k\frac{(x-x_k)^{d_k}}{d_k!}$ that's supposed to be involved with a by-hand construction of interpolatory Lagrange-Hermite polynomial bases, but somehow I'm fumbling the details & failing at Google all over again.

[pine jetty]

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

[wide spear]

https://epubs.siam.org/doi/pdf/10.1137/S0036144502417715

Does this explain

[crisp lake]

I think it only mentions things once.

[halcyon thunder]

is numerical analysis mostly for pde’s or does it often show up in other problems as well?

[wide spear]

Numerical linear algebra is another big subfield

Also strictly speaking, numerical methods for odes are not pdes

And then you also want to do things like integration

And interpolation

[crisp lake]

There are numerical analysis things for special functions, solving nonlinear equations, numerical differentiation exists, stochastic ODEs & PDEs & a wide variety of statistical things fall within its purview, there are integral equations, KKT optimisation, variational optimisation even apart from Euler-Lagrange & Hamilton etc. equations, also discrete optimisation, linear programming is enormous & it keeps going on.

[sacred mural]

I am an undergraduate mathematics student interested in pursuing a thesis in numerical analysis. Initially, I chose a topic from fluid mechanics, aiming to apply numerical techniques, but my supervisor’s expertise lies in transforms, not fluid mechanics. Since my interest is in numerical analysis while my supervisor specializes in transforms, I am unsure whether to change my topic within the numerical field to align with my supervisor’s expertise or to continue with my supervisor's field transforms If I choose transforms or numerical analysis, what trending and suitable topics are? Alternatively, if I stick to numerical analysis outside transforms, what topics are currently relevant and suitable for an undergraduate thesis?

[heady mesa]

Lots of stuff to look at for both.

Transforms pair well with Spectral methods, which is undergrad level and you can find lots of good resources for them.

Outside idrk

But there is always lots of stuff, it's a massive field

[molten knot]

scott-brenner mathematical theory of finite elements is the standard book

i think if you really want to understand fem in terms of error estimates and analysis, you would need a bit of background in sobolev spaces (brezzis is a good book)

but for implementation purposes it’s not needed too much

[vast nexus]

How do I solve this?

I don't really understand the formula

What do I plug in for a? And what is p equal to?

[wide spear]

p is the limit of p_n

a is what you are trying to find

[vast nexus]

Ohh, how do I find it?

[wide spear]

By looking at the structure of the limit

[vast nexus]

is it something like this?

[wide spear]

I feel like you're overthinking this

[vast nexus]

hmmm most likely

[crisp lake]

Hmmm. There's a weird thing where Gelfgren's error bound is pertinent to a complex domain containing the interval.

[molten knot]

p=lim 1/n^2 so you just plug in
pn=1/n^2, p=1/n and test a=1,2,etc until you stop getting a constant for lambda

[crisp lake]

Hmm. Comparing convergence behaviour of $C^d\left[x_0,x_0+m\cdot\Delta x\right]$ piecewise polynomials & rationals with knots at $\left{x_0+j\cdot\Delta x\right}_{j=0}^m$ where each polynomial piece has degree $2d-1$ & each rational piece degree type $(d,d-1)$ is awkward for me. Gelfgren's 1975 paper is tough for me to decipher.

[pine jetty]

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

[worthy geyser]

could someone point out where i've gone wrong in this 4-point ifft flow diagram?

[wide spear]

What is the answer supposed to be

[worthy geyser]

According to https://scistatcalc.blogspot.com/2013/12/fft-calculator.html the answer is {0.175 + 0.175i, 0.175 + 0.375i, -0.075 + 0.625i, 0.125 - 0.175i}

[wide spear]

Is that using the same normalization

[worthy geyser]

normalization?

you mean the division at the end? or the twiddle factor? either way it should be the same by definition of idft

psi(0) = 0.7/4 + 0.7/4i is consistent with the online calculator. But then it goes to shit with the others

[wide spear]

I will double check your math in a bit

[wide spear]

Hmmm I tried computing the FFT myself and got a third answer

😵‍💫

[worthy geyser]

I took a lil break and got back to it. Turns I forgot to do the bit reversal even though I had it written down

[wide spear]

Ah ok

[worthy geyser]

so omega(1) and omega(2) should be swapped

[wide spear]

Well I did have omega(1) and omega(2) swapped

And I did not get the correct answer

[worthy geyser]

so you did omega(1) = 0.8 - 0.5i?

and then omega(2) = -0.2 + 0.6i?

[worthy geyser]

Also this is inverse not forwards

[wide spear]

Yes doing the inverse

[wide spear]

Yes and yes

[worthy geyser]

what do you have for your twiddle factors>

second stage that is

[wide spear]

Anyways you should still do it yourself with the swapped omegas

[worthy geyser]

that's what im doing now

[wide spear]

It's just 1s and i's

[worthy geyser]

I got all but 1 imaginary part the same 🥲

[wide spear]

😵‍💫

[worthy geyser]

nevermind... I just can't read my own writing

I need to invest in a tablet asap

FFT is unbelievably genius work

wish i came up with it

[crisp lake]

Gauß was good.

[unreal birch]

Gaub

[crisp lake]

ß is a ligatured form of sz, expanded as ss in modern German.

I don't think I even understand how to derive polynomial or piecewise polynomial bounds or rational bounds without the rational being piecewise.

Chebyshev polynomials have nice DCT algorithms for multiplication.

[crisp lake]

Piecewise constants seem to use the mean value theorem.

[crisp lake]

Hmm. The references I can whip up for whatever reason presume $C^d([x_0,x_N])$ spline interpolates are interpolating $d$ function derivatives vs. the derivatives just being continuity conditions & only the 0th-order derivatives / actual function values interpolating anything. At that point, I'm clueless as to what the errors of spline interpolation are & they don't seem to be dealing with a different problem from 2-point $d$-order Hermite interpolation repetitively packed end-to-end.

[pine jetty]

Nadia

[crisp lake]

I'm reading the Oliveira-Takahashi paper on ITP and am not convinced it's really doing useful things.

[crisp lake]

Of course one can just compare Halley, Newton, secant, and bisection probe points and beyond that, and since you went through the trouble of evaluating the function at all these points, you can examine all of those points that are admissible and pick the tightest opposite sign bracketing interval. It also had the not so small issue of not verifying satisfaction of residual tolerance or signalling the failure to satisfy such.

[crisp lake]

I can’t figure out how to extend the Rolle’s theorem argument for error bounds on a Lagrange interpolant to (C^d) piecewise polynomials or rationals.

[crisp lake]

I’m pretty sure there are known linear combinations of $\prod_{0\leq k\leq K}\frac{\left(x-x_k\right)^{m_k}}{m_k!}$ where $\left(\forall k\right)\ 0\leq m_k\leq d_k$ to produce a basis for Hermite interpolants. I think the only relevant ones have $m_k=d_k+1$ for all but one index. But that still doesn’t really help estimate errors for $C^d$ piecewise polynomial or rational interpolants.

[pine jetty]

Nadia

[meager sky]

Hi guys, can you help me to do this exercise using this algorithm 7.1, I've been doing this over and over and always got the coefficient zero.

[wanton dawn]

Set up a linear system by requiring the stencil to be exact on polynomials 1, x, x^2, ..., x^d, where d is as high as you can go. The error big O follows from this and Taylor theorem

[meager sky]

thank you for the ideaa

[crisp lake]

I'm still having trouble figuring out some way to get bounds on the errors for rational interpolants that are of more than minimal degree with the other degrees of freedom being used for e.g. continuity conditions with neighbouring intervals. I can't figure out how to get bounds on rational interpolants period. Or why do they even have to be interpolants? What if I do something like LOESS with C^d piecewise rationals? I can't figure out how to get any meaningful bounds on really any problems of the kind, not even for polynomials or anything beyond piecewise linear.

Is there a different / better channel to get ideas about these things e.g. analysis channels?

[hot otter]

hello

[hot otter]

wow what a nice problem you read that again i am sure i am confused i have no idea what you talking about but you should make your problem specific meaning make different different piece out of it make it descrete then try to solve it

[crisp lake]

Alright, let's try an example.

Suppose one has a rational $r(x)=\frac{p(x)}{q(x)}$ of degree type (3,2) on an interval $[a,b]$ so that it would interpolate a function $f(x)$ at $a$ and $b$. Its other degrees of freedom would be absorbed by setting $r'(x)$ an $r''(x)$ at $a$ and $b$, presumably with the usual implicit differentiation ideas involving constraining coefficients of $p(x)$ and $q(x)$ by considering $p^{(k)}(x)=\frac{d^k}{dx^k}\left(q(x)\cdot r(x)\right)$. What kinds of bounds can I get on $\left\vert f(x)-r(x)\right\vert$ on $[a,b]$ in these sorts of situations?

[pine jetty]

Nadia

[crisp lake]

@hot otter Does that clarify what I mean?

I chose degree type (3,2) thinking that it's sort of like the piecewise linear case for rationals.

[crisp lake]

Does trying to do KKT optimisation even make sense?

How even are strict inequalities dealt with for KKT?

[wide spear]

#optimization

[crisp lake]

I think the answer is that optimisation is not how one derives error bounds so it does not make sense.

[crisp lake]

I have some vague idea that $\left\vert R(x)\right\vert$ where $R(x)=f(x)-r(x)$ that gets extremised at $\xi\in(a,b)$ might have $R(x)=R(\xi)+\sum_{k=0}^\infty(f_k-r_k)\cdot\frac{(x-\xi)^k}{k!}$ and then getting from that to any sort of upper bound on $|R(x)|$ mystifies me, where $f_k, r_k$ are the Taylor coefficients about $\xi$ of $f(x)$ and $r(x)$. At the point that leaves me, it doesn't even get a bound in terms of the function.

[pine jetty]

Nadia

[crisp lake]

Milne-Thompson (1933) has something on error bounds for Thiele interpolants but it's a little ugly.

[crisp lake]

Basically the denominator of the interpolant is still part of the bound, something like:
$\frac{\prod_{k=1}^n(x-x_k)}{n!\cdot\left(q(\xi)\right)^2}\cdot\frac{d^n}{d\xi^n}\left(f(\xi)\cdot\left(q(\xi)\right)^2\right)$ or similar.

[pine jetty]

Nadia

[lusty elm]

This question is regarding numerical stability and higher order (convergence) methods.
I have implemented a Runge-Kutta 4 method as part of an N-Body simulator and have measured the convergence order empirically but it appears to be order 5 instead of the theoretical 4.
I can't seem to figure out why or even if this is theoretically plausible. The truncation error is O(h^5) so that may be related but I can't piece it together
I can provide my data if needed
My Explicit Euler and Leapfrog methods are yielding the expected orders 1 and 2 but I can't explain the one for RK4

[wide spear]

Funky

[stone drum]

if the method is 'stable' (i. e. small perturbations in input cause small changes in output) then the truncation error carries over to the convergence error, which should explain why the order is 5

could you provide more details about the method, e. g. specific weights?

in general, there are (implicit) runge kutta method with higher convergence order than number of points, see gauß, radau or lobatto formulas

[lusty elm]

I think I'm using the standard explicit RK4 method computed as 1/6 * (k1 + 2*k2 + 2*k3 + k4)

This is the graph I've been able to generate (I know its information dense, page count limit is ridiculous so I've had to make compromises 😅 )
Also I'm a CS student so a lot of the in-depth maths went over my head

[stone drum]

well, superconvergence does happen sometimes and it depends on the regularity of your differential equation. but I don't know if I can say much about it here

[lusty elm]

hmmm

would you say my implementation has an issue then? I don't really know about regularity.

[stone drum]

hard to say, but I don't think there's a bug or anything since the approximations seem to converge just fine

[lusty elm]

im conflicted on whether I should keep this rk4 in my report or not since I can't seem to explain it

I might just email my professor for help since all of this was entirely optional

[stone drum]

yeah, an email sounds good

[wide spear]

Very unlikely

You're doing n body right

Try doing it with a 2 body system

You should be able to perform a theoretical analysis directly for the 2 body solution

[lusty elm]

this is how I have arrived at my convergence measurements

[wide spear]

Oh

Then you should be able to do a theoretical analysis for the 2 body system in particular

Some terms may magically cancel

[lusty elm]

I took the simple example of two bodies of mass 4, separated by 2 units with velocity 1 in opposite directions to form a circular orbit

[lusty elm]

i do not know how to do this

[wide spear]

Try doing it for a noncircular solution

[lusty elm]

whats the intuition for that

[wide spear]

When you have very symmetric things, you get better behavior

[lusty elm]

such as?

the two cases i can think of is either bodies which just oscillate or circular

[wide spear]

I suspect that the symmetry of the solution is improving your solution accuracy

Can't you do an elliptical solution

Alternatively, do it for a n body system where n>2

[lusty elm]

how would i measure convergence for that? dont i need to compare error against analytical solution?

[wide spear]

Then the reference solution, you compute using something like RK5 or RK6 or a higher order method

[lusty elm]

oh

hmm

might just need to leave this as an exercise for the reader (my professor)

[cold wasp]

Hey

[magic spruce]

why does something need to be well posed/well conditioned to obtain numerical stability. im trying to find an intuitive explanation or proof for this

[wide spear]

Presumably you mean that the underlying problem needs to be well posed for a numerical solution to be stable?

[magic spruce]

yes

[wide spear]

You are familiar with floating point right

Essentially, whenever you are solving a problem numerically, you can only approximately solve the problem

So you find the approximate output for an approximate input

If the underlying problem is not well conditioned, then even for two inputs that are close together, the solutions can be very different

In which case, your numerical solution is not of any use

[wanton dawn]

See machine epsilon: https://en.m.wikipedia.org/wiki/Machine_epsilon. Double precision is most common in non deep learning computing. This says that when approximating a real number by a floating point number, you may incur relative error up to machine epsilon. Thus if your mathematical function amplifies relative errors in the input numbers, your numerical implementation of it may have large relative error.

[wary oriole]

I have a stupid question. so I did fea solver for helmholtz equation. and now trying to compare it to fundamental solution. for fundamental solution do I have to numerically solve the roots? It doesn't make sense to me cause Im trying to compare my numerical solution to an exact solution but for exact solution I have to numerically solve it? than Im compare numerical solution to numerical solution? If I root finding there will be errors, so my comparison is not exact right?

I'm a bit confused here.

[wide spear]

Where is the root finding coming in

[wary oriole]

i need fundamental solution's roots. my goal is to calculate the wavelength of two solutions

btw, how should I calculate the wavelength since the wavelength are different further out?

$F(r)=\frac{i}{4}H_0^{(1)}(kr)$

[pine jetty]

yehuihe

[wary oriole]

for 2d

also how do i get roots of my numerical solution? it's just 200 points

[deft sedge]

Let's say I have a succession of values $x_{n}$ which tends towards $z$. The absolute error is majored such that $ \z-x_{n+1}| \leq | x_{n+1}-x_n |$. I'm trying to find a similar expression but for the relative error. can someone help?

[pine jetty]

lazynjnja

[wide spear]

Divide both sides by z

[deft sedge]

yeah but i dont know z

[deft sedge]

... because I'm told to use this inequation so that the program halts when the relative error becomes less than epsilon, so I probably need to find some bound for the relative error that I can compute

[wide spear]

Ok

If you have an iteration scheme and are deciding when to stop there's a few things that you can do

Using the absolute value abs(x_{n+1}-x_n) of the update isn't uncommon

You can also use the relative value of the update

So abs(x_{n+1}-x_n)/abs(x_n)

[mystic jacinth]

Would anyone be able to provide me with resources that have butcher tableu's for Dormand-Prince methods? I currently use a RKF4(5) method to integrate my system as it's the only one ive been able to find sufficient material on. I also have no clue how i'd get the truncation error from just the butcher tableu and if how its used changes between what I currently have implemented to find the new timestep if the error is above/too far below my tolerance

[mystic jacinth]

https://www.uni-muenster.de/imperia/md/content/physik_tp/lectures/ss2016/num_methods_ii/rkm.pdf seems to go over the process of finding the optimal timestep :D

pg 9-11, im not sure if that applies to all methods?

[wide spear]

Is the dormand prince paper not satisfactory

[mystic jacinth]

sometimes I feel like giving up 😭 I was searching dormand-prince methods , dormand prince numerical integration, the list goes on, but just searching 'dormand prince' gave me what I needed

[stone hill]

i am confused on conditions for this to converge

so far: p(M) < 1 where x0 and b are eigenvectors for M and, b = 0 when M is idempotent are the only ones i have gotten so far

I was told that these conditions are not enough to cover all cases in a help channel given that eigenvalues could be complex or smth, so i do not know what to put

additionally, there is a second part that asks how M must be defined for the iterations to converge to a solution to Ax=b. with some basic algebra I got

x = Mx+b
x-Mx=b
(I-M)x=b

therefore I-M=A => M=I-A

is this right or am i forgetting smth

[mint helm]

Has anyone seen functions that graph like this? What field of study is this sort of thing? (these are all the same class of function with different class parameters over the same input range. I've never seen this before).

[stone hill]

From some googling, I found that the only conditions necessary for convergence is that p(M) < 1 and that's it. Why is that?

Sorry for skipping your message

[rancid smelt]

the easiest way should be using Banach fixed-point theorem
you have to show that $\varphi(x) = A x + b$ is a contraction, e.g.
$|| \varphi(x) - \varphi(y)||_2 \leq k ||x-y||_2$ with $k < 1$
In your case:
$|| Ax -b - (Ay-b) ||_2 = ||A(x-y)||_2 \leq ||A|| ||x-y||_2$

@stone hill and if $\rho(A) < 1$ then $\varphi$ is an contraction and you can apply Banach fixed-point theorem

[pine jetty]

malakai

malakai

[stone hill]

Oh i see. I guess i overthought it. Would not have thought of using that theorem

Or actually more like under thought it

I feel bad for being frustrated with that guy last night now lol, but to be fair i did not have any of this pre req knowledge.

[rancid smelt]

with more experience it gets easier

[stone hill]

Ok so it turns out i overthought the conditions for convergence even

[wide spear]

[stone hill]

all i had to do was x* = Mx*+b => Ix*-Mx*=b => (I-M)x*=b => x*=(I-M)^-1b$

[mint helm]

it looks like logarithmic curves that start non-differentiable and smoothly become differentiable

anyone seen something like that before?

[rancid smelt]

that's not true. Look at the following example:
$ x_{n+1} = \begin{pmatrix} 5 & 0 &0 \ 0 & 5&0\ 0&0&5 \end{pmatrix} x_n$
An unique solution exists x = 0.
However your iterative alg. will not converge

[pine jetty]

malakai

[rancid smelt]

for $x_0 \neq 0$

[pine jetty]

malakai

[stone hill]

x* is defined as whatever the thing converges to, but that's true oof. So i wasn't overthinking at first 😭

Oh well :(

I guess I got that one partially wrong

[wide spear]

If it converges it converges to x*

[stone hill]

Yeah

but they're right in that that would break my x* definition

Oh wait no

It wouldnt

Cuz if b is zero then it would just make x* zero and that would work actually

[wide spear]

?

[wide spear]

No the problem with this is that it sometimes implies that convergence happens when it doesn't

[rancid smelt]

you have differentiate between the solution of the problem x*=0 itself and the convergence of your algorithm. Just because your problem has a solution does not mean that the algorithm will find the solution.

[stone hill]

Oh ok

[stone hill]

Doesn't the spectral radius condition already cover this?

The sequence will converge in that situation but it's a kinda pointless convergence

[wide spear]

Well it diverges if x_0 is nonzero

[rancid smelt]

yeah the spectral radius condition covers the convergence

[torpid sandal]

this is an interesting question, because you can just linearize it to see the answer

[stone flax]

how do you guys get to this point?

how much time have you put in math?

[torpid sandal]

$\begin{bmatrix}x_{n+1}\1\end{bmatrix} = \begin{bmatrix}M&b\0&1\end{bmatrix} \begin{bmatrix}x_n\1\end{bmatrix}$

[pine jetty]

Trivial Lemma

[stone hill]

The skill for this channel really varies

[torpid sandal]

this reduces the problem to "when does x_(n+1) = Ax_n converge"

hope this helps @stone hill

[stone hill]

Oh I already submitted

It was a really dumb answer

[torpid sandal]

ah, sorry just looked at it 😅

[stone flax]

why is the value of e about 2.7?

[stone hill]

Nah i think it's cool anyway

[torpid sandal]

because 1 +1 + 1/2 + 1/6 + 1/24 is about 2.7

[torpid sandal]

I think you have convergence if the matrix norm is < 1

[stone hill]

We never learned with matrix norm was

[stone flax]

you in harvard or smth

[stone hill]

But I'm curious what this would be useful for

[torpid sandal]

which in this case doesn't work because it will probably be 1 now that I look at it

[torpid sandal]

you can use it to find a closed form for x_n

by diagonalization

$\begin{bmatrix}x_{n}\1\end{bmatrix} = \begin{bmatrix}M&b\0&1\end{bmatrix}^n \begin{bmatrix}x_0\1\end{bmatrix}$

[pine jetty]

Trivial Lemma

[stone flax]

@torpid sandal what was your GPS

GPA

[torpid sandal]

I'm not American

I saw the equivalent GPA score for my bachelors after I graduated

and it was about 3.9 iirc

[stone flax]

whats your chess elo

?

[torpid sandal]

I haven't gone to tournaments since High school

it was about 1700 fide

[stone flax]

what opening should i play?

@torpid sandal

[torpid sandal]

idk

I found openings boring

[stone hill]

Yeah honestly I dunno where I'd go from here. But it's neat thinking there's a ton of ways to go about this

[torpid sandal]

so I just played what I knew

[torpid sandal]

Do you know when x_(n+1) = Ax_n converges? for A a matrix?

[stone hill]

p(A)<1 right?

[torpid sandal]

that does ensure convergence yes

and it converges to 0

[stone hill]

Oh there's more?

[stone flax]

are the millennium problems solvable?

[torpid sandal]

in this case the spectral radius is at least 1

however, having an eigenvalue of 1 isn't an issue

[stone hill]

Ig it also converges if A is idempotent

[torpid sandal]

your thing will converge if p(M) < 1 (in general)

regardless of b

there are other solutions too

[stone hill]

Like a condition that would cover more cases?

[torpid sandal]

Yeah

if the generalized eigenspace corresponding to 1 is equal to the eigenspace and all the other eigenvalues of M are < 1 (in absolute value)

this will be regardless of x_0

depending on x_0 there are several conditions

[stone hill]

Why would this guarantee it?

[torpid sandal]

you can essentially separate x_n in a part there and a part outside of it

and look at two separate sequences

[stone hill]

I'll try thinking ab it

[lost cedar]

I'm having a hard time understanding maximum possible error. I've found the Taylor approximations for t_1(x) and t_3(x), but I'm not sure how to move forward. I know the setup is max -1<=x<=1 | f(x) - t_n(x) |, but how do you move forward? I'm using Elementary Numerical Analysis by Kendall Atkinson Weiminhan, followed one of the examples in the book so looks like they just plugged in 1 into f(x) and t_n(x), but I'm gonna assume that's not always true.

[torpid sandal]

prove that t1 is an upper bound (on (0,1))

and t3 is a lower bound

i.e. tan^(-1)(x) < x for 0 <= x <= M (and M >= 1)

[mint helm]

Nobody? Is this like, a new class of curves? I have simple functions that generate these. Is that interesting?

I'm not involved in academia in any way, so, a bit hard to gauge if I've got something interesting.

[fair forum]

on the power method, am I correct in saying that the error is bounded by O(|λ2/λ1|^k), and so the rate of convergence is O(|λ2/λ1). and therefore we would have for the inverse power method that error is bounded by O(|λn/λ(n-1)|^k) and so the rate of convergence is O(|λn/λ(n-1)) ?

[rancid smelt]

It looks like a series of step functions that converges to a root function as the step width gets smaller and smaller. The interesting question is what the generating function looks like.

[rancid smelt]

you consider $g_i(x) := f(x) - t_i(x), i=1,3$ and then you search for local extrema in ]-1, 1[ and additionally, you have to consider the existence of boundary extrema. Therefore, you compute $g_i(-1) $ and $g_i(1)$. The last step is to compare the local extrema and the values at the boundary.

[pine jetty]

malakai

[rancid smelt]

yeah, that should be correct, since the inverse power method is only variant of the power method

[fair forum]

great thank you

[mystic jacinth]

is there a way to estimate global error from accumulated local error?

like just number of iterations * truncation error ?

and something to do with the order maybe idk

[molten knot]

yes! if let’s say you have an iterative method and $u_n$ is your current solution, $y_n$ is the true value. Also define
$u_{n+1}^=y_n+F(y_n)$, where $F$ say is the method and $u_{n+1}^$ is the approximation using the true solution ($y_n$), then you have
$$|y_{n+1}-u_{n+1}|\leq |y_{n+1}-u_{n+1}^|+|u_{n+1}^-u_{n+1}|$$

[pine jetty]

ansh

[molten knot]

the first term on the right is the error incurred in a single step while the second term is the error over all old steps, i.e, starting your method n+1 steps ago vs starting your method at the previous step

[mystic jacinth]

is there a way to create a bound on what the global error may be if you dont have access to the true solution?

[molten knot]

usually both of those quantities are easily boundable, like for euler’s method the first term becomes $h\tau(h)$ where $\tau$ measures the bound of the second derivative

[pine jetty]

ansh

[molten knot]

and with a lipschitz condition you can bound the second term also with h

and then you just end up getting some recurrence relation in terms of the one step errors which you can solve for the global error by repeatedly plugging into itself

[molten knot]

you could apply your method on a really refined mesh/small step length

and treat that as a “true” solution

Also for analytical purposes you don’t necessarily need the true solution

like here you just denote it as something and know what it satusfies

[mystic jacinth]

yeah I think I did read in papers that some people just use a higher order method with smaller timesteps to create a 'true' one

[mystic jacinth]

@molten knot

hey man not sure if youre around

I have no idea if this is representative of anything

but inside the integrator I take note of the timestep and the truncation error

this plots the cumulative sum for each

so if I assume it adds linearily , could I just extract this relation for my other plots to estimate a global errors on the final values?

it shows error vs seconds simulation rn

after a year simulated

[molten knot]

i think so, it’ll depend a lot on consistency and stability of your method

what sort of problem are you working on?

[mystic jacinth]

its numerical integration surrounding the n body problem

[molten knot]

also stiffness could always be an issue, some problems might be ill conditioned/stiff

[molten knot]

quadrature?

[mystic jacinth]

the orbits im considering right now are fairly circular and I dont think they are considered stiff

[molten knot]

or more like integrator in the sense of a diff eq solver

[mystic jacinth]

diff eq solver

its the dormand prince

rk(4)57m

[molten knot]

i see

[mystic jacinth]

I didnt see a point in going for a symplectic integrator for my 'short' timeframe things

[molten knot]

hmm id may compute the global error on the current step length

and a step length half the size

and measure their ratio

or the log of their ratio

and keep halving it and doing that repeatedly

see if you can get an order of convergence from that

maybe try changing parameters or the problem up and see if the order stays the same

if the method is consistent it should

[mystic jacinth]

I added a shitton more bodies and its still showing linear

I could try the log thing youre mentioning

I think it'd kinda be a pain in the ass to setup

[molten knot]

nah it’s not bad

just do like
Err = []
for n in range(1,10):
h=2^n
compute global error(h)
Err.append

take ratios of everything in list and logs

plot them with respect to h on loglog

assuming you can easily compute global error

[mystic jacinth]

I cant 💀

https://rnowling.github.io/math/2016/11/19/leapfrog-global-error.html

linear scaling of global error doesnt seem completely impossible

[rancid smelt]

Is your local error at the last timestep also linear?
Actually, you should expect a better error

[mystic jacinth]

I took the norm up to element i of the local errors and found the global error to be proportional to the root of the time integrated

Not sure if that makes sense

[wide spear]

What solver and ode

[mystic jacinth]

rkdp(4)5

acting on this

[wide spear]

Where are the derivatives

[mystic jacinth]

these?

[wide spear]

Yes

Ok

I'm not super surprised tbh

Have you tried for something like y'=y

[mystic jacinth]

no I havent, the implementation is kinda finnicky and focused around a specific thing so modifying it isnt super easy

[mystic jacinth]

is there an explanation for the sqrt(t) proportionality?

[wide spear]

Well if you do the taylor series analysis then you can probably derive something

[mystic jacinth]

[mystic jacinth]

these hitches appear if moments in the simulation arent that smooth and the time step has to be reduced

[rancid smelt]

what is the min value for r_21? could be a numerical instability to the singularlity at 0

[mystic jacinth]

theres a case for 0 to make sure it doesnt diverge

https://gyazo.com/1e3f48e6c7bc7651a3c2487395f043ef

the highly eccentric stuff is responsible for the step size having to be reduced

[rancid smelt]

OK, you implemented a proper singularity removal strategy

[mystic jacinth]

yeah

honestly ill just go with the demonstration that its growing proportional to sqrt(t)

it should be sufficient

[rancid smelt]

if you have time, you could try an rk for stiff odes

and compare the errors

[mystic jacinth]

oh I have a sarafyans fehlberg implented too

thats a nice idea

dont think its for stiff scenarios but its another solver to get a different perspective

[rancid smelt]

good point

[brave crypt]

hi

[wide spear]

Hi

Do you have a numerical analysis question

[molten knot]

are there any papers on efficiently computing sobolev norms? i have an iterative scheme where i want to compute say the H^3 norm of my solution (vector evaluated at grid points) at each iteration

doing quadrature + finite difference for H^1 is fine but anything higher kinda gets tricky in terms of speed

ideally if the result extends to W^p,q would be even better

[wide spear]

Can you elaborate on the speed concern

[molten knot]

like if i’m doing 1000 iterations, id want to compute the H^k error as quickly as possible

L^2 is fast and H^1 is fastish

but with H^2 it starts slowing down a lot since i’m applying central difference and then quadrature to it

and doing that 1000x

[wide spear]

What quadrature are you doing

Midpoint rule?

[molten knot]

yeah

for H^1 at least i am

[wide spear]

Ok I have a meeting in 6 minutes but I will think about this after

[molten knot]

it would be nice if i could just have like loss(p,q) and pass in a positive integer p or q and get the corresponding sobolev norm

and quickly since it’s computed at each iteration

but even just having an extension to H^2 would be great

[molten knot]

thanks

[wide spear]

Ok this is how I would do it

In 1d

Grid points $x_0, x_1, x_2, \ldots, x_n$ with grid spacing $\Delta x$ and function values $f_0, f_1, f_2, \ldots, f_n$

[pine jetty]

Angetenar

[wide spear]

Then the $L^p$ norms of $\mathbf{f}$ are [\norm{\mathbf{f}}{L^p}=\left(\frac12f_0^p\Delta x+\frac12f_n^p\Delta x+\sum{i=1}^{n-1}f_i^p\Delta x\right)^{1/p}]

[pine jetty]

Angetenar

[wide spear]

Ok now for higher derivatives

First, we have $f_i'=1/(2\Delta x)(f_{i+1}-f_{i-1})$ and $f_0'=1/(2\Delta x)(-3f_0+4f_1-f_2)$ and $f_n'=1/(2\Delta x)(f_{n-2}-4f_{n-1}+3f_n)$

[pine jetty]

Angetenar

[wide spear]

So [\norm{\mathbf{f}}{W^{p,1}}=\norm{\mathbf{f}}{L^p}+\left(\frac12(f_0')^p\Delta x+\frac12(f_n')^p\Delta x+\sum_{i=1}^{n-1}(f_i')^p\Delta x\right)^{1/p}] which can be simplified to \begin{multline*}\norm{\mathbf{f}}{W^{p,1}}=\norm{\mathbf{f}}{L^p}+\\Bigg(\frac{1}{2^{p+1}\Delta x^{p-1}}((-3f_0+4f_1-f_2)^p+(f_{n-2}-4f_{n-1}+3f_n)^p)+\\frac{1}{2^{p}\Delta x^{p-1}}\sum_{i=1}^n(f_{i+1}-f_{i-1})^p\Bigg)^{1/p}\end{multline*}

[pine jetty]

Angetenar

[wide spear]

This makes sense right

Next, for second order

[\norm{\mathbf{f}}{W^{p,2}}=\norm{\mathbf{f}}{W^{p,1}}+\left(\frac{\Delta x}{2}((f_0'')^p+(f_n'')^p)+\sum_{i=1}^{n-1}(f_i'')^p\Delta x\right)^{1/p}]

[pine jetty]

Angetenar

[wide spear]

And so on

This is how I would code it in C++

double loss(const int p, const int q, const std::vector<double>& fvals, const double deltax) {
    int size = fvals.size();
    std::vector<double>& derivs (size, 0);
    std::vector<double>& temp (size, 0);
    double val;
    double loss = 0;
    for (int j = 1; j < size-1; j++) { // points 1 through n-1
        loss += pow(fvals[j],p)*deltax;
    }
    // points 0, n
    loss += pow(fvals[0],p)*deltax*0.5;
    loss += pow(fvals[size-1],p)*deltax*0.5;
    
    // derivatives
    derivs = fvals;
    for (int i = 1; i <= q; i++) {
        // derivatives of order i
        temp[0] = 1.0/(2.0*deltax)*(-3*derivs[0]+4*derivs[1]-derivs[2]);
        temp[size-1] = 1.0/(2.0*deltax)*(3*derivs[size-3]-4*derivs[size-    2]+derivs[size-1]);
        for (int j = 1; j < size-1; j++) {
            temp[j] = 1.0/(2.0*deltax)*(derivs[j+1]-derivs[j-1]);
        }
        // accumulate to sum
        loss += pow(temp[0],p)*deltax*0.5;
        loss += pow(temp[size-1],p)*deltax*0.5;
        for (int j = 1; j < size-1; j++) {
            loss += pow(temp[j], p)*deltax;
        }
        derivs = temp;
    }
    return pow(loss, 1.0/p);
}

This is how I would implement it in c++

[molten knot]

thanks i’ll try to implement it tmrw and see how it holds

i just wonder how fast doing that central difference calculation is

especially for like d=2 or d=3

[river thicket]

could someone please help me solving

A is complex, but not \beta or Y

and i think that i've seen that you take the derivative of $\hat{\beta}^T$ sometimes

[pine jetty]

Jonathan

[wide spear]

What are you trying to solve for

[wide spear]

Huh

You want to find the min of L(beta hat)?

[river thicket]

yes indeed

[river thicket]

but i'd like to solve this one

[wide spear]

Are you struggling with taking the derivative?

[river thicket]

yes i am

[wide spear]

https://www.matrixcalculus.org

[river thicket]

thank you!

[wanton dawn]

use differentials to get the gradient vector with respect to beta

[neat trench]

Heya! Anyone is interested studying together Statistics and Programming languages??

Do message me soon

[wide spear]

Not channel appropriate

[royal crag]

how do i solve something like this?

[wide spear]

Have you done anything similar in class

[wanton dawn]

Expand $\frac{F(x) - F(0)}{F(0)}$ as a Taylor series around $0$.

[pine jetty]

L

[wide spear]

Maybe

This is from a numerical analysis text book

So I think that there may be some additional required consideration of floating point error

[wanton dawn]

Maybe that can be handled assuming we can evaluate $\hat{F}(x)$ where $\hat{F}(x) = F(r(x))$, with $|\frac{r(x) - x}{x}| \leq \varepsilon$.

[pine jetty]

L

[wanton dawn]

not sure what the typical assumption is for numerical anal

[royal crag]

so i have to replace F(x) with its taylor series expansion and then show that what i get is bounded?

[rancid smelt]

Lipschitz continuity is enough

[wide spear]

[rancid smelt]

[wanton dawn]

The question was kind of vague. I suggest Taylor series because it gives full information for $(F(x) - F(0))/F(0)$. You may use the expansion up to second order with the lagrange remainder to get an explicit bound for the relative error in terms of $\delta$. It's gonna be on the order $O(\delta^2)$.

[pine jetty]

L

[vale mango]

How are trancendential functions computed to arbitrary precision

I know computers usually use minimax polynomials

for float computations

but how would i compute sin(x), e^x, etc to arbitrary precision(there are a few families of generalizatons of these functions that do not have nice taylor series, but that I want to implement fast arbitrary precision algorithms for)

[wide spear]

It'll depend on the specific function

NIST has a lot of information on this: https://dlmf.nist.gov/

[vale mango]

What'd be the most computationally efficient way to solve a system of 2 nonlinear differential equations with the following very nice properties to arbitrary precision? As in, the error having a maximum error less then a given arbitrary value?

There is a conserved quantity(though it is not energy)
We know the solution is periodic, and we know the exact period.
We have a closed form solution in terms of an inverse integral(but computing the solution this way is numerically unstable and only works over part of the solution, hence why normal numerical methods for odes are prefered)
Our end goal is a fourier series for the function computed to again, arbitrary precision(I should be able to specify an error bound and it should compute fourier coefficients to the solution within that error bound).
A dense output would be preferable.
The solution over the real numbers is bounded between -1 and +1

Is there a way to directly find the fourier coefficients here?

actually, are there any books on numerical ode simulation in general?

Going into theoretical physics probability i have a hunch that will be useful

[wide spear]

Ummmm

If you could share the actual problem that would be helpful

[vale mango]

y' = x^3, x' = -y^3, x(0) = 1, y(0) = 0. Note y^4+x^4 = 1 is the conserved quantity

[wide spear]

Is RK time integration not satisfactory?

[vale mango]

I have a feeling a better method exists

I want the solution to be accurate up to machine precision for all times t

i'm trying to make a high fidelity table for some functions decribed in terms of solutions to these des

[wide spear]

If you take a fourier transform, you should end up with a nonlinear system of equations that you can solve for all the fourier coefficients

[vale mango]

the fourier transform of the equations themselves or of the numerical solution to those equations?

[wide spear]

The fourier transform of the equations

[wide spear]

Iserles is a classical text

[vale mango]

alright, this turns into -iwY(w) = F(x^3), iwX(w) = F(y^3). Then what? Sorry if this is trivial

F(f(x)) is the fourier transform of f, same as X(w) is the fourier transform of x

god i wish there was a chain rule for fourier transforms haha

[wide spear]

I'm busy now, I'll respond in an hour

[wide spear]

Ok so the idea is

The functions are periodic so instead of a full fourier transform, a fourier series is sufficient

Assume that $y=\sum_ny_ne^{int}$ and $x=\sum_nx_ne^{int}$

[pine jetty]

Angetenar

[wide spear]

Then $y'=\sum_niny_ne^{int}$ and $x'=\sum_ninx_ne^{int}$

[pine jetty]

Angetenar

[wide spear]

And $y'=x^3$ so $\sum_niny_ne^{int}=\left(\sum_nx_ne^{int}\right)^3$

[pine jetty]

Angetenar

[wide spear]

Now for the right hand side, we have that [\left(\sum x_ne^{int}\right)^3=\sum_n\sum_a\sum_b\sum_c\delta_{a+b+c=n}x_ax_bx_ce^{int}]

[pine jetty]

Angetenar

[wide spear]

Think about why this is true

delta here is the kronecker delta, not dirac

Splitting about the modes, we then have that [iny_n=\sum_a\sum_b\sum_c\delta_{a+b+c=n}x_ax_bx_c]

[pine jetty]

Angetenar

[wide spear]

And likewise, for $x_n$, [inx_n=-\sum_a\sum_b\sum_c\delta_{a+b+c=n}y_ay_by_c]

[pine jetty]

Angetenar

[vale mango]

this looks like a convolution. Can we use the FFT to speed this up at all?

[wide spear]

What are you going to take the FFT of

Anyways

To solve this numerically, truncate the range for n, a, b, and c

a, b, c, n in [-N, N]

x and y should be quite regular so y_n and x_n should decay quickly away from 0

[vale mango]

could you give an example computation here for say x_1? I'm having a bit of trouble wrapping my head around this

again, thanks for this though

[wide spear]

This is a nonlinear system of equations

[vale mango]

yeah, could you show how to write out that system explicitly i guess is what i'm asking

[wide spear]

Let $X=[y_{-N},y_{-N+1},\ldots,y_{N-1},y_N,x_{-N},x_{-N+1},\ldots,x_{N-1},x_N]$

[pine jetty]

Angetenar

[wide spear]

Define $F(X)=inw_i-\sum_{a=-N}^N\sum_{b=-N}^N\sum_{c=-N}^N\delta_{a+b+c=i}z_az_bz_c$

[pine jetty]

Angetenar

[wide spear]

Apply newton's method

The notation is kind of scuffed

But it shouldn't be too hard to figure out

Anyways I wouldn't recommend doing this

Directly doing RK-n will be much easier

[vale mango]

This is probably a stupid question, but truncation error won't matter with enough steps right?

As in there is some nonzero stepsize i can employ such that the total error over the entire interval will be less then any given error bound right?

[wide spear]

This is a bit of a tricky question

There is a small enough time step such that the total accumulated truncation error will be < epsilon

However, if you're doing this numerically, you should be careful of floating point error as well

Which is on the order of 10^-16 for double precision and will acumulate linearly with the number of time steps

If you need smaller errors, then quad precision has floating point error around 10^-34

[vale mango]

now that i think about it because i'm only going to be doing this computation once to make a table and to find fourier coefficients, i could just use multiprecision floats

yes it's more expensive but it's not like i need to do this over and over again

[wide spear]

When you actually do the fourier coefficient part, you should be careful about the period of the functions

[vale mango]

My goal is to just find a fourier series that is accurrate to within 10^-34 for all inputs(so that this series is equivilent to the actual function for quad floats)

[vale mango]

the period is known 100% exactly. It's a messy thingy involving the gamma function but it's an exact value

[wide spear]

Have you seen this

[vale mango]

glancing over this, yes there is a closed form solution for the integral of sin_4(x)

It was both listed in a textbook i found on the subject and I found it independently

the integral of sin4(x) iis 1/2*sin4(x)^2 * 2F1(1/2, 3/4; 3/2; sin_4(x)^4) +C

[vale mango]

where 2F1 is the hypergeometric function

[vale mango]

@wide spear Actually, I know the function exactly at a few other points. Does this make it any easier?

[wide spear]

It does not, no

[vale mango]

1. There is actually a simpler form of the equation, x' = (1-x^4)^3/4, x(0) = 0. You can do something similar here right? Just raise both sides to the fourth power basically.
2. This also works for just the sine part of the transform, sicne we know this function is odd, right?
   3, doing this for the sine part, we get.

b_1^4*w^4*cos(wx)^4 = (1-(b_1sin(wx))^3)^4

we can expand the right, but then what?

do we plug in x=0?

also for n>1 wouldn't that system not be uniquely determined

[vale mango]

wait nvm

i get this

I think i might have found a better way

I reduced it to solving a weird dft equation for a sequence basically

[vale mango]

how do you just get rid of e^int here?

[wide spear]

For this step?

[vale mango]

yeah

[wide spear]

You can separate each fourier mode

If $\sum a_ne^{int}=\sum b_ne^{int}$ then $a_n=b_n$

[pine jetty]

Angetenar

[vale mango]

alright

[vale mango]

also, i think i have a faster way to do this

,tex
$(\sum^{k}_{n=0} a_n)^{p}$ is, essentially, taking a convolution of the sequence $a_n$ with itself p times. So use FFT fast convolutions instead of a triple sum.

[pine jetty]

Colonizor48

[wide spear]

How are you going to take the convolution of a sequence that you don't know the values of

[vale mango]

,tex
Rewrite $(\sum x_ne^{iwnt})^3$ = $(\sum \mathscr{F}^{-1}(\mathscr{F}(x_ne^{iwnt}})^3)$. where ${x_n}^3 = {x_n^3}$

[pine jetty]

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

[vale mango]

Do the same to the other side

,tex
So we get $iw\sum\mathscr{F}^{-1}({\mathscr{F}(n*y_ne^{iwnt})$) = $(\sum \mathscr{F}^{-1}(\mathscr{F}(x_ne^{iwnt}})^3)$

where here, the sum operator simply sums over a finite sequence

note how the composition of the sum operator and F^-1 operator is a linear operator

therefore, we can factor the sum F^-1 operator out of both sides and get

$iw\mathscr{F}(ny_ne^{iwnt})$ = $\mathscr{F}(x_ne^{iwnt}})^3$

[pine jetty]

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

[vale mango]

now this is reduced so solving for the sequences x_n and y_n

which can probably be done using newtons method, but could maybe also be done symbolically or through some other method

[pine jetty]

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

[vale mango]

Also I think you can cancel out e^iwnt here right

[wide spear]

How

[vale mango]

Note that the fourier transform of the sequence can also be interpreted as a multivariate vector function

so we can start with guesses for x_n and y_n and use netwons method or something

or there could be a way to solve analytically\

[vale mango]

Do this first, cancel the modes, then undo the step and rewrite the sum as a sum cubed again without the exponential term, then apply the fft optimization

start with an initial guess for the series of fourier coefficients, then use some sort of method to solve

@wide spear Does this hold up?

[wide spear]

I'm not sure what this means

[vale mango]

can someone direct me to resources on good general linear methods?

[fair forum]

is order of convergence the same as order of accuracy

[wide spear]

They're used in different contexts

Order of convergence is used for iterative methods

[fair forum]

ok thanks

[wide spear]

Write your own c++ code

[wary oriole]

    """
    Rayleigh Quotient iteration.
    Examples
    --------
    Solve eigenvalues and corresponding eigenvectors for matrix

             [3  1]
        a =  [1  3]

    with starting vector

             [0]
        x0 = [1]

    A simple application of inverse iteration problem is:

    >>> a = np.array([[3., 1.],
    ...               [1., 3.]])
    >>> x0 = np.array([0., 1.])
    >>> v, w = rayleigh_quotient_iteration(a, num_iterations=9, x0=x0, lu_decomposition="lu")

    """
    x = np.random.rand(a.shape[1]) if x0 is None else x0

    for k in range(num_iterations):
        sigma = np.dot(x, np.dot(a, x)) / np.dot(x, x)  # compute shift
        x = np.linalg.solve(a - sigma * np.eye(a.shape[0]), x)
        norm = np.linalg.norm(x, ord=np.inf)
        x /= norm  # normalize
        if verbose:
            print(k + 1, x, norm, sigma)
    return x, 1 / sigma

hi all, im trying to implement rayleigh_quotient_iteration here.
but I don't get this graph of calculation by my own hand calculation tho

so I set x0 = [0, 1],
a = np.array([[3., 1.],
... [1., 3.]])

then I do hand calculation, first sigma is indeed 3.000, but after solving x, the next vector, I got [1., 0.]
how the hell the book got [0.333, 1.0]?

where is this k=1 line from?

I did hand calculation, after first step x_k is wrong. x_1 = [1., 0.] after normalization it's still [1., 0.]

[robust oar]

c

[wide spear]

[royal crag]

what should i get after first iteration?
i got (8/3, 4/3) but im not sure if its correct

[vale mango]

how do i derive new general linear symplectic integrators

[wide spear]

Are you familiar with how existing symplectic integrators are derived

[vale mango]

No learning about that would be useful

[vale mango]

@wide spear

[fathom rain]

[ 1.95, 1.75] (it asked for 2 iterations)

[wary oriole]

I found this notes in the Trefethen book. seems industy standard like matlab and LAPACK has better
Stopping Criteria than regular things we write ourselves. Does anyone know what they usually uses?
Is there some paper on stopping criteria?

[rancid smelt]

normally you have a while loop with the stopping criteria $|| v^{(k)} - v^{(k-1)}||_2 > \varepsilon$ with $\varepsilon > 0$

[pine jetty]

malakai

[wary oriole]

Yeah I know this. This is very blend off the text book solution. I think what this meant is there are advanced
industrial practice on how to implement stopping criteria

[rancid smelt]

sorry to disappoint you, there aren't any advanced ones at least for the power iteration

[wary oriole]

I think it meant not just for power iteration, but all iterative algorithms

[wide spear]

This is a very common stopping criteria

I'm not sure what you are looking for in terms of advanced industrial practice

A lot of it is very ad hoc

[rancid smelt]

i suppose that yehuihe believes that there are some special ones only used in advanced industrial practices which are not teached in generall

[wide spear]

If you have a specific problem, sometimes you can use some problem structure for a different stopping criteria

[wary oriole]

Yes the way this book sentence framed pretty clear right? termination conditions are one of the major differentiator for my own solution and top-quality software. Anyway

[wide spear]

You can see what stopping conditions lapack uses, it's all open source

But tbh, I suspect that the stopping conditions that lapack uses are not that fancy

I took a class with one of the lead devs

[rancid smelt]

booth Matlab and lapack use the standard ones for each algorithm.
fancy ones you only have in codes adapted for special problems

[wary oriole]

Ok got it. seems these are for generic robust algo in the library. Thank you

[wary oriole]

oh one more question. In inverse iteration, we suppose to solve A-mu*I. If we choose a shift so close to an eigenvalue then we will converge very fast. If we choose mu to be an eigenvalue then should take one iteration. This suppose to make inverse iteration very useful if we know what eigenvalue already to find out eigenvectors

nearly singular ill-conditioned shift close to an eigenvalue cause no trouble indeed as explained in exercise. But if I use exact eigenvalue then shifted A is a singular matrix, which cannot be solved. How should I modify this in my code for this case?

for example A=[[3, 1], [1, 3]], eigenvalue 4 and 2. If I choose shift to be 4.000001, it converge in two iteration to 4 and find the eigenvector. It works

but if I choose shift to be exactly 4, then I can't solve this equation, therefore failed.

[wide spear]

Well if you can't solve the equation, then you can terminate the iteration immediately

[rancid smelt]

it depends, if you are only interested in the eigenvalue, you are fine. If you want to know the eigenvector and as in your case \mu = \lambda, then you can restart your algorithm with \mu = lambda + \varepsilon in order to compute the eigenvector

[wary oriole]

Yeah the point is if shift is exactly the eigenvalue then you can get eigenvector in a single iteration. if I add a epsilon then it will take two. But I get you. I think there is no nicer way to get it except adding an epsilon

[hallow flint]

hey, can anyone pls recommend me a good yt playlist for learning the mentioned numerical methods syllabus?

[wide spear]

Consider reading a book

[hallow flint]

can you pls suggest me some good numerical methods books?

[wide spear]

I think this roughly corresponds to burden and faires

[hallow flint]

thank you

[crystal badge]

hi guys

ive this course for my 4th sem
and i don't understand it

can someone suggest me tips to study

[sour needle]

whats about exactly?

[crystal badge]

i tried studying from the suggested readings, i just don't get how to apply the methods

like lu decomposition, guass jacobi seidel etc

[sour needle]

oh i see

do you have to do proofs or programming or practical exercises?

[crystal badge]

all, but rn the test is of exercises

where you gotta do iterations and stuff

[sour needle]

in paper?

i suggest you do the basics and know very well when to apply it or not

if the matrix is positive-definite etc

[crystal badge]

yeah from where do I do the basics

[sour needle]

i would follow a book

but my advice is so trivial im sorry

i do think thats the basis for finite differnces

[wide spear]

What is a specific question that you are struggling with

[wide spear]

far? aside?

Well being conservative is pretty useful

What qualifies as remarkable to you?

Well conservation is pretty numerically useful

That's the primary benefit

It's also easy to formulate for unstructured meshes

From a math perspective, it's also interesting

There is the notion that finite volumes correspond to an alternative weak formulation

Finite element corresponds to the standard weak formulation right

Finite volume methods correspond to an alternative weak formulation where you formulate the pde in terms of its average value

[tired depot]

how the proof for local convergence work

idk how this starts like this

quite werid

[brave crypt]

Its the thing where the solution isnt smooth or whatever

Actually idk

I just know its the thing where you multiply by test function and integrate

And use integration by parts

Ok im wrong sorta

[wide spear]

Ok so this is the most common notion of a weak solution

But there are others

Like the viscosity solution

I remember reading something about averaging as another form of weak solution but I don't remember where and I've been trying unsuccessfully to find a reference

[brave crypt]

Yeah no worries

[pine jetty]

Mandy Moose

[wide spear]

The variational formulation corresponds to finite elements

Distribution weak solution and test function weak solution are the same

[brave crypt]

My confusion is how do you write solvers for coupled diffeq in general

Ok bad question

Ig my problem is confusion with update steps and stuff

[wide spear]

Start with y''=-y

[brave crypt]

Yeah but that isnt what i mean ig

Sorry im on phone so i cant actually ask question that well

But more like I only know how to write solvers for linear higher order odes

[wide spear]

If you write out a model problem I'll take a look tomorrow

[brave crypt]

Thanks ill do exactly that

Ill just make some shit up

And show my attempt

Ig i have confusion with that + nonlinear solving but I have to learn more about nonlinear myself me thinks

Because ive only done nonlinear solving for second order diffeq

Nah im lying lemme stop being disrespectful and give you full details when im ready

[wary oriole]

I have a conceptial question. In the Trefethen book when he mentioned he only restricting to symmetric matrices for power iteration, inverse iteration, and rayleigh iteratioin. But all computed rayleigh quotient in the algo. Does all three works for Non-hermitian matrix as well?

I asked chatgpt it says not for rayleigh quotient iteration but works for inverse iteration. Which is strange cause in inverse iteration also calculate rayleigh quotient

[rancid smelt]

@wary oriole i suppose the issue is that the eigenvalues of real symmetric matrices are real, too. Therefore, you don't have to consider the complex case.

[wary oriole]

These methods works in general for hermitian matrix. Not just real symmetric ones

My question is whether it works for non hermitian

[rancid smelt]

for power iteration and inverse iteration, obviously

and rayleigh, too.

[rancid smelt]

and if you want you can use M = A^H A instead of A
for the sign of the eigenvalue you can use one step of the power itertion used on A

[hollow belfry]

In some situations they do and in some they dont. It depends on the matrix.
However, they lose most of the nice properties that are guaranteed for hermitian matrices.

For the non-hermitian case there are more specialized methods like Arnoldi Iteration and LOBPCG

[wary oriole]

Yeah seems like it. like cubic convergence.
but does power iteration and inverse iteration the same? I thought inverse iteration works in general sense

[wide spear]

Power iteration should work for general matrices

You can prove the convergence using Jordan form

[hollow belfry]

Power iteration and inverse iteration are the same exact method. Inverse iteration is just power iteration applied to an inverse matrix in order to ensure that the smallest eigenvalue becomes the largest

[brave crypt]

Hey

[mellow island]

Hello.

[wide spear]

Hello

[mellow island]

I kinda misunderstood Hermit method + Tchevechyv.

Hang on,

[wide spear]

Chebyshev?

Persumably Hermite also

[mellow island]

I always misname it but thanks.

Do you understand french or I need to use chatgpt?

[mellow island]

Hello?

[wide spear]

Ummmm

I can try

[mellow island]

Page 82 - 84 Hermite methode.

[wide spear]

Ok hermite interpolation

[mellow island]

So, since I was in Newton and lagrange methode, it was direct.

The main formula and error formula.

Since I started diving into Chebyshev and Hermit, I couldn't keep up with the author...

[wide spear]

Do you understand the idea behind hermite interpolation

[mellow island]

More accurate interpolation?

Than the rest of methodes since it uses derivatives?

[wide spear]

Well so with newton/lagrange interpolation, you construct the interpolant purely using function values

With hermite, you also use derivative information

So it should be more accurate

[mellow island]

Exactly.

[wide spear]

Yes

[mellow island]

But.

[wide spear]

Yes?

You don't need to translate them

I can read this french

[mellow island]

Oh alright.

So, it confused me with that example...

Why do I need to know it is linear independent?

[wide spear]

So the H_0, H_1, K_0, and K_1 need to be linearly independent

This is sort of a fundamental property when doing interpolation

[mellow island]

Those are the base of hermite, right?

[wide spear]

You consider a basis for a function space

Yes, the hermite basis functions

[mellow island]

Alright, what is H and K?

Nvm

I mean this.

K is associated with derivatives.

[wide spear]

$K_0$ is constructed to have the properties that $K_0(x_0)=K_0(x_1)=K_0'(x_1)=0, K_0'(x_0)=1$

[pine jetty]

Angetenar

[mellow island]

This what I don't understand.

[wide spear]

$K_1$ is constructed so that $K_1(x_0)=K_1(x_1)=K_1'(x_0)=0, K_1'(x_1)=1$

[pine jetty]

Angetenar

[mellow island]

When it equals to 1 or 0?

[wide spear]

H_n equals 1 at x_n

K_n' equals 1 at x_n

[mellow island]

What?

[wide spear]

Yes

[mellow island]

Do you mind explaining it to me please?

[wide spear]

I'm not sure what you don't understand

[mellow island]

I'm a bit slow in these stuff (even in lagrange nominator)

So,

k indicates what exactly?

if i indicates which point..?

[wide spear]

i and k are indices

[mellow island]

Yes, i follows the point, k follows what?

[wide spear]

n is the number of interpolation points

k indexes the basis polynomials

[mellow island]

n is the max number.

basis polynomials?

[wide spear]

H_0, H_1, etc..., K_0, K_1, etc....

[mellow island]

So, how on earth i doesn't equal to k?

[wide spear]

H_1(x_0)=0, H_1(x_1)=1

[mellow island]

alr alr alr.

But how is that here?

[wide spear]

What do you mean

[mellow island]

Do you agree that the interpolation of hermit give that polynome?

[wide spear]

Yes, that is the hermite interpolant for two points

[mellow island]

Exactly.

I can literally see the indice k in that polynome, but how does x change?

If it doesn't bother you, can we join a call, please?

[wide spear]

x is between x_0 and x_1

Sorry I can't call right now

[mellow island]

It is your right, mate.

Anyway,

As I understood,

We have 2 points for example.

p(2n+1)= y0H0 + z0K0 + y1H1 + z1K1, right?

[wide spear]

Yes

[mellow island]

alright.

Is it the following?

y0(H0 + H1) + z0(K0 + K1)...ect?

[wide spear]

No

y_0 and y_1 are different

[mellow island]

I still cant reach/find the case where I need use this:

[wide spear]

You use those to check that the p(x) defined on page 83 satisfies 4.12 and 4.13

[mellow island]

OHHHHH

I see I see.

So when you get the polynome (final result)

Instead of calculating H and K (exhausting) you just replace with their values directly, yeh?

[wide spear]

At x_0 and x_1 you can just replace with the values

[mellow island]

like p(x0) = y0(1) + y1(0) + z0(1) + z1(0) which means p(x0)= y0 + z0 right?

[wide spear]

$p(x_0)=p_0H_0(x_0)+p_1H_1(x_0)+p_0'K_0(x_0)+p_1'K_1(x_0)=p_0$

[pine jetty]

Angetenar

[mellow island]

p0 + p0'

[wide spear]

K_0(x_0)=0

[mellow island]

Since in K both i and k are identical in the case of i=k=0.

Oh yes, thanks for reminding me.

Silly question but

I can use either H to verify OR K to verify, right? (It all went down hill after BAC)

[wide spear]

?

You need to use both

[mellow island]

Bro..

What is H'?

[wide spear]

H"(x_0)=H"(x_1)=0

[mellow island]

They don't even exist in the formula...

[wide spear]

You can explicitly compute H_0' and H_1' from the given formulas on page 82

[mellow island]

Bro... why do I need H' and K'?

I hope you bear with me because These 2 methods are the only pain in the ass for me.

[wide spear]

You need K' to capture the derivative information correctly

[mellow island]

But I still can't find them in the polynome formula...

[wide spear]

This formula?

[mellow island]

Yes.

[wide spear]

They are not used in the definition of p

But if you check that $p'(x_0)=p_0'$ and $p'(x_1)=p_1'$ then you need to compute $H'$ and $K'$

[pine jetty]

Angetenar

[mellow island]

hang on.

I GOT IT.

So, finding the polynome is the easier part.

But verifying it is the hardest part, we need to have the full equation of H and K and H' and K', this is why we have those propities to avoid all that wor, right?

[mellow island]

Am I right?

[wide spear]

I guess so

[mellow island]

Thanks man.

You got tired or shall I continue?

[wide spear]

You can continue

[mellow island]

I'm in page 84, cas generale.

[wide spear]

Ok

[mellow island]

Alright.

So, why m is important?

[mellow island]

What it is highlighted?

[wide spear]

m is the number of derivatives that you use