numerical-analysis

[mellow island]

Number of derivatives or how much we deriv?

[wide spear]

?

m=1, you use information about f and f'

m=2, you use information about f, f', and f''

[mellow island]

So,

Besides the 1st method.

We use this to determine H, right?

H isn't the base??????

Bro what?

[wide spear]

H is the interpolating polynomial now

[mellow island]

I see.

Why do I need this now?

Hang on.

So, I use this process instead of... just using lagrange methode

[wide spear]

Well

You can

It's an option

[mellow island]

Depends on the exo, got it.

I think I understood.

The last part.

[wide spear]

Ok

[mellow island]

Btw.

The more points I have , the more m grows larger right?

[wide spear]

More points = larger n

[mellow island]

N = the sum of m + 1.

But stillm how to identify the value of m?

[wide spear]

You choose m

[mellow island]

WHAT?

Example please?

[wide spear]

Well on pages 82 and 83, we were working with m=1

But you can choose larger m if you want

[mellow island]

Alright, here, let suppose m=10, what is m0,m1...ect?

Do they have to have a equidistinct values like 1,2,3...ect or doesn't matter?

[wide spear]

Well

You can choose m non uniformly

So m0 is the number of derivatives that you know at x0

m1 is the number of derivatives that you know at x1

[mellow island]

Alright.

I think this it.

Thanks a lot man.

You helped me a lot.

[wet vector]

Hey! I've been trying to find integrate the analytical solution to this particular Burgers problem, but I couldn't find a way to actually solve it so I tried numerically integrating the problem with scipy integrate. For some reason, I can't get the exact solution which was found in the paper, does anyone happen to know what I might've missed?

The last one I clipped some points

[wide spear]

All 3 plots are titled exact solution

[vapid plume]

<@&268886789983436800>

Hello. Can anyone recommend a textbook for C/C++ in the context of numerical methods?

[wide spear]

I would never recommend a textbook to learn a programming language

[molten knot]

i agree with angetenar but this might be kinda related to what you're looking for

https://yangpl.wordpress.com/wp-content/uploads/2017/12/parallel-scientific-computing-in-c-and-mpi-2003.pdf

its a good book that happens to use c++

[vapid plume]

Thank you, I'll try it out.

[wary oriole]

numerical recipe is very standard

[wide spear]

?

[neat grail]

..?

[brave crypt]

?

[tall solar]

Does anybody here know how chebyshev systems occur in numerical analysis?

[wide spear]

https://encyclopediaofmath.org/wiki/Chebyshev_system

This chebyshev system?

[tall solar]

Yeah exactly do you know about these? For example, I’m curious how relevant these are to numerical analysts and how they use them. I have an idea that these chebyshev systems allow for a unique best approximation. Do people spend time constructing these for particular problems? I don’t know if this is a niche topic in numerical analysis.

[fathom rain]

never heard of it (and probably rather niche, 1 in arxiv in 25, 0 in 24, and 3 in 23 using standard spelling)

[tall solar]

Thanks I didn’t think of checking that

[wide spear]

Not relevant

Optimal polynomial interpolation in R^n is solved

Optimal interpolation on uniquely shaped domains is impossible

[tall solar]

How do chebyshev systems relate to optimal polynomial interpolation in Rn? I am pretty sure they have a unique closest point propery but I don’t see how that’s the same as polynomial interpolation. I could have a chebyshev system with functions that are not polynomials. I probably don’t understand how they aren’t relevant.

[wide spear]

I don’t know how chebyshev systems relate

I can tell you how to do optimal polynomial interpolation though

[tall solar]

That’s fair okay I see. Optimal polynomial interpolation is solved. That is helpful thanks

[wide spear]

On some domains

[tall solar]

And also, maybe optimal polynomial interpolation is only solved for some system of polynomials? I don’t know optimal polynomial interpolation but I imagine there are chebyshev polynomials P1,…,Pn and for a given function f you since a best approximate c1P1+…+c1Pn. But maybe the result you had about optimal polynomial interpolation is only solved for that choice of basis Pj. Unless optimal polynomial interpolation doesn’t just mean interpolation by chebyshev polynomials but instead a more broad family of basis functions

[long lake]

Wait, so how do you do it? And optimal in what sense?

[wide spear]

It doesn't matter what the basis is, they're all different ways of writing the same polynomial

[wide spear]

Use function values at chebyshev points

Optimal in the sense of minimal lebesgue constant

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

[long lake]

Oh interesting! I'll check this out thanks

Out of curiosity, did you learn about this in a class or on your own at some point?

I should probably take a class on numerics

[wide spear]

Classes on numerical analysis will probably not touch on the issue of optimality

[long lake]

Ahh

[tall solar]

I’m confused about how you use a basis here. I thought if I had two different basis p1,…,pn and q1,…,qn then for a given function f there will be two different best approximations a1 p1+…+an pn and b1 q1+…+bn qn. My point was that just because it’s hard to compute a1…an for f doesn’t mean it’s hard to compute b1,…,bn for f. Maybe for any given basis there will exist a function it is hard to approximate. Then it would make sense to me that this sort of approximation is irrelevant.

[wide spear]

Well

Translating between different polynomial bases is easy right

[tall solar]

Well suppose I had a basis consisting of sin(t),cos(t),sin(2t),cos(2t),…sin(10t)cos(10t) for t in [0,1]. I don’t think linear combinations of these functions gives the same space as the space generated by linear combinations of chebyshev polynomials. Do you know if it’s hard to approximate with this basis of trig functions? Maybe the result is different than for “polynomial interpolation”

[wide spear]

With trig functions you use a FFT to get the interpolation coefficients

[tall solar]

Okay so for some basis it’s possible. This is a chebyshev system

[wary oriole]

How to solve a 3d laplace equation with finite difference? This is a from mary bioas's mathematical methods book, a standard BVP in cylindrical coordinate

analytical solution are well known. But I'm wondering how to solve it numerically

[wary oriole]

so far I only done plain old 2d laplace with FD on rectangular domain. I know how to discretize that. But never done in curvilinear coordinates or in 3d

[wide spear]

What is the laplacian in cylindrical coordinates

[wary oriole]

$\Delta u=\frac{1}{r}\pdv{r}\left(r\pdv{u}{r}\right)+\frac{1}{r^2}\pdv[2]{u}{\theta}+\pdv[2]{u}{z}$

[pine jetty]

yehuihe

[wary oriole]

this

[wide spear]

Ok

Can you discretize this

[wary oriole]

no. never done that. i only done in rectangular coordinates

[wide spear]

They're just derivatives

You can do it

[wary oriole]

like this? assuming all $\Delta x=h$ are the same

[pine jetty]

yehuihe

[wide spear]

Please put it in latex

You'll probably have a Delta r, Delta theta, and Delta z

[wary oriole]

$\frac{1}{r}\frac{u_{i+1,j,k}-u_{i-1,j,k}}{2\Delta r}+\frac{u_{i+1,j,k}-2u_{i,j,k}+u_{i-1,j,k}}{(\Delta r)^2}+\frac{1}{r^2}\frac{u_{i,j+1,k}-2u_{i,j,k}+u_{i,j-1,k}}{(\Delta \theta)^2}+\frac{u_{i,j,k+1}-2u_{i,j,k}+u_{i,j,k-1}}{(\Delta z)^2}=0$

[wide spear]

Not sure why you have u_{i,j,k} in denominators

1/r is not 1/u

[pine jetty]

yehuihe

[wary oriole]

this?

[wide spear]

Closer

Why is there a 1/theta^2

[wary oriole]

$\frac{1}{r}\frac{u_{i+1,j,k}-u_{i-1,j,k}}{2\Delta r}+\frac{u_{i+1,j,k}-2u_{i,j,k}+u_{i-1,j,k}}{(\Delta r)^2}+\frac{1}{r^2}\frac{u_{i,j+1,k}-2u_{i,j,k}+u_{i,j-1,k}}{(\Delta \theta)^2}+\frac{u_{i,j,k+1}-2u_{i,j,k}+u_{i,j,k-1}}{(\Delta z)^2}=0$

[pine jetty]

yehuihe

[wary oriole]

this. a lot of mistakes haha

[wide spear]

Yes that looks better

Anyways now that you can discretize, you know how to solve right

[wary oriole]

by forming a linear system? this is not that trivial tho. Laplace equation in rectanglar domain is very standard, this is horrible

[wide spear]

You can still form the linear system right

[long lake]

(Deleted an above message since it wasn't about numerical analysis)

[fluid cliff]

Okay:

i want to know if a non diagonalizable matrix is convergent

What do i do

Because nornally id check the spectral radius of the matrix but idk if that also works for non diagonalizable matrices, chatgpt says it does but i need a source

I dont need a true proof, just a reliable proof stating it so i can just do it with the spectral radius method

OR IF CHATGPT is wrong
Please tell me asap hahaha
And if possible help me find the correct method

Oh and if this isnt numerical analysis also pls tell me because i dont really know what this is haha sorry mb 😭

Sorry i notice this msg seems a bit rude ive been talking to chatgpt to much lately hahaha, my apologies

[wide spear]

What do you mean by convergent

[hasty ferry]

Hello

[wide spear]

Are you familiar with recurrence relations

And the techniques and methods used to "solve" them

[hasty ferry]

Not really

[wide spear]

Ok

Well

Hmmm

Can you derive T_n(cos(x))=cos(nx) from this?

Then, use the complex exponential definition of cos and expand

[hasty ferry]

Yes i showed it in this way

[wide spear]

Ok

If you want to avoid cos

[hasty ferry]

but isnt it only for x in [-1,1]?

[wide spear]

Then use induction I guess

Prove this for T_0(x), then assuming that it's true from T_n, prove it for T_{n+1}

[hasty ferry]

so from that we get this expression but only for x in [-1,1], am i right? Or is there any way of expanding it

[wide spear]

Yeah I suppose so

[hasty ferry]

Hmm actually thats almost the only way i know for proving these things 😄

[hasty ferry]

im stuck here

but one can not see it

[wide spear]

Oh my

[pine jetty]

Sentinel

[hasty ferry]

yeah it doesnt look nice

[hasty ferry]

ok i got it now 🙂

thx

[river current]

Hello, I am looking for insight on how to approach question 4. One of the first thoughts that comes to mind is whether I can simplify this expression further to hopefully make the calculation more meaningful. After some simplification, I arrived at the form

2x^2 * sqrt(9 + 7/x^2))

However, this form still doesn't seem to work well with the large x=76543217654321.

Any insight on how to correctly solve this problem would be greatly appreciated!

[wide spear]

Find a taylor series for sqrt(9x^4+7x^2) centered on 3x^2

[river current]

Is a first order taylor series expansion enough?

[wide spear]

No I don't think first order will be sufficient

I also don't think that the answer is correct

[fathom rain]

it is correct

(tested in julia)

[wide spear]

Is it

Mathematica disagrees

[fathom rain]

f (generic function with 1 method)

julia> f(BigInt(76543217654321)//1)
-1.16666666666666666666666666662794724165303930075401550573758460966712273705370935088660952474621511024478206717357759205978144501315472210802934002873285482```

[wide spear]

Hmmm probably mathematica being funky

[fathom rain]

or i have some minor bug? because if i increase accuracy it does not tend to -7/6 (maybe the answer is not that)

i do not see anything contaminating the types

it is 17576592506630274028349913123-sqrt(308936604144131499563240274478770188749226284147797410416)

and that is not -7/6

[wide spear]

Manually setting the precision fixes it

[brave crypt]

Why do you use rational int

[fathom rain]

To avoid Float64 contamination of BigFloat computations (most of iit not needed but I usually overdo it because it is so annoying when in happens)

[brave crypt]

Ok

[brave crypt]

I have to compute the comatrix of $I + D^2u_n$ where $D^2$ is the hessian matrix with an example being $(D_{12}u){ij}=(u{i+1,j-1} - \dots + u_{i-1,j-1})/(4h^2)$

My question is how is $I+D^2u_n$ a matrix? If $u_n$ is a lexicographic ordering of the grid points then I dont really understand this expression to make a matrix

[pine jetty]

bongoboy44

[brave crypt]

So far ive just written finite difference matrices for each D11,D12=D21,D22

But I cant really interpret the D2un expression as a matrix

For reference

I interpreted the expression $M_{ij}$ as $I + D_{ij}$, but clearly this is wrong because it somehow depends on $u_n$

[pine jetty]

bongoboy44

[brave crypt]

They go on to say this

But its not the finite difference scheme in used to hearing about

Is this a pointwise vs global fd scheme thing?

[fathom rain]

Don't know about the syntax but isn't that just central 2nd order fd first derivative in 2 dimensions? So, kron(D,I(n))+kron(I(n),D) where D is just -1 on one offdiagonal and 1 on the other?

[brave crypt]

The one they talk about in paper would be kron(D,D)

Or you can write how you wrote it and put a multiplication in middle instead of addition

But yeah thats the right D form

But this is for a global derivative defined on a boundary

Im curious if they mean locally computing derivative, but this would confuse the fuck outta me

But i cant really make sense of the constructing comatrix part

[fathom rain]

I downloaded the paper, I will not parse that syntax 😛 French people always have to make stuff unreadable and unnecessary complex 😄

[brave crypt]

Lol i need to tho 😢

[fathom rain]

fun fact glowinski (the editor) is my scientific grandfather 🙂

[brave crypt]

Like I dont think they want you to compute 512 x 512 grid

Comatrix for 512 x 512 matrix sounds insane

[fathom rain]

i dunno, sorry

[brave crypt]

Or ig i could get around that

Ye no worrie

Ill just try like 6 things out

[wide spear]

I’m pretty sure that the matrix is supposed to be I+D^2

[brave crypt]

Then how am I supposed to make sense of the algorithm update comatrix step?

The part where it says assembling M^n

Like I+D^2 u_n

They can’t mean thats just I + D^2

But if they mean its (I+D^2)u_n like applied then this makes less sense to me

[wide spear]

Apply the comatrix to u_n?

Is this a solution for monge ampere

[brave crypt]

Yes

[wide spear]

Surely there are better written papers

[brave crypt]

I can link paper

Yes there are

But i chose this one

And sadly I have to stick with it

[wide spear]

Unfortunate, sorry for your loss

[brave crypt]

I want to find something using similar method though

Because interpreting this makes no sense

Like just to check im not crazy

But computing the comatrix of the result of a FD approximation to D^2 would have no use

Let alone multiplying that by another matrix

Like I would understand computing to comatrix of D2 not this though

Holy shit

They could mean take M_ij to the nth power

Ill kms if this works though

[wide spear]

I hope not

[brave crypt]

No that makes especially no sense

[brave crypt]

Been trying to understand what it could be for last 2 days tbh

My conclusion is that I am more confused than I was initially

[brave crypt]

@fathom rain @wide spear I found another paper discussing the algorithm update steps https://arxiv.org/pdf/1009.6039

This definitely seems cleaner

[wide spear]

This doesn't answer the question though does it

What is det(I+D^2u_n)

[brave crypt]

It doesnt lol

Like im pretty sure its needs to be a matrix

But I dont see any possible consistent interpretation for that

[wide spear]

Ummmm

Maybe D^2u is the matrix of second derivatives

[brave crypt]

Yeah they tell us its the hessian

But it cant be where D11 represents dxx at each point on grid

Because then the dimensions for the equation wouldnt make sense I dont think

[wide spear]

Urgh

What if it's det(I+D^2u_n) at each grid point separately

[brave crypt]

I thought about this too

But what is D2 theta then?

I had assumption D2 theta could just be the normal fd matricies

Oh wait lol

[brave crypt]

This is how I calculated the determinant

I was fine with that

And using same assumption I was letting my M_ij be the comatrix of normal 2D hessian

So like $I + D^2$ having
$\begin{pmatrix}
\partial_{yy}+1 & -\partial_{xy}\
-\partial_{yx} & \partial_{xx}+1
\end{pmatrix}$

[pine jetty]

bongoboy44

[brave crypt]

And then letting the M_ij be the respective blocks of this

But when I multiplied MijDij and solved for rhs I was getting results blowing up after 2 iterations

[wide spear]

What tau did you use

[brave crypt]

1

[wide spear]

Oh

That's not damped

Try tau=10 or something

[brave crypt]

The original paper used 1

[wide spear]

Oh

[brave crypt]

Ok imma try 10

[wide spear]

Why do they call it damped then

[brave crypt]

Lol because they have graphs that vary it

Also my blow up was crazy so I had a feeling i was doing something deeply wrong

But the determinant i calculated was convering to the true determinant

But it takes like 500+ grid size to have reasonable L infty error

Meanwhile the paper is using like 16 grid size to get their convergence

[brave crypt]

Lol

Tried something dumb

It still makes no sense to me

But i pretty much multiplied each entry of the D^2 matrix row wise by grid evaluations of comatrix evaluated on the grid

Then turned damping to 10

Got the error to be under 10 and decreasing lol

[brave crypt]

Ok

In their masters thesis they wrote about it

This is the stencil they gave lol

Im not sure where it comes from

Im confused

Cuz this is a fourth order

[brave crypt]

Yea ill have to give up );

[brave crypt]

The paper lied

I used a fourth order scheme and got convergence

[wide spear]

[brave crypt]

The condition number of the matrix was better ig?

Before it was reading like 10^16 and stuff

Now its consistently 40

[stone hill]

Can someone explain how we go from third to last to second to last? This is just a derivation of simpsons rule but I got no clue how we go between those two steps. Everything else is intuitive enough

Prof explained in class but I was doodling in the upper right corner the whole time so I forgot to note it down.

Test is tomorrow and we need to have this memorized(which i do, but id still like to know what happens there)

[wide spear]

Which two lines are you referring to

[stone hill]

Third to last and second to last

[wide spear]

Ok there are typos

It should be Aa^2+Ba+C and Ab^2+Bb+C

If that explains the confusion

[stone hill]

Ah yeah thats it thanks

Such a silly thing lol

[stuck oak]

Hello

I'm learning at university with a book called Numerical Analysis 9th edition, can I asked is there any where to find full solution of the book?

[wide spear]

If there are specific questions that you need help with, feel free to ask here

[mortal pier]

say that i have a python function Matrix(l1, l2) that creates a 2x2 matrix with eigenvalues $l_1, l_2$. Now i implemented the power iteration method and theory says that if the absolute values of $l_1 \text{ and } l_2$ are close to each other i will get slower convergence. Why do I get much slower convergence for Matrix(-(N+1), N) than Matrix(N+1, N) even though $|\frac{-(N+1))}{N}| = |\frac{N+1)}{N}|$ (here N is some natural number)?

[pine jetty]

dompa

[steady lantern]

Can someone help me in discretization of an ODE?

I can't seem to get it into set of linear equations.

It has quadratic terms.

[wide spear]

Post the ODE

Show what you've tried

[steady lantern]

Not sure how to work with quadratic terms

@wide spear

[wide spear]

By quadratic do you mean second derivative terms

Or is the (d phi/dx)^2

[steady lantern]

uh no, I get phi^2 terms

which I cant get it in a linear equation form to solve with gaussian elimination

diffusion coefficient varies with specific quantity, its annoying to solve

[steady lantern]

ok a little bit of research I did

and its called richardson iteration

to solve such systems

[fading glade]

hello, i have trouble with Allan Variation. After analyzing the variations and identifying the noises, I calculated the Bias instability of 0.06 deg/h, when according to the Datasheet it should be about 6 deg/h

Allan Var for MEMS gyro for all XYZ axis

[brave crypt]

Is your 0.06 measurement from averaging the x y and z plots

Think need more context from you regarding datasheet too

[wide spear]

[long lake]

Please go to one of the discussion channels for casual conversations

And please don't post the same message in multiple channels

[echo ferry]

I'm trying to determine (lol) what I can about the determinant of an $n\times n$ matrix (indexed from $0$ to $n-1$) whose entry in row $j$ column $k$ is $\sin(x_k)\cos(x_{j-k}) - \sin(x_{j-k})\cos(x_j)$.\
My advisor suggested somehow decomposing the matrix into the difference of two matrices, one with $\sin(x_k)\cos(x_{j-k})$ for the entries, and one with $\sin(x_{j-k})\cos(x_k)$ for the entries.\
Here, $x_0, ..., x_{n-1}$ are arbitrary real numbers corresponding to the arguments of a biunimodular sequence $(e^{ix_0}, ..., e^{ix_{n-1}})$.\
I'm not sure how to relate the determinant of the whole matrix to the smaller parts, since determinants tend to respect multiplication but not addition.

[pine jetty]

Tiborculosis

[wide spear]

You can rewrite that using trig identities right

This feels very close to a DFT matrix

[echo ferry]

Yes; I started with $\sin(x_k-x_{j-k})$ and was encouraged to expand and separate.

[pine jetty]

Tiborculosis

[wide spear]

Oh

I wouldn’t expand

[echo ferry]

It does feel like a DFT matrix. I'm trying to study the finiteness of circulant Hadamard matrices of given dimensions via an isolation argument, and trying to ascertain what conditions on dimension and $x_0,\ldots,x_{n-1}$ lead to isolation.

[pine jetty]

Tiborculosis

[wide spear]

How different from a DFT matrix is it

Is this a discrete sine transform matrix?

[echo ferry]

My advisor suggested expanding, which leads to a difference of matrices that have constant factors along the columns corresponding to $x_{k}$ for column $k$. It's very related to DFT matrices; there's a theorem that says that a circulant matrix $M$is a complex Hadamard matrix if and only if the first row of $M$ forms a biunimodular sequence, via conjugation by the DFT matrix.

[pine jetty]

Tiborculosis

[wide spear]

Expanding doesn’t tell you anything about the determinant

[echo ferry]

I can't see what one gains from doing so, but my advisor said something about using the multilinearity of the determinant to break it down further into two from here.

[wide spear]

I will experiment a bit later

[fathom rain]

how are $x_{1-n},\ldots,x_{-1}$ defined? $x_{-j}=x_j$?

[pine jetty]

se

[echo ferry]

Like so

Far as I'm concerned though, just real numbers between 0 and 2pi

I'm sure the other aspects of their definition are good for something

but i'm not sure what that something is.

[fathom rain]

so x_j are just random numbers, [0,2pi]? no ordering or anything?

[fathom rain]

if so, then I do not think there is much to do

[fathom rain]

where does this come from?

[still forum]

I'd like a hint on the following problem which is supposed to be easy but I blanked and it's 1am here

For functions $f \in C[0,1]$ with the usual supremum norm, define $E_n(f) = \operatorname{dist}(f, \mc{P}n)$ as the best approximation error to the function $f$ with degree $n$ polynomials. Prove that if $f^{(n)} > 0$ on $[-1,1]$ then $E{n-1}(f) > E_n(f)$

[pine jetty]

jamiecjx

[wide spear]

What have you tried

[still forum]

Ok, after waking up I still didn't figure it out, but I tried mostly combinations of Chebyshev equioscillation theorem together with some properties of divided differences in order to use the condition on $f^{(n)}$

[pine jetty]

jamiecjx

[still forum]

I think I figured it out, but it's longer than I expected and I hope there's a smaller solution

Let $p \in \mc{P}{n-1}$ such that $E{n-1}(f) = E_n(f)$. Then $p$ is the best approximation to $f$ in $\mc{P}n$ as well so there exists a strictly increasing sequence of $n+2$ points $x_1, ... ,x{n+2} \in [-1,1]$ such that

\begin{equation}
f(x_i) -p(x_i) = (-1)^i \gamma \qquad |\gamma| = E_n(f)
\end{equation}

By considering divided differences we have $(f-p)[x_1,..., x_{n+1}] = f[x_1,..., x_{n+1}] > 0$ by the MVT for divided differences (since $f^{(n)} > 0$). Similarly $(f-p)[x_2,..., x_{n+2}] = f[x_2,..., x_{n+2}] > 0$.

But we also have the explicit formulas

\begin{equation}
(f-p)[x_1,..., x_{n+1}] = \sum_{j=1}^{n+1} \frac{(-1)^j \gamma}{w_1'(x_j)} \qquad w_1(x) = (x-x_1)...(x-x_{n+1})
\end{equation}
\begin{equation}
(f-p)[x_2,..., x_{n+2}] = \sum_{j=2}^{n+2} \frac{(-1)^j \gamma}{w_2'(x_j)} \qquad w_2(x) = (x-x_2)...(x-x_{n+2})
\end{equation}
Note that since $w_1'(x_2) < 0$ and $w_2'(x_2) > 0$ and the signs of $w_1'$ and $w_2'$ alternate between each root, it's not possible for both of these sums to be positive, hence contradiction

[pine jetty]

jamiecjx

[fathom otter]

Hello guys, I'm new here. I have a question regarding an error estimate for the explicit Euler method for a 2D heat-equation. It's on channel #help-4 . Can anybody give it a look?

[wide spear]

Please repost your question here

[fathom otter]

Sure, here goes:

Hello all, I'm new in the server!

I just posted this on Computational Science Stack Exchange, but decided to try my luck here as well.

I am solving the heat-equation $\frac{\partial u}{\partial t} = \alpha\nabla^2 u$ on the domain $\Omega = [0,1]^2$ and interval $t\in[0,1\times10^{-2}]$ with homogeneous Dirichlet boundary conditions in 2D. I am using a central finite-difference method to discretize in space and the forward Euler method to discretize in time, and end up with a scheme like:

\begin{equation}
u^{n+1} = u^n + \delta t\alpha L u^n,
\end{equation}
where $u^n$ are the nodal values of the numerical solution and $L$ is the discrete analog of the Laplacian operator in 2D.

If I compute the relative global error as $e \colon = \frac{||u^T-u^\text{exact}(t_f)||_2}{||u^\text{exact}(t_f)||_2}$ and plot it against the mesh spacing $h$, I get the following plots for $\delta t = 5\times 10^{-4}$ and $\delta t = 1\times 10^{-4}$ (plots in the attached figures)

The red dots correspond to simulations done outside the CFL stability condition $dt \leq \frac{1}{4\alpha h^2}$.

I have \textbf{two questions} about these graphs.

\section{Question 1}
Why isn't the error blowing up sooner (for larger $h$) in the case of $\delta t = 5\times 10^{-4}$, since the points in red all violate the CFL condition? Might it be a sign of a problem in my implementation?

\section{Question 2}
I find this to be weird behaviour since I would expect the error to be a monotonic increasing function of $h$, instead it seems that there is an optimal $h$ for which the error is minimal.

An explanation I found for this behaviour is that, as $h$ decreases, the condition number of $L$ increases, namely $\frac{\lambda_\text{max}}{\lambda_\text{min}} \approx h^{-2}$. This causes the ODE system to become increasingly stiffer, and thus the local truncation error for an explicit method, such as the Euler method, becomes increasingly large. Is this the correct explanation for this behaviour? If so, can someone give me a more quantitative explanation? In specific, can I predict where the "inverted peak" will happen?

[pine jetty]

Rui Martins

[wide spear]

I might have time to answer you today

[fathom otter]

Thank you @wide spear !

I got several replies from the Computational Science Stack Exchange and already understand what is happening. Here is the link of the question:

https://scicomp.stackexchange.com/questions/45072/euler-method-error-estimate-for-heat-equation-for-constant-time-step

It's basically the leading terms of the local truncation error cancelling out, for a particular choice of $N$, that casuses this behaviour.

[pine jetty]

Rui Martins

[wide spear]

Ok nice

[slender fiber]

Do you happen to be familiar in this topic?

I've been having issues with a model converging near the center of a circular radius, tends to blow up with a singularity and I'm trying to find numerical/mathematical methods to correct it without using blind assumptions.

This paper is certainly useful, though it tickles my brain to read its math rhetoric.

[brave crypt]

My code is not working. I just want to plot the solution to my ODE. This seems like the exact same thing as in the ODE45 documentation. What's wrong with this?

[wind patrol]

Can anyone help me solve or at least find a good general approximation for this system of differential equations with initial conditions $v_{i;x}(0)=\cos{\alpha}$ and $v_{i;y}(0)=\sin{\alpha}$ : $\frac{dv_x}{dt} &= -\frac{k}{m} \sqrt{v_x^2 + v_y^2} \cdot v_x \$, $\frac{dv_y}{dt} &= -g - \frac{k}{m} \sqrt{v_x^2 + v_y^2} \cdot v_y$? Thanks!

[pine jetty]

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

[vapid crystal]

ok i believe this is the right channel
so we are learning about approximating roots of functions using a bunch of different types of methods
one of the exercises asks me what the convergence order of the newton method is
in the text book it says p=2 if it is a singular root. So im wondering what it is if it is a multiple root

[wide spear]

More details

[wide spear]

#computing-software

[wide spear]

#odes-and-pdes

[wide spear]

Have you looked at the proof of convergence for newton's method

[vapid crystal]

ah thank you there is in fact a proof abt a multiple root a few pages before where i was looking

[wide spear]

Nice

[vapid crystal]

what is the order of a perturbation?
there is a sentence in the book that says
"Since the perturbation of $f$ is of the order 0.001. We do indeed expect that $\Delta x = 0.001$

[pine jetty]

you_are_me

[vapid crystal]

We have the function $f(x)=\sin (x)$ and the perturbation function is $g(x)=0.001(1+\sin(x))$

[pine jetty]

you_are_me

[vapid crystal]

and the condition number is 1 btw

[wide spear]

Order = size in this context

[vapid crystal]

aha so like if i make the number 0.003 instead of 0.001 then it is of order 0.003?
that makes sense

ok and so to find $\Delta x$ u just do the order times the condition number?

[pine jetty]

you_are_me

[vapid crystal]

i think that makes sense in my head

nvm that doesnt work in the next exercise 💀

[wide spear]

Well, usually by order you mean a power of 10

0.003 is still of the order 10^-3

[vapid crystal]

ooooh yeah obv
my brain is deepfried ngl

ok then everything works out yipie

Thanks btw before i forget to say that

[wide spear]

Glad to help

[lime smelt]

who here is a finite element method chad

[slender fiber]

I can provide a bit more when I sit down at my office computer. Short story is that I have a momentum balance model that iteratively calculates the balance of forces from an outer radius to the inner point of a circle. Each iteration, it numerically calculates some (difficult to solve explicitly) integrals, of which are a function f(y)= Annular area / annular outer perimeter and f(y) = 1 / annular outer perimeter. Both cases, especially the latter, approach infinity as each iterative calculation approaches the centerline, ie singularity.

I know what the final solution profile z = f(y) should be - it should approximately fit a 1/7th velocity power law profile except with a 0 slope at the centerline.

I was curious about numerical methods to a) identify when and if a singularity effect begins/is significant and b) how to stop it from happening numerically.

So far, I’ve been just finding when the solution z = f(y) second derivative changes signs because that would keep the (slope) first derivative from approaching zero at the centerline, and then approximating a polynomial fit. But, it doesn’t feel verifiable/tied to the model.

[lime smelt]

that was my question

regards

[slender fiber]

As per, I'm using shear from radial distance y = d to y = R, where d is the radius of the film and R is the inner radius or centerline of the circular cross-section. See attached image. We assume angular symmetry. Marching from y=d to R, we iteratively calculate the velocity from a shear balance.

The velocity equation is shown in an image from one of my model sketchbooks. The integral terms in the equation likely have to be solved numerically, as the 1 / mu_g_t term essentially results in a VERY nasty integral. But, the problem lies in the A_H / P_H,y and 1 / P_H,y terms in the integrals. In my case, A_H = pi ((R - d)^2 - (R - y)^2) and P_H,y = 2 pi (R - y).

These integral terms blow up, particularly the integral of 1 / P_H,y, as y approaches R simply due it becoming the integral of infinity. Are there numerical methods to avoid such a singularity, or to help identify when the singularity becomes significant?

Also, the colored graphs shows the terms of the equation calculated in terms of [m/s] velocity contributions, i.e. you can see them blowing up as they approach the centerline radius R.

[wide spear]

Is radial distance the distance to R?

Can you latex the velocity equation

[slender fiber]

\begin{align*}
u_{g}(y) &= -\frac{dP}{dz} \int_{d}^{y} \frac{A_{H}}{\mu_{gt} P_{Hy}} , dy

• \tau_i P_{fd} \int_{d}^{y} \frac{1}{\mu_g P_{Hy}} , dy
• \Gamma_{fg} u_i \int_{d}^{y} \frac{1}{\mu_g P_{Hy}} , dy \
  &\quad + \frac{d}{dz} \left( \left[ \int_{d}^{y} \rho_g u_g^2 P_{cy} , dy \right]
  \left[ \int_{d}^{y} \frac{1}{\mu_{gt}} P_{Hy} , dy \right] \right)
  \end{align*}

[pine jetty]

kryle

[slender fiber]

Here’s the equation!

[slender fiber]

It’s not in the image, but y is the radial distance from the outside wall, from y=0 at the exterior to y=R at the centerline

[green willow]

regarding numerical solutions to partial differential equations:

i hate numerical solutions to partial differential equations. that is all

[grave spoke]

why?

[green willow]

a dozen little inequalities and bounds chained together in a not-entirely-intuitive for the simplest elliptic equation on a square domain with as much continuity as possible; and then as soon as you want to consider a less smooth function it gets three times worse, and if the domain isn't square then it's even more terrible, and if the equation isn't simple then god help you

and all that's just for elliptic equations, for parabolic or hyperbolic you want a completely different set of tricks and inequalities and you have to worry about stability analysis and forwards and backwards schemes and oh my god

if nothing else the sheer amount of writing is going to kill me, the equations take three lines and every time i find a mistake somewhere on page 1 another tree dies

[grave spoke]

Are you talking about the finite difference method?

[green willow]

yes?

[wide spear]

Oh

Well simply do not care about error bounds

[green willow]

good idea, i'll try that in the exam

[slender fiber]

Go ahead, find another explicit solution

[lime smelt]

solved laplaces equation in 2D with boundary condition on an arbitrary mesh with finite element method, the height of a point is equal to the value of the function, I think I cooked chat

[wide spear]

Ok

[lime smelt]

[wide spear]

Have you checked for convergence

[lime smelt]

yeah I did I calculated the solution analytically aswell for some simpler meshes like square shape and f(x,y) sin(x)*e^y and the average error per node was like 0.002 or something like that if I compared my approximation u with boundary condition f and f itself at those nodes

like if thats my boundary condition then I would expect u = f on the interior aswell

[lime smelt]

are you proud of me?

[wide spear]

That's not what it means to check for convergence

[lime smelt]

what does it mean

[wide spear]

Change the spatial resolution

Compute the relative L^2 error for each

Check if it converges 1st order or 2nd order or whatever

[lime smelt]

what range of mesh size should I use? 0.1 to 0.0001 ? how many zeroes?

I need to refactor my code to use sparse matrices because now its just slow, I'll try checking afterwards thanks

[wide spear]

Sure

A wide range of resolutions

[molten knot]

i like to just loop through my solver, divide by 2 in the mesh size each time, compute the ratio of the previous error and the current error, and the log base 2 of it

if you have half the mesh size and the error reduces by a factor of 4 then you’ll have second order convergence

[radiant salmon]

I am currently learning about Spectral Methods to solve PDEs numerically and I am would like to ask about it’s accuracy/stability

Spectral methods consider the points/variables globally and does something akin to interpolation to get the function/model of best fit. From my experience w/ normal interpolation, I know methods like cubic splines often are far more stable

I know Spectral methods have exponential like convergence already but how would the method fair when mixed with splines?

I’m still learning all this as I’m asking and its been awhile since i learned interpolation so apologies if my question seems rather vague

[wide spear]

There are spectral element methods where you do something local

Frequently it's just with a standard set of polynomials but you can use splines as well

[radiant salmon]

Could you elaborate?

In what part of the model would the splines take place? Would it be after you found the basis functions?

[wide spear]

The splines would form a local basis

Instead of a global one

You can look into spectral element methods

[rare schooner]

ayo hi let me poke my head in here. i got a for-fun algorithm that i want to analyze for numerical precision

im thinking that maybe i can expect a worst case error of about log_10(1.00005) per computational step? which means 10 times that in all, so expecting 4dp is a bit strenuous...?

here's an example of a computation with this algo

[wide spear]

log_10(1.1) rounded to 4 decimal places is 0.0414 so you aren't that far off

Absolute error of 0.0004 so relative error of 0.0097

So for 12 steps, you introduce roughly 0.0008 relative error at each step

[wide spear]

<@&268886789983436800>

In many channels

[vapid crystal]

so the question asks
the newton cotes formula iwith even degree n=2k is as accurate as the newton cotes formula with uneven degree n=2k+1

im looking at the error upper bound formula for n=2 and n=3 and i see in both a fourth derrivative and h^5 but for n=3 it has a term 3/80 while for n=2 it has a term 1/90

so i conclude the upper bound for n=3 is in fact higher? thus meaning that the statement is wrong and in fact n=2k+1 is less accurate than n=2k

but im seriously doubting that so can someone tell me if i went wrong and how i went wrong i i did go wrong?

oh wait i think i see where i went wrong h is smaller for 3 than for 2 right? so 3 would be more accurate

[brave crypt]

how would i solve this numerically

[brave crypt]

You can use Python for it

[brave crypt]

can you do it for me pls i dont have it installed rn

but i need to do it by hand as well

[wide spear]

[brave crypt]

no?

like what method do i use

and how

idk how to use python for this

[wide spear]

Newton's method

[wide spear]

Yes

[fathom otter]

You can use Newton's method as @wide spear . Other non-linear solvers will also work.

I don't use python, but I suppose you first have to write your system as F(x) = 0, where x is (b_1,b_2) in your case.

Something like this:

$F \left(b_1,b_2\right) = \begin{bmatrix}\cos(b_1)-0.72\cos(b_2) \tan(b_1)-\tan(b_2)+0.8697\end{bmatrix}$

[pine jetty]

Rui Martins

[fathom otter]

I think you can use fsolve from scipy.optimize, but there might be a better option. If you just want a solution for this problem, any nonlinear solver will do really.

[pine jetty]

No selfroles matching not studying.
See ,selfroles --list for the list of valid selfroles.

Jun 18, 2025, 2:25 PM

Please select the desired selfroles! (Use ,iamnot to remove your current selfroles.)

1.   Helpers
2.   Not Very Ppl
3.   Talks
4.   Role that does nothing
5.   she/her
6.   they/them
7.   any pronouns
8.   ask pronouns```
Please type the numbers corresponding to your selection, separated by commas, or type `c` now to cancel. (E.g. `2, 3, 5, 7, 11`)

Removed the studying! role from you.

Selfrole selector timed out, your roles were not updated.

Jun 18, 2025, 2:27 PM

[wide spear]

<@&268886789983436800>

[mystic umbra]

Was something perhaps deleted?

[fluid sedge]

I'll take a look, there is a bit of a delay

Just some bot command spam it seems. @unkempt warren please use #bots for bot spam and never delete messages once a mod ping is made.

[late kettle]

This Galerkin approximation is using the same basis for test functions?

[wide spear]

Yes the phi_i in 1.5.5 and 1.5.6 are the same as in 1.5.1

[late kettle]

Thank you!

[next garden]

I am doing a thing with legendre polynomials, for that I need the first derivative. I found this SO question:

https://math.stackexchange.com/questions/4751256/first-derivative-of-legendre-polynomial

Which states that the first derivativeof the basis element n is

$P'n = \frac{n(P{n-1} - xP_n)}{1 - x^2}$

And that cannot be right.

First basis element is 1 second is x

So derivatives are 0 and 1 respectively. Even assuming that if you assume that $P_{-1} = 0$ neither derivative fits the statement

For x the formula would tell you $(-nx / (1- x^2))$ which trivially cannot give you 1

[pine jetty]

Makogan

[next garden]

So what am I doing wrong?

[wide spear]

The formula holds for n>=2

[next garden]

do you know of a place that has the proof in a straightforward manner?

[wide spear]

The MSE post has the proof in the answer

[next garden]

It assumes background context I don't have, such as teh generating function

[wide spear]

What does that even mean

[next garden]

exactly!

[wide spear]

I don't know what you mean

[next garden]

Idk where this is coming from

[wide spear]

Read any standard treatment of legendre polynomials

[next garden]

The answer in MSE relies on something called the generating function. The definition of Legendre polynomials I am working from is a recursive definition

[wide spear]

The recurrence relation is derived from the generating function

[wide spear]

This is very normal

[next garden]

Yes but I have no formal training in any of this and no one to ask direct questions to. For example, from the wikipedia entry:
https://en.wikipedia.org/wiki/Legendre_polynomials#Definition_via_generating_function

That statement seems arbtirary.

The first two terms in the sum are $P_0(x) t^0 + P_1(x) t^1$

How the text jumps from there to asserting that $P_0 = 1$ and $P_1 = x$ doesn;t seem like a given based on the defnition.

I don't doubt it;s correct, but I don't see why it is

[wide spear]

Do you know how generating functions work

[next garden]

All I have is the definition from the wikipedia link. As mentioned I have no formal training in any of this.

[wide spear]

You are, of course, welcome to ask direct questions here

[next garden]

I do not know how generating functions work

[wide spear]

To get $P_n$ from the generating function, take $\frac{1}{\sqrt{1-2xt+t^2}}$ and perform the Taylor series expansion in $t$

[pine jetty]

Angetenar

[wide spear]

First 5 terms for example

The coefficient of $t^n$ is $P_n(x)$

[pine jetty]

Angetenar

[next garden]

taylor series centered where? Mclaurin series?

[wide spear]

Centered at 0

[next garden]

ok so the Mclaurin series, thank you, let me think a bit

follow up question, why is the generating series the generating series?

i.e. why is that the specific expression that was picked to generate the legendre polynomials?

[wide spear]

https://en.wikipedia.org/wiki/Legendre_polynomials#Expanding_an_inverse_distance_potential

Gravity

[next garden]

This answer has a lot of weight

(Thank you)

[wide spear]

[next garden]

Follow up, trying to understand how they are deriving this step from the prior:

why are those two expressions equal to each other

oh I think I get ti

the equality follows from algebraically manipulating the generating funciton and its derivative

so that they look the same

[wide spear]

Yes

[halcyon nacelle]

It arises naturally in calculating the integral kernels of the orthogonal projection operators onto eigenspaces of $-\Delta$ on $S^n$ (Legendre is $n = 2$). See page 130 of https://link.springer.com/book/10.1007/978-3-031-33700-0.

[pine jetty]

L

[next garden]

Sorry for the spam. Stuck on this step:

I don't see how that derivation follows form the antecedent

[wide spear]

Distribute the 1-2xt+t^2, equate equal power of t

[next garden]

Ah by the fact that polynomials are euqal iff their coefficients are euqal

[next garden]

Still stuck here. I managed to get three of the terms, not the last one.

I get where the RHS comes from.

I get where most of the lehs equality comes from.

But it seems that I should get a non zero residual?

i.e. I distribute (1-2xt+t^2) giving me 4 sums.

Two of those sums have matching powers of t^n+1, equality is trivial.

Then one sum is overt t^n+2 and one over t^n

So most terms in the polynomial over t match by shifting every two. That gives (P'_(n+1) + P'(n-1)) t^(n+1) to equate to the rhs.

But then there's some elements in the lhs with no match?

oh this is abusing the fact the sum has countably many elements?

Right that;s why I seem to have a residual

This is just hilbert's hotel

yes?

[next garden]

Final question tha tis a bit tangential to all of this.

I know the quadrature rule for Legendre polynomials:
https://en.wikipedia.org/wiki/Gauss–Legendre_quadrature

But it feels that if I am doing the inner product, i.e.:

$\int P_n(x) f(x) dx$

For each basis polynomial (i.e. extracting the coefficients of the hilbert space) there has to be a simplification with the $w_i$, no?

[pine jetty]

Makogan

[wide spear]

You're doing this on [-1,1] right

If you do a method where you break up [-1,1] into many sub intervals, then no simplification occurs

If you do it on the entire interval, then if P_n matches the quadrature degree chosen, the integral for P_n vanishes

[next garden]

why would it vanish? is it not dependent on what f actually is?

say for example taht f = 1/x. that would lower the degree of each monomial, that should yield a different answer, no?

[wide spear]

If you are doing a degree $n$ Gauss-Legendre quadrature, then you approximate $\int_{-1}^1 P_n(x)f(x)dx\approx\sum_{i=1}^nw_iP_n(x_i)f(x_i)$ and $x_i$ are the roots of $P_n$ so $P_n(x_i)=0$

[pine jetty]

Angetenar

[next garden]

Oh fuck, yeah you are right...

How on earth do you do that projection then?

i.e. what is an effective way of projecting f into the legendre basis?

All the paper I have on hand says is this:

Which is not terribly useful

The goal is to extract the coefficients of the Legendre basis:

I am getting more and more tempted to just do the naive sum and delta

but that's gonna be highly inaccurate

[wide spear]

Yes quadrature degree 4p

So that you don't have the above problem

You can also do something like a midpoint rule or simpson's rule and break up the interval [-1,1] into many subintervals

[next garden]

if I use simpson's rule then it;s no longer a gaussian quadrature is it?

[wide spear]

Do you need to use a gaussian quadrature

[next garden]

idk, the paper uses that and I am reticent to deviate

this is part of an algorithm that;s going to be approximating very complex functions

[wide spear]

Well the paper is using a gaussian quadrature of degree 4p

[next garden]

so I want to be both accurate and fast

[wide spear]

While you only compute legendre coefficients up to degree p

So you don't have the cancellation problem

[next garden]

I am probably miss understanding both you and the paper.

What I was trying to do is $\sum ^n_{i=1} w_i f(x_i) P_n(x_i)$

with n = 4p

[pine jetty]

Makogan

[next garden]

which as you pointed out is wrong

[wide spear]

At any rate, I think that adaptive quadratures are better for complicated functions

You compute $<P_n,f>=\int_{-1}^1P_n(x)f(x)dx$ for $n=1,\ldots,N$. The paper approximates this with $\int_{-1}^1P_n(x)f(x)dx\approx\sum_{i=1}^{4N}w_iP_n(x_i)f(x_i)$ where $x_i$ are the roots of $P_{4N}$

[pine jetty]

Angetenar

[wide spear]

Be careful about the indices

[next garden]

You are probably correct and that would also lead to better handling of discontinuitiees. I am probalby going to go with Simpson;s rule due to simplicity of coding.

Let me read

Oh!

I said nothing

I get it

Thank you

Final question, this is going to be approximating highly complex functions, like the ones int he image. And I need to try to update the coefficients as fast as possible.

Simpson;s rule is nice in that it's dumb simple to code, but I worry about eneding a large number of sub intervals to get good accuracy.

[wide spear]

Use adaptive intervals

Instead of fixed size intervals

[next garden]

how do I pick the interval?

[wide spear]

Start with a base spacing Delta x

In each interval of size Delta x, estimate the error of the quadrature method

If the error > some threshold, split the interval in 2

[next garden]

Hmmm, that might be too computationally expensive. I will read a bit more into other numeric integration methods and see what I get.

Thank you enormously for patiently helping me btw

I appreciate it

[next garden]

Follow up, I know that there exists a tridiagonal matrix whose eigenvalues yield the nodes of the quadrature rule of order n.

But I have been unable to find a resource stating what that matrix is

does anyone know?

[wide spear]

Golub-welsch algorithm

[next garden]

thank you!

[next garden]

Am I correct that, in the Legendre basis, $b_j$ is 0?

[pine jetty]

Makogan

[next garden]

For the Golub-Welsch algorithm?

The resulting matrix would have 0's in the diagonal if that's the case, so I am doubting myself

[fathom otter]

In their article, they relate $b_j$ and $\alpha_j$ with:

Then $\alpha_i$ can be computed with 4.3:

And $r_ij$ are the entries of the cholesky decomposition matrix of what they call the Gram matrix

For bj to be zero, alphaj needs to be zero. But for me its not obvious why this would be the case, at least looking to this problem from this perspective. What makes you think that the bj coefficients are zero?

[fathom otter]

There also might be a simpler way to compute the roots of the Legendre Polynomials. I think I did this once and computed the roots with newton's method and the weights have a formula on wikipedia.

[next garden]

I trust the eigenvalue method more. It doesn't ned an initial guess

[next garden]

I am using this writeup since it's more condensed:
https://paulklein.ca/newsite/teaching/cookbookintegration.pdf

but if I am reading things right the b_j here would be 0

[fathom otter]

For the Gauss-Legendre quadrature, in section 4, indeed it seems so that bj =0, but I'm not sure why.

If I have some time tomorrow, I can take a look at this.

But im guessing you'd need to take expression 4.3 of Golub's article and work out why the two terms of the lhs cancel each other

By using the definition of r_ij for the case of the Legendre Polynomials.

[next garden]

at least we are in agreement

[next garden]

thank you so much either way. i appreciate the help

[next garden]

I can confirm, experimentally that the bj are 0

[wide spear]

We are glad to announce that the University of Pisa has launched a new interdisciplinary industrial PhD program in High-Performance Scientific Computing (HPSC) with 4 fully funded positions open for applications (EU and non-EU), starting November 2025.

This innovative doctoral program is designed in collaboration with the medical technology company Sordina IORT Technologies, and brings together 8 university departments and several national research institutes. It focuses on the development of advanced computational methods and HPC-based solutions for complex problems in science, engineering, and medicine.

Specific research topics for the funded scholarships include:
 - Iterative Methods and Preconditioners for Sparse Linear Systems on Pre-Exascale and Exascale Architectures
 - HPC Software for Design and Analysis of Electronic Devices using Advanced Materials and AI/ML
 - Computational Models for Flash Radiotherapy and Radiobiological Applications (2 positions)
 - Reduced-Precision and Matrix Units on AI GPUs for Efficient Solution of the Wave Equation

The program offers a unique interdisciplinary environment where PhD candidates will develop expertise in:
 - Numerical algorithms and scientific computing
 - HPC architectures and parallel programming
 - AI/ML techniques applied to physical and biomedical systems
 - Application domains such as climate modeling, advanced materials, and smart medical technologies

This is an industrial PhD: selected candidates will have the opportunity to work closely with Sordina IORT Technologies on translational research projects in medical computing and radiotherapy, ensuring both academic excellence and real-world impact.


Application deadline: July 18, 2025


Program website and full call: www.dm.unipi.it/phd-hpsc

We welcome applications from outstanding Master's graduates in mathematics, computer science, physics, engineering, chemistry, and related fields.```

[wide spear]

Hello

Is anyone familiar with clenshaw-curtis quadrature

I've found two online sources which present two different sets of weights

https://numbersandshapes.net/posts/high_precision_clenshaw_curtis/

https://mathsfromnothing.au/clenshaw-curtis-quadrature-rules/?i=1

These weights are not the same!

Oh wait

Lol

These are for different sets of nodes

Lol

Nevermind

[wide folio]

Hi
Would like some advice or inputs on how I can tackle this issue I've been facing

Now I need to solve the 2d viscous burgers equation using a semi lagrangian method for the material derivative and fem for the viscous term.

Now the problem I've been facing is this: I do not have a velocity function for calculating the departure points. I have only the initial and boundary conditions available to me and I need the velocity function for me to be able to use rk4 or a suitable time integration scheme to calculate the departure points.

How do I do this😔

[wide spear]

What do you mean you don’t have the velocity

[wide folio]

😃 I've figured it out pls ignore

[wide spear]

[wide folio]

I was the problem. It was me 😩😔

[wide spear]

[next garden]

This is perhaps a dumb question but my brain is fried.

How can I generate all tuples of integers representing the basis of a 3dimensional space of polynomials of degree p.

i.e. a tuple of three integers such that i + j + k < p.

For exmaple for p = 2

1, x, y, z, xy, zx, yz, z^2, x^2, z^2

[wide spear]

Three nested loops

[next garden]

This right?

        for i in 0..=degree
        {
            for j in 0..(degree - i)
            {
                for k in 0..(degree - i - j)
                {}
            }
        }

[wide spear]

Yes

[next garden]

ty

[wide spear]

Which discrete cosine transform does this matrix correspond to

[wide spear]

Disregard this matrix is written incorrectly

[wide spear]

Sobbing

Unable to use a jacobian to do a spherical quadrature

[vapid plume]

Hello, currently studying FEM and looks like I need to study some things in detail.

Does anyone have any book recommendations for:

1-) Explanation of the mesh generation process using triangular/quadrilateral elements in 2 & 3D,

2-) FEM using MATLAB which makes use of mesh generator functions and solving of problems with large elements so that I may better understand the whole process (boundary conditions, mesh generation, etc).

The book I'm following only solves problems for 1-2 elements so everything is done manually. It is great for a beginner in FEM though but I believe it is time to read other books as well.

Thank you.

[wide spear]

1. https://onlinelibrary.wiley.com/doi/abs/10.1002/pse.135
2. Is this more about the implementation

Anyways also maybe read Finite Element Analysis: Method, Verification and Validation by Szabó and Babuška

[topaz vale]

could someone explain what analysis is about 💀

[wide spear]

#math-discussion

[vapid plume]

Yes to the implementation part, and thank you for the reply.

[wide spear]

https://finite-element.github.io/

[hardy oak]

This is from Brenner & Scott, which lacks an errata

Is this problem posed right? This definition of phi-epsilon just doesn't seem right

phi^epsilon has completely unrelated support?

and this support changes over epsilon

further, I believe that claim of the identity is untrue

but I can't verify it

[next garden]

Equation 5 in this paper is wrong right? the Gauss-Legendre basis is merely orthogonal , not orthonormal. So that term needs to be normalized, right?

[wide spear]

They are working with a normalized basis

Sentence before eq 4

[next garden]

I should learn to read more carefully

[fathom otter]

About 2, search for this book:
The Finite Element Method: Theory, Implementation, and Applications

This is the book I used when I was learning FEM and I liked it.

They give code snippets in matlab

[vapid plume]

@fathom otter Thanks.

[wide spear]

https://arxiv.org/abs/2506.21654

[bold yew]

Does anyone have a good reference on stability for the Newmark-beta method? In the engineering texts they are only stating the stability without any proofs.

[next garden]

Question, is there a relatively simple way for me to test the numeric coefficients of projecting a torus SDF onto a multidimensional Gauss-Legendre basis?

I coded that exact thing and the results I am getting are not right. So I need to figure out what the problem is. But if I have to code something very complex to compare coefficients, might as well just debug the code I already have.

like, I have the definition of everything math wise, I just want to know if there is a sympy based library or something of the sort I can quickly spin to check if I get the same numbers

[wide spear]

Hand compute some values

[next garden]

oof...

That has to fall under some clause of the Geneva conventtions

[wide spear]

Lol yes

The worst hand computation I've had to do was proving RK4 for vector functions

[wanton skiff]

Can someone help me with the exercise 14.3? I couldn't do it.

[minor osprey]

What are you having trouble with in particular? What have you started, what have you tried, etc.

[wanton skiff]

I’m looking for the answers for exercise 14.3, the solutions for it.

[minor osprey]

This isn't really the place to just get answers to problems

which is why I'm asking about what you're stuck on so I can help

[wanton skiff]

Can you do it?

[minor osprey]

Sure, I can solve it, but I won't just give you the solution

There's nothing to be learned there, but also I'm lazy and don't want to write the code associated with it on my phone

[wanton skiff]

I don’t trust you, you will use ai

This is university mathematics, specifically used for master degrees or bachelor.

[minor osprey]

I'm sorry you feel that I would do that

I actively avoid llms in my daily life lol

[minor osprey]

but to get back on topic, what about this problem is causing you difficulty?

[vapid plume]

@wanton skiff Best to take a 3x3 matrix and do everything by hand, you'll notice how you use the function (which you have to code) over and over again

[wanton skiff]

Thank you

[wanton skiff]

Will do, of course.

[shut widget]

Hey what is the difference between numerical analysis and real analysis or what is numerical analysis only

[minor osprey]

Numerical analysis is essentially studying how to practically compute objects often seen in analysis and linear algebra while understanding how error propagates from approximation.

Many problems in analysis have solutions that are continuous, but computers only work in a discrete setting so ways to get approximations are fairly important

[buoyant fiber]

Hello. I have a system with roughly 35 variables that need to be determined. For now I am using the monte carlo markov chain to try to get them. The problem is, there are first six conditions that must be satisfied and if they arent the code return a likelihood equal to -inf. I know that this is a high dimensional problem, but I already got millions and millions of attempts, more than 40Gb of data and all the values return with a likelihood of -inf. How can I change the way I am doing this to improve the results or I just need to wait for the code to find one point where it actually passes through all 6 conditions?

[wide spear]

Are you trying to solve a system

Or optimize

Is it linear or nonlinear

What form do the constraints take

[celest inlet]

Hello i need help

Im using a modified newton's method to find a root for initial values p0 = 0.35, and p0 = 0.95
for 0.35 i found the root to be 0.5, and 0.95 i found the root to be 1

however after asking AI i got confused since they mixed answers for p0 = 0.95 because 1 can or cannot be a root or they called it a false root because it would make f'(x) = 0 im looking for some explanation for this

[hard flare]

I probably wouldn't recommend asking AI for help with this. What do you mean by "1 can or cannot be a root"?

[celest inlet]

i was asking deepseek and chatgpt to compare and both gave differing results.
chatgpt said p0 = 0.95 does not diverge
deepseek said p0 = 0.95 converges to 1
so i asked chatgpt about it, it said it is converging to a false root

[hard flare]

Yeah I mean, you've experienced firsthand the reason I wouldn't recommend asking AI - it's not a search engine and it's not gonna give you consistent/good answers. Your Newton's method told you there's a root at 1, and plotting the original function shows that's the case. It happens to also be true that f'(1) = 0, but that doesn't change the fact that f(1)=0

I guess I'm just a little confused on what you need help with - were you trying to get the AI to do the convergence analysis part for you, or just verify that your code worked?

[celest inlet]

Both kinda, i tried to get the AI to do convergence analysis but i also wanted to see if my code worked (also made by ai)
I kinda dont know what to think whether p0 = 0.95 should be converging or not or diverging or whats suppose to happen. I guess im kind of trying to like check my work

[hard flare]

So you're using AI to check the work that was also done by AI? Lmfao

I think you would probably be able to understand the assignment better if you didn't outsource all of the work for it

[celest inlet]

yep

[hard flare]

Gotcha. Well, if you try it again yourself I'll be happy to help you figure it out, but I don't really want to help ChatGPT do your homework for you. Good luck though

[celest inlet]

Thank you for trying to help

[next garden]

Very naive question. I ahve a quadrature rule defined on [-1, 1]. And I am using it to project a function onto a function basis.

Say $\int^1_{-1} f(x) \phi_i(x) dx$

[pine jetty]

Makogan

[next garden]

Since my basis is defined on that interval, my quadrature rule becomes:

$ \sum w_j f(x_j)\phi_i(x_j) $

Now, say that f is instead defined on some arbtirary interval [a, b] and I want to do the same procedure

How should I modify my quadrature values?

[wide spear]

w_j becomes (b-a)/2 w_j

x_j -> (b+a)/2 + (b-a)/2 x_j or something

[next garden]

for both?

or only for f?

sorry replied to the wrong comment

x_j -> (b+a)/2 + (b-a)/2 x_j

Do I normalize the nodes only for f or do I also normalize them for phi?

[vapid plume]

@next garden Is this in FEM context? just curious. Btw yes

[next garden]

Surprisingly no, although that;s what I usually do. The actual algorithm is trying to simplify an SDF by decomposing it into an oct tree of Legendre basis

So I am not solving a differential equation or anything of the sort

[vapid plume]

Ah ok.

Whats SDF?

[next garden]

I am just changing the represnetation of an SDF into a nicer one for a computer toe valuate

[next garden]

Signed distance function

It's an implicit representation of a surface

[vapid plume]

I see

[next garden]

The english definition is, given a surface, the SDF gives you the euclidean distance from any point on the space to its closest mathc on the surface

except that the distance is positive if the point is outside the shape

and negative if the point is inside the shape

It;s extremely popular for physics simulations

paritcularly fluids and raytracing

[vapid plume]

I was just reading on how meshes are generated for FEM and I did come across sdfs, okay. Forgot the abbreviation.

[next garden]

because you can easily comput things like penetration depth and pushback forces

I am however not gettign what I expect from my implementation

I tried projecting the function 1 into the legendre basis of [-1. 0]^3

and I am getting non zeros for the coefficients of the non constant legendre vectors

so clearly

my implementaiton is borked

:C

[vapid plume]

what do you mean by this? replacing x in the function with the transformation?

[next garden]

So, by default, the legnedre basis is defined in [-1,1]

You can extend that to multiple dimensions by making a tensor product

that gives oyu a basis in [-1,1]^3

i.e. the cube centered at teh origin with sidelength of 2

that part works

I tested that basis with many toy polynomials and the projection gives me the same as the ground truth

But if I try any other domain (and trying to shift and re-scale)

so e.g. if I try to define eveyrhting to projet properly on [-1, 0]^3

It no longer works

despite me trying to scale things appropriately

[next garden]

So follow up on this tuff because I am still kinda stuck:

I am trying to solve the following problem. Assume we have defined the one-dimensional, normalized, Legendre basis as:

$$ \sqrt{\frac{2 \rho + 1}{b-a}} L_\rho (x)$$

Where $\rho$ is an integer denoting the degree of the polynomial.

This is an orthonormal basis defined on $[-1,1]$. Now, take the associated quadrature rule:

$$ \int^{1}_{-1} f(x) dx\approx \sum w_i f(x_i) $$

What we seek to do is project a function onto the basis, i.e. we seek to do:

$$\int^b_a f(x) L_\rho(x') dx$$

For each basis vector, where $x'$ is shifted and normalized from $[a, b]$ into $[-1,1]$.

[pine jetty]

Makogan

[next garden]

We further seek to extend this to multiple dimensions, in particular, three. So define:

$$ \int^{b_1}{a_1}\int^{b_2}{a_2}\int^{b_3}{a_3} f(x', y', z') L{\rho_1}(x)L_{\rho_2}(y)L_{\rho_3}(z) dxdydz $$

i.e. we have extended our one-dimensional basis into three dimensions by multiplying together, tensor style, the one dimensional polynomials. And we are now projecting our function onto that basis.

Now we want to come up with a quadrature rule for this setting. If we assume that our bounds are $[-1,1]^3$, this is not so hard, we can just do:

$$ \sum_{(i, j, k) \in [0, n]^3} w_i w_j w_k f(x_i, y_i, z_i) L_{\rho_1}(x_i)L_{\rho_2}(y_i)L_{\rho_3}(z_i)$$

Where $x_i, y_i, z_i$ are the points of quadrature.

I have code that implements the above and it works, tested it.

The part I am struggling to adapt is getting the correct values for when $[a,b]$ are not $[-1, 1]$.

What I thought would work is:

$$ \sum_{(i, j, k) \in [0, n]^3} w_i w_j w_k f(x_i', y_i', z_i') L_{\rho_1}(x_i)L_{\rho_2}(y_i)L_{\rho_3}(z_i) \frac{b_1 - a_1}{2} \frac{b_2 - a_2}{2} \frac{b_3 - a_3}{2}$$

i.e. grabbing the nodes $x_i, y_i, z_i$, shifting and scaling them to fit their respective bounds, then evaluating the Legendre polynomials on the standard nodes, the actual function on the shifted and scaled nodes, then re normalize by scaling down the result due to the change of bounds for the integrals. This is not working.

For example I tried to simply project the function $f(x, y, z) = 1$ and I got non zero coefficients for almost all projections. When, due to orthogonality, the expectation is that only the first one, corresponding to the constant basis function, should be non-zero.

[1.0000000000001648, 5.1961524227072, 5.196152422707172, 5.196152422707172, 46.95742752750177, 27.000000000001915, 27.00000000000192, 46.95742752750183, 27.00000000000194, 46.9574275275018, 484.1724899248873, 243.99795081108888, 243.99795081108883, 243.997950811089, 140.29611541308253, 243.99795081108905, 484.1724899248869, 243.99795081108888, 243.99795081108877, 484.1724899248873, 5355.000000000783, 2515.834056531001, 2515.8340565309986, 2205.0000000002524, 1267.8505432424427, 2205.000000000252, 2515.834056530996, 1267.8505432424427, 1267.8505432424429, 2515.8340565309945, 5355.000000000786, 2515.8340565309973, 2205.000000000254, 2515.834056530996, 0.0]

I am stuck, I was hoping someone could lend me a hand to understand what I did wrong.

[pine jetty]

Makogan

[vapid plume]

@next garden Sorry I'm still a bit confused. Your main task is to find $\int_{a}^{b} f(x) L_{p} (x) \dd{x}$ for each $L_{p} (x)$?

[pine jetty]

Sup?

[next garden]

the main task is doing that but in 3D

so L_p is really a polynomial in three vairables that is constructed as a tensor product of Legendre polynomials in one variable

[vapid plume]

@next garden Couple of questions just for my clarity lol.

The L_p here, is it the original L_p defined on [-1,1] or is it the one that you have shifted to [a,b]? If L_p is the original, then the integral is only defined if [a,b] is within [-1,1]. If this is the case, then you can easily transform from your [a,b] (which is within [-1,1]) to [-1,1] for the quadrature by applying the transformation to both f and L_p and changing the limits of the integral and also scaling dx, dy, dz, etc appropriately. If L_p is already shifted to [a,b], then the integral is defined so you can just apply quadrature directly.

[next garden]

I guess teh notation there is poorly done.

Assume on 1D our basis is actually:

here L_\rho is the stnadard [-1, 1] Legendre polynomial. x' is a normalized coordinate, where x \in [a,b] and x' in [-1, 1].

b and a are just the bounds of a one dimensional interval.

[next garden]

the one dimensional quadrature jsut requires the scaling factor from this comment

so you grab that expression I just shared and do quadrature on f(x) times it

the part I seem to have made a mistake is teh generalizaiton to 3D

[next garden]

Finally!!!

left experiment right ground truth

jesus that was an annoying bug to track

[wide folio]

Not even sure if this is the right place to ask this, but has anyone here used Fenicsx before ?

[jade surge]

yeah I tried to and honestly don't there are much simpler tools

[late kettle]

Simpler under the terms of what? The package is well conceived in the math community as far as I know. It hides the internals of the FEM implementation (for example) and focuses on the variational formulation itself

[grave spoke]

I was not particularly fond of the FEniCS documentation back then

Maybe they updated it by now, but before, there were so many things that were outdated

[buoyant fiber]

I have a theory with a large number of variables to be determined. I am using markov chain monte carlo. Is there a better method?

[summer hornet]

no

unless you discover one

[buoyant fiber]

Hmmmm okok.

Thx

[stone leaf]

what if he discovers one

[summer hornet]

fields medL

[steel shard]

could I have a hand with this exercise? I'm really struggling to come up with basis functions for this space

here's the piecewise linear case

but im not too sure how to extend this
i thought about maybe cooking up functions phi_j that were just -(x-x_{i-1})(x-x_{i+1}) and then rescaling the height but it doesnt really work

[vapid plume]

@steel shard I can recommend a book (there probably are many, I happened to read this one) which derives the basis functions.

If you want to do it yourself, good.

Work with local coordinates. We have an element which now needs three nodes and three functions. Each function is a quadratic function, so the original function, say V, is a linear combination of the three quadratic functions on that element. V^(e) (x) = V_1 phi_1 + V_2 phi_2 + V_3 phi_3, where the phis are our quadratic funtions and V_is are the values of the original function at the nodes. Obviously summing quadratics gives us a quadratic, so we have V^(e) (x) = a0 + a1 x^1 + a2 x^2.

Now we use the definition of the basis functions. We have that V^(e) (0) = a0, V^(e) (h/2) = a0 + a1 (h/2) + a2 (h/2)^2 and then carry on similarly for V^(e) (h). Solve this system of equations for a0, a1, a2. Substitute the a_i back into V^(e) (x) = a0 + a1 x^1 + a2 x^2 and then rearrange so that you get V^(e) (x) = (some function 1)V_1 + (some function 2)V_2 + (some function 3)V_3. Remember that V_1, V_2, V_3 = V^(e) (0), V^(e) (h/2), V^(e) (h). Then you get your basis functions from some function 1, some function 2, some function 3.

The Finite Element Method: Basic Concepts and Applications with MATLAB, MAPLE, and COMSOL, Darrell W. Pepper, Juan C. Heinrich.

[steel shard]

ohh ok very interesting, thank you!! that middle node was the key piece i was missing
that's very cool
thanks again!

[next garden]

Does anyone know if there is a theory for defining orthogonal polynomial basis functions for arbitrary convex polyhedra?

[wide spear]

No

[wide spear]

Seems tough, you might want to look into the finite element literature

[thin rover]

Anybody here particularly apt at numerical linear algebra? Just curious, I'm doing research on it and have yet to really run into anyone has done alot of it before

[fathom rain]

I do, but in a very niche subject

[thin rover]

what subject?

[fathom rain]

structured matrices mostly

[thin rover]

so wait is that like, tridiagonal, hessenburg and such?

Am I think of the right thing, google was unhelpful lol

[fathom rain]

for example yes

matrix sequences in general

[thin rover]

You ever do much with krylov subspace methods?

[fathom rain]

no

[thin rover]

rip lol

I am finding them interesting, but quite hard to actually code and to find information on

it also seems like alot of them just aren't done being researched

[fathom rain]

there is quite a bit of research in that direction

[fathom rain]

https://householder-symposium.github.io/schedule.html

or https://ilas2025.tw/files/ILAS2025-program.pdf

[thin rover]

thanks for the info, I'll look into it more, but it almost instantly went over my head lol

[fathom rain]

a lot of it for me too 🙂

[wide spear]

If you have specific questions just ask

[thin rover]

have you ever dealt with arnold or lanczos iteration?

[wide spear]

Yes

[open tusk]

I'm starting a numerical analysis course this fall, does anyone have advice for how to prepare or work through the course? Its being taught by a notoriously strict professor.

[hard flare]

I'd probably need more info about how the class is structured and what you'll actually need to do? Or do you mean in terms of material to review ahead of time

[open tusk]

I mean in terms of material to review ahead of time

I currently do not know how the class is going to be structured

[hard flare]

No problem. I'm assuming its a first course? Probably linear algebra and real analysis (if applicable)

If you happen to know what book you'll be using I could probably say something more specific but it'll generally boil down to those two things. And you wouldn't need to review anything terribly difficult from either of them, iirc

[open tusk]

The professor is a bit of an interesting person, and provides his own materials for his classes, so we don't have access to the materials until he gives them to us lol so we most likely won't have a textbook

[cunning raven]

what i safly discovered in practice is that for larger matrices the simpler krylov methods are too unstable and require far more iterations than supposedly required 😔

i suspect they are only really effective for very large sparse matrices

they are cool though

[wide spear]

Did you check the condition number of your matrices

[cunning raven]

they werent too bad iirc, they were just 8x8 dense matrices

other methods did show better results too

although they were not as lightweight as krylov 😔

[celest elk]

Lmk what you find

[wanton dawn]

Gemini excerpt: Key Names and Terminology
This is an active area of research in numerical analysis and approximation theory. If you want to read more, look for these terms and researchers:

Dubiner Polynomials: The foundational basis for triangles and tetrahedra.

Koornwinder Polynomials: A general class of polynomials on simple domains.

Warp and Blend Mappings: A specific technique for constructing the geometric maps.

Spectral/hp Element Method: The primary field of application. Key figures who have developed this theory include George Karniadakis, Spencer Sherwin, and Robert Kirby. The open-source code Nektar++ is a major implementation of these ideas.

[next garden]

Thank you enormously for these resources

II will try taking a look

[crisp birch]

Depends a lot on what you will be specifically doing during said course, but I can tell you that linear algebra was consistently pretty useful during my two Numerical Analysis courses.

Some calculus knowledge also helped.
If you also know on which software you will be coding, you could familiarise yourself with it in advance, so you don't need to lose time learning how it works.

[late kettle]

Does anyone feel satisfied with the Nitsche penalty method

[open tusk]

Thanks!

[next garden]

is the hp element method essentially a tree subdivision of your space with varying degrees of polynomials?

[wide spear]

Yes

[deep olive]

idk if this is the right channel for it, but does anyone know and understand how CORDIC works?

[wide spear]

What part do you not understand

[cinder topaz]

Rotation mod is fairly simple

You start with (1,0) and multiply it by a set of rotation matrixes
with each iteration getting closer to the desired angle

The angles are chosen in such a way, that multiplying by rotation matrix is very simple in binary representation

[deep olive]

interesting

[next garden]

Question. I have two functions, an approximating function, let;s call it $\phi$ and a ground truth function, let's call it $f$.

I want to know how much $\phi$ deviates from $f$, analytically that's just a simple integral of a difference $\int || f - \phi|| dx$

The problem is that $f$ is ridiculously expnseive to evaluate, so quadrature methods that attempt to give me high levels of accuracy take a lot of time. What altenratives do I have, if any?

[pine jetty]

Makogan

[fathom otter]

I can't think of any options really. Any measure of the deviation of f from \phi will involve evaluating f.

Depending on f and on the integrand norm(f-phi), there should be an optimal quadrature rule that will give you the best order of accuracy with the least evaluations of f.

How important is that you measure the difference with the integral and not with another measure, say |sup(f) -sup(phi)|?

[next garden]

not important. The use case is for an algorithm to decide if it should do more work or not. So really I just need a way of estiamting the deviation with as few calls to f as possible

[wide spear]

What do you know about the regularity of the function

Do you know anything about the derivatives

[next garden]

that it's an SDF, that;s about it

I can make little assumption about even smoothness in the strict sense

[wide spear]

I think your best bet is a trapezoid rule quadrature

So that at every iteration step you can reuse all the previous function evaluations

[fathom otter]

Is that not true for any quadrature rule? The quadrature points should be constant across iterations right?

[wide spear]

Not all quadrature rules reuse points when shrinking the grid spacing

[fathom otter]

So the grid spacing would be changing across iterations, I see.

Im not familiar with this problem

Why is it so expensive to evaluate f anyway. You just need to evaluate the signed distance of the quadrature points to the boundary, right?

[wide spear]

f can be expensive to evaluate for a number of reasons

This is not too uncommon

[fathom otter]

Hum

I see

Yours is a good suggestion

[cinder topaz]

Romberg might help too

[next garden]

well for example say your f is actually a high resolution mesh. For each evaluation of f you are doing a tree traversal, checking for bounding box inclusion and potentially even backtracking up the tree due to near misses. The worst asymptotic time for a BVH queary is O(n), and youa re doing that for each individual query. Not fantastic...

[next garden]

What is romberg?

[next garden]

I will try that out then, ty

Tangentially related but I realized a limitation of I think any quadrature rule?

Say my basis is an orthonormal polynomail basis,

Assume p_i is one element of my basis with i really big.

then for ALL finite basis with j < i, the projection of p_i onto te basis is 0 (due to orthogonality)

but that doesn't mean my method has converged

Anyway to try to estimate what I wouldbe without testing until I find the right index?

This is more personal curiosity I don;t think it will be a problem in practice

[fathom otter]

Right thanks for the clarification

[wide spear]

Also check out Richardson extrapolation

[wide spear]

Use a lower order quadrature like midpoint

[cinder topaz]

Romberg integration is Richardson extrapolation + trapezoid rule

[cinder topaz]

Now when I think about it

Doesn’t Richardson extrapolation rely on function having continuous derivatives?…
It might not work in your case

[wide spear]

The integral should have a continuous derivative

[tall solar]

Why is the best uniform approximation problem particularly interesting? For example, say I have a function f on a real interval J with real values and I want to minimize |f-p|_inf for some polynomial p where |.| denotes the L inf norm.

For context, I found an entire book about these problems that sure seemed to suggest this is a fundamental problem in approx theory but I can’t tell why this “norm” is particularly interesting.

[crisp lake]

I suppose I should mention adaptive Gauss-Kronrod-Patterson in passing. But what I myself am trying to do is characterising polynomials positive on the compact interval not representable as convex combinations of Bernstein polynomials of matching degrees. But I can’t even produce an example of such a polynomial.

[crisp lake]

It means the error can never exceed that norm in absolute value.

[wide spear]

Are you familiar with why people care about L inf error bounds

[tall solar]

No I’m not I mean to ask that

[crisp lake]

Should I be asking my Bernstein polynomial question in a different channel?

[crisp lake]

I tried $\int_0^1\left(\sum_{k=0}^n(p_k-c_k)\cdot B_{n,k}(x)\right)^2,dx$ and got closed forms for the integrals of $\int_0^1 B_{n,m}(x)\cdot B_{n,k}(x),dx $ but they didn't really give any insight into whether the Bernstein sum $\sum_{k=0}^n c_k\cdot B_{n,k}(x)$ with $c_0, c_n > 0, (\forall k) 1\leq k < n,,c_k \geq 0$ could successfully reproduce a degree $n$ polynomial $p(x)$ positive on $[0,1]$.

[pine jetty]

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

[wide spear]

I'm not quite sure what your question is

[wide spear]

Do you know why analysts care about L inf bounds

Why uniform convergence is important

[tall solar]

I have no idea why analysts care about L inf bounds rather than say, L2 or L1 bounds. I know the definition of these norms and that there is some duality between L inf and L1 but I have no ideal why analysts care about L inf in particular.

[wide spear]

You should think about this

Why is uniform convergence better than pointwise convergence

[tall solar]

I know that L inf is a stronger quality of convergence.

[wanton dawn]

I don't think they really care that much, but to have a pointwise error bound guarantee for discontinuous functions is very good. See https://en.m.wikipedia.org/wiki/Gibbs_phenomenon

[crisp lake]

Characterising polynomials positive on the unit interval not representable as convex combinations of Bernstein polynomials.

[tall solar]

Ah okay. I have been studying haar spaces and was hoping to see some applications of these. A Haar spaces is certain finite dimensional vector subspaces H of the real valued continuous function C(J,R) on a compact real interval J. Haar spaces are characterizing as exactly being the finite dim vector subspaces of C(J,R) which have a best uniform approx property. So for example any g in C(J,R) has a best uniform approx in span (cos(t),sin(t),cos(2t),sin(2t)).

Maybe it’s a niche topic but maybe it’s useful in signal processing for all I know.

[wide spear]

The convex combination restriction is pretty harsh

For example, consider f(x)=1

This is only expressible as a convex combination for deg 0 bernstein polynomials

f(x)=c for c>1 is not expressible as a convex combination of bernstein polynomials

[crisp lake]

Hmm. Given that what I'm hoping is for a non-vanishing denominator of a rational function, it may be that I bungled by not restricting the polynomials or modding out factors of the sums of absolute values of coefficients or similar. I think because the Bernstein basis spans the polynomials of the relevant degree, I think it may be a difference of $\sum_{k=0}^n\left|c_k\right|=1$ with and without the additional restrictions of $c_0,c_n>0$ and $(\forall k)0<k<n,,c_k\geq 0$. This is when considering a polynomial $\sum_{k=0}^n c_k\cdot B_{n,k}(x)$ where $B_{n,k}(x)=\binom{n}{k}\cdot x^k\cdot (1-x)^{n-k}$.

[pine jetty]

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

[crisp lake]

@wide spear So it seems to be coming down to whether $p(x)>0$ on $[0,1]$ is possible when some of the $c_k<0$ for $0<k<n$. Clearly $c_0=0$ or $c_n=0$ yields a zero at one of the endpoints of $[0,1]$, failing the strict inequality criterion and $c_0<0$ or $c_n<0$ yield negative values at the endpoints of $[0,1]$ which trips over intermediate values needing to go through a zero before becoming positive anywhere.

[pine jetty]

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

[wide spear]

Ok so you have some rational function? And you want the denominator as a Bernstein basis and want conditions on the coefficients for no singularities?

Or are you trying to approximate another function by a rational function?

[crisp lake]

Conditions on the Bernstein form of the denominator for no singularities within the interval, or at least to understand the way it fails to cover the whole space.

A sharp condition is always welcome.

Ultimately it feeds into trying to get sharper bounds on the error of a $C^d$ piecewise rational approximation on an interval than Gelfgren got in 1975 and 1978 by using real analysis vs. his complex methods.

[pine jetty]

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

[wide spear]

Well each bernstein basis polynomial is positive except potentially at the endpoints right

So I think that c_n > 0 for each n should be necessary and sufficient

[crisp lake]

I think the interior coefficients might be able to go negative.

[wide spear]

Hmmmm

I think there might be something in partition of unity theory

[crisp lake]

$\min_{x\in[0,1]}c_0\cdot(1-x)^n+c_n\cdot x^n$ is when $n\cdot(c_n\cdot x^{n-1}-c_0\cdot (1-x)^{n-1})=0$ so $\left(\frac{x}{1-x}\right)^{n-1}=\frac{c_0}{c_n}$ so $\frac{x}{1-x}=\left(\frac{c_0}{c_n}\right)^\frac{1}{n-1}$ so $x=1-\frac{1}{1+\left(\frac{c_0}{c_n}\right)^\frac{1}{n-1}}$ which with $c_0,c_n>0$ falls within $(0,1)$, at which point I think there may be room for negative $c_k$ for $0<k<n$.

[wide spear]

https://arxiv.org/abs/2309.10675

Have you seen this

[pine jetty]

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

[crisp lake]

Checking it out.

[crisp lake]

That's great.

[vocal shell]

Any references that you recommend for numerical PDE ?

[wide spear]

Any specific pde or general?

Iserles for general

[vocal shell]

Thanks you for that .
I think something that had quantum
mechanics in consideration would be great

[wide spear]

First become familiar with the various ODE time stepping methods

[vocal shell]

Aight

[crisp lake]

I think even numerical methods for SPDEs are based partly on ODE stepping methods.

[balmy merlin]

which book?

[tall solar]

Topics in approximation theory by Harold S Shapiro

[next garden]

I know Taubin's matrix is used to extract curvature (normal, principal directions, curvature values...) information from discrete spaces like poitn clouds and meshes. I remember I once read a resource that related all curvature information to the PCA of taubin's matrix.

So for example, if you wanted the direction of minimum curvature, it was the eigenvector associated with the x eigenvalue ordered by magnitude, things like that. Is anyone here familiar with that?

I can't find it anymore

[cinder topaz]

Not sure if that's the right channel to ask (channel description doesn't mention differential equations explicitly, but #odes-and-pdes doesn't mention numerical methods either)

I have a system of first order differential equations, some of them are wrt time, and others are wrt space

I tried a few ways to solve it numerically, and now I'm trying to frame it as a system of differential-algebraic equations
space derivatives are treated as algebraic equations ( 0 = (u[z+dz] - u[z])/dz - u'[z+dz] )
I use Rosenbrock method, and it works relatively well

but it seems there're some instability in region below 200, and with each step the error amplifies (picture below shows one of the components alone the time at different space points)

So the question is...
Is there a way to smooth the solution/any other ways to deal with that instability?

And if you find that approach is a dead-end, please, do tell me, how to handle these problems properly.
Thank you in advance.

[wide spear]

Have you tried using a centered difference

[cinder topaz]

No.
I'll report back with the results 🫡
Thank you!

[ruby fox]

hey guys i was working on some functional analysis , and found a paragraphe about the Bidual space E** of a space E , i can't understand why would we introduce such thing and what is it's importance in FA and EDP ( i wanna do a carrer in quant reaserch )

[glossy cipher]

its an important object of study in pdes and optimization theory, e.g. fenchel duality. reflexivity, which is the property that the canonical injection into the bidual is surjective, is also an important characterization of topological vector spaces (its involved in many central results e.g. kakutani, eberlein smulian)

[ruby fox]

thanks

[zenith jetty]

The bidual is essential in defining the completion of normed space in a beautiful and succinct way. Also without the bidual we wouldn’t be able to say that taking adjoints preserves the norm, meaning $|| T ||=|| T^* ||$

[pine jetty]

DazaiOsamu

[zenith jetty]

Actually the isometric embedding into the double dual is just so impressive whichever way you look at it

[ruby fox]

very interesting thank you so much

Monsieur

[brave crypt]

What is a good way to numerically find roots of polynomial of 5th degree since there is no formula for solutions in radicals

[wide spear]

Have you tried anything

Like bisection?

Or newton?

[hard flare]

yeah polynomials are friendly there are lots of great options!

[wide spear]

There’s a whole Wikipedia article about this

[jolly widget]

Hello,I have question about numerical methods.When my proffesor do two or three iteration manually he give some table with like 9 iterations.Do I have to do all that iterations manually or is there some method to calculate with table?

[wide spear]

You can use a computer

[hardy oak]

Kinda stuck trying to show that these nodal variables determine P_5

I've shown that if p vanishes under the nodal variables on the edges, that p = q * L_1^2 * L_2^2 * L_3^2

but q thus must have degree k - 6

Which I can't think of how to try to deal with those terms with the nodal variables in the interior

maybe using the largest line through the interior points but I can't seem to show it vanishes "enough"

[minor vigil]

guys, it is known that the matrix pseudoinverse can be used to solve least squares

for instance, $AA^\dagger b \equiv proj_{Im(A)}(b)$

[pine jetty]

jeca tatu

[minor vigil]

but is it possible to interpret geometrically what $A^\dagger b$ alone means?

[pine jetty]

jeca tatu

[wide spear]

What is the dimension of P_5

How many degrees of freedom do you have

How many data points do you have with that triangle

[hardy oak]

I meant to say P_7 sorry

But the issue is that we have one degree remaining which I’m not sure how to kill off the linear term

[hardy oak]

I just now realized why the psuedoinverse matters

It’s related to the classical projection defn

So like, let’s say we have a linear map between Hilbert spaces T: A -> B right

Then we have by projection that:
b = T(a) + r for some a in A right, with r being orthogonal to the image of T, i.e T(x) • r = 0 for each x in A

However that is equivalent to x • T^t(r) = 0 for each x in A, thus T^t(r) = 0.

So T^t b = T^t T a

If we can INVERT T^t T, then we have

(T^t T)^-1 T^t b = a, so
T (T^t T)^-1 T^t b = T(a), we’ve projected b onto the image of T via T (T^t T)^-1 T^t

The Penrose psuedoinverse is (T^t T)^-1 T^t, the inverse of T mod the orthogonal complement

In general we can just say T^+(b) = a where b = T(a) + r

Actually wait does this mean that the Penrose psuedoinverse exists for any continuous linear map into a Hilbert space H

No the image needs to be closed, shit

[lean barn]

hey what is the formula used here?

or how to derive this result

[glossy cipher]

thats neat

ive always just thought of it as the svd thing

[hardy oak]

Me too lmfao

I am going to try to keep myself on a regimen

At least 4 exercises every 2 days

[glossy cipher]

good pace, im gonna try to do the same

[late kettle]

It's all consequences of the minimization problem which is usually proven by orthogonal projection. You want to minimize some rejection vectors (that's why they're called "least" squares) and generally SVD and SVE in Hilbert space come guarenteeing a regularity condition for both linear systems and operators (the resulting pseudoinverse rely on this framework)

[lean barn]

pls guys i'm going insane

it's some divided difference shenanigans

[tidal slate]

Looks like a barycentric Lagrange interpolation.

I'm not at a computer to write TeX right now, but if you still need help beyond that tip I'll be able to later.

[minor vigil]

Hi

Anyone can explain why in Inverse power method with shifts one considers a vector x0 out of R(B - lambda I) and not out of R(B - 1/(lambda - alfa))?

[late kettle]

Repeated operations project the vector towards the eigenvector associated with the eigenvalue closest to the chosen shift. The primary requirement is that the initial vector has a non-zero component in the direction of the target eigenvector. And therefore, the method converges in this way as an advantage.

It mainly care about the properties of the shifted matrix rather than initial vector you begin with

[amber briar]

What's a good reference for a proof of Gaussian quadrature

I know nothing about numerical methods so idk any textbooks

[minor vigil]

someone can explain this?

I dont see the relationship between QR method and inverse power method

[late kettle]

As usually being said a complete discussion of Gaussian quadrature requires a substantial background in Approximation theory. For beginners the book by James F. Epperson (the proof avoids a lot of the need of approximation theory and is based on an idea due to Roy Mathias) and the classic one by Stoer & Bulirsch. For further treatment in the theory I'd recommend Narayan Kovvali

[lean barn]

What's stopping me to do this

and this would be 0 since it's for n+2 points and the polynomial x^a has degree <= n

[lean barn]

nvm I am completly wrong, found some other formula

[tidal slate]

Out of curiosity, what did you find?

[lean barn]

I've scroled through another book of my proffesor's and there was an exercise to prove the equality between the devided difference of any arbitrary funciton f(t)/(x-t)

which equals L(x0,x1,...,xn,f)(x)/(x-x0)...(x-xn)

tried the proof before again and it was pretty trivial

But now I once again am stuck on another proof

i was thinking generalized Rolle?