numerical-analysis

[lean barn]

but I need a good auxiliary function that = 0 in the points x0,x1,x2,...,xn

and simplifies nicely when taken the nth derivative

[lean barn]

nvm it's the Lagrange mean value t

th*

[obsidian leaf]

In the 2 phases method how do you solve the table?I know simplex algorihm.Is 2 phases the same?

[late kettle]

Well, it's an alternative of the Big M method. You utilize the simplex method through the multiple phases in the absence of some readily feasible solution. If you start with a BFS solution that has some constraints not satisfied, then introduce a large enough artificial variable modifying the original constraints x0 that arbitrarily satisfies the constraints. At the end you want the variable x0 to be equal to 0 to make it a leaving variable and thus non-basic, which allows you to have a feasible dictionary equivalent to your original LP problem.

This is all happening in phase I where you force all artificial variables to be zero. Solving the first phase will yield to sometimes varying cases. We're mostly intrested in this case where either minimizing the sum of all artificial variables or maximize(-x0) which is our objective function (not the original objective function in LP) that we're intrested in and particularly when its optimal value is equal to zero. This yields the phase II where we apply simplex method and the resultant optimal solution in this phase is the optimal solution (at least one) to our original LP

[mellow linden]

yo guys will need help with my numerical analysis end of semester exams

we are to solve analytically and code the answers using python

if you can be of help kindly send me a text but will also post the questions here

oh dont worry there will be a reward too

for your deed

[wide spear]

Unfortunately we cannot help you with exams

[mellow linden]

oh thats very unfortunate. oh im not trying to cheat tho. its a mock exams should have been clear tho.

[wide spear]

Ah ok

If you send a question, we can try to help you reason it through

[lavish scroll]

hi would anybody tell me if my process and/or answer is incorrect? im in numerical analysis but havent had experience with most of the topics before

let me know if i can be any help expanding on what i did

[wide spear]

Looks reasonable

[lavish scroll]

thank you

[ember nexus]

Do we have a single general formula for the order of finite difference for a given order of derivative? I searched through the internet for some PDFs and similar resources, but to no avail.

[wide spear]

I'm not quite sure what you mean

But

If you have function values at a set of points at {a1, a2, a3, ...., an}*delta + x and you want to numerically approximate the derivative of order n, then you can find the finite difference coefficients by solving a matrix equation

[ember nexus]

You can approximate the n-th derivative with any order. I meant this.

[ember nexus]

Do you have any resources for this?

[wide spear]

https://web.media.mit.edu/~crtaylor/calculator.html

Read the "how does it work" section

[ember nexus]

Thanks a lot!

[rich galleon]

My white whale recently has been trying to take a crack at fast, precise algorithms for calculating the Bessel functions and their zeroes in constant time relative to the size of the input, the order of the function, and, for the zeroes, the order of the zero

I’ve made only so much progress (I can use one of the integral forms to calculate a reimann sum for J_a for a not negative which seems to be precise to almost 20 or so floating point bits and runs in a few tens to hundreds of microseconds)

But are there people here who are also interested in this?

[wide spear]

If you come up with something nice then there will be people who are interested in it

[rich galleon]

This is what I’ve come up with so far. If you are generous, you might call it semi decent (and despite what I claimed at the time, it does seem to be the integral that is suffering in precision)

[wide spear]

What precision are you targeting

[rich galleon]

I would like to hit 64 bit float precision

This is a pipe dream without a lot of work

[wide spear]

If quadrature accuracy is holding you back, you might check out curtis clenshaw quadrature

[rich galleon]

But there’s a reason I’m calling a white whale

[wide spear]

Have you looked at the NIST recommended algorithms

[rich galleon]

The recommended algorithm from numerical methods (which I think gets copied over into NIST, I’ll have to check) requires a recursion

[wide spear]

Are you targeting 10^-16 relative error or absolute error

[rich galleon]

Ideally, relative error

See above on “pipe dream”, especially near the zeros

[wide spear]

Actually I don't understand how you want constant time relative to the size of the input

As z increases, cos(zt) becoems more oscillatory and the integral requires more points to accurately resolve it

Unless you can use fourier transform tricks somehow

[rich galleon]

I’m cheating and only calculating cos(zt) one time for the smallest non zero value

[wide spear]

Oh I see is that the chebyshev recursion you mention

[rich galleon]

For the others, I can get basically exact values with the chebyshev recursion

Yes

[fathom rain]

just so you know, double precision is typically referred to as 53 bit precsion.

[wide spear]

[rich galleon]

I’m less familiar with the specifics of the nomenclature than I would like, so thank you

[wide spear]

I would disagree with that

[rich galleon]

I mean, it’s an argument about whether you’re referring to all the bits in the computer or just the bits that provide precision

[fathom rain]

the significand precision is what defines the precision in a reasonable sense

[wide spear]

At any rate it doesn't really matter for the problem

Hmmm ok so you need more accurate quadrature

[rich galleon]

Anyway, the biggest disadvantage of the chebyshev recursion approach is that it constrains me to a uniform spacing for the quadrature

[wide spear]

But you don't want to move away from uniformly spaced points for the cosine recursion

Simpson's rule instead of trapezoid?

[rich galleon]

Yeah, you hit the nail on the head

[wide spear]

Have you tried using simpson's rule

[rich galleon]

For the whole interval?

Or do you mean cutting the interval into pieces that I use Simpson’s rule on?

[wide spear]

Well you're using subintervals for the trapezoid rule right

[rich galleon]

Yeah

[wide spear]

So you can use the same points for simpson's, just with different weights

[rich galleon]

Hm…

I might give that a whirl

I’ll give it a try when I get to my computer proper, but it might not work as well as we might hope

Cosine is fairly oscillatory

But it’s a better idea than what I’ve had

[fathom rain]

why not use high precision data types if you want something accurate?

[wide spear]

Quad precision won't help if the integral error is O(2^-20)

Runtime of hundreds of microseconds is already kind of pushing it for a single function evaluation

[fathom rain]

that seems very slow (dunno what language and benchmarking is used)

using BigFloat (notoriously slow) and SpecialFunctions.jl it takes about 2-3 microseconds

this is with 256 bit precision

[rich galleon]

For the Bessel functions?

[rich galleon]

I’m using Rust: partially so that I can do weird pointer shit if need be, and partially because I would like the tooling to exist for rust

[fathom rain]

I ran besselj maybe you do something else?

https://specialfunctions.juliamath.org/v0.1/

[rich galleon]

Was it consistent for nu and z?

[fathom rain]

pick any you like and i try

[rich galleon]

Can you try all the combinations of nu from 0 to 10 by 0.1 and both 10 to 100 and 10.1 to 100.1, each by 1, and then z from 0.0 to 100.0 by 2.0, please?

[wide spear]

It does look like the julia code is calling C code in the background

And also requires nu to be between 0 and 1

[rich galleon]

I would have been much more surprised if it wasn’t, tbh

[fathom rain]

there is also Bessels.jl in pure julia

[wide spear]

Unless I'm misreading

[rich galleon]

Like I’ve implied, the Bessel functions are a notorious pain to calculate

[fathom rain]

100 100 took 43 microseconds

is this the solution you want for that 0.09636667329586155967431402487040184831175541982502185591797352259414559484434024

10,10 benchmarks to 5.500 μs (I use BenchMarkTools)

[wide spear]

Last 2 digits differ from mathematica eval

But still much more than 10^-16 accuracy

[rich galleon]

That’s basically irrelevant though

[rich galleon]

Yeah

[fathom rain]

well all floating point computations will have noise in kast digits

[rich galleon]

What algorithm is Julia using?

And is it open source?

[fathom rain]

https://github.com/JuliaMath/SpecialFunctions.jl/tree/master/src

[rich galleon]

I might try porting it if it’s an MIT license or some such

[fathom rain]

line 641 in bessel.jl is the call for different data types

uses mpfr i guess bessel is built in (like Angetenar said)

[fathom rain]

https://www.mpfr.org/algorithms.pdf

[rich galleon]

Thank you!

[fathom rain]

seems there are wrappers for mpfr in rust

[rich galleon]

I’ll check them out, but I also want to implement this in native rust

[fathom rain]

You can also check out https://github.com/flintlib/flint/ might be faster (often is)

arb was merged into flint

[rich galleon]

I’m trying to avoid using arbitrary precision arithmetic if possible

[fathom rain]

then i would probably look at this https://github.com/JuliaMath/Bessels.jl/tree/master/src/BesselFunctions

[cursive heron]

Anyone have any book recommendations for numerical analysis? In particular I'm interested in learning how the convergence of Lagrange interpolation of a function f is related to the poles of f when extended to the complex plane

[wide spear]

Surely trefethen has some papers about this

[cursive heron]

https://people.maths.ox.ac.uk/trefethen/papers.html that's a lot of papers 💀 most of those papers are way too advanced for me, but I found Numerical Analysis by Sauer, which might contain what I want

or maybe Approximation Theory by Trefethen

[wide spear]

Trefethen writes in a very approachable way I'd say

[clever siren]

Hey there,
I'm a CS major and I thought of taking a course in numerical analysis from the math department, does anyone know if it will help with or related to CS, or is it overkill?

(I think my domain will probably be in computer vision and graphics)

here is its description (it was translated):
"Introduction: Various uses of numerical computation. Representation of numbers on a computer, numerical stability and computational errors. Methods for solving equations: the intersection method, the point of intersection method, Newton's method, geometric methods. Order of convergence of a method. Interpolation by polynomials, Lagrange, Newton and Hermite methods. Error blocking. Spline functions and Bernstein polynomials. Numerical differentiation and integration. Numerical algebra: Finding eigenvalues, solving a system of equations using the Jacobi method and the Gauss-Seidel method."

[wide spear]

Hmmmm

It'll have some applications

But

If you can take a course on numerical linear algebra, I think that would be a lot more applicable

[clever siren]

i can either choose from calculus 3 or numerical analysis
or i can take formal verification or algorithmic game theory from the CS department but they are a bit too far from my field

[wide spear]

I think multivariable calc is more important

If there is a fourier analysis class available, I would recommend taking that in the future

[clever siren]

really? i think thought that numerical analysis is more related to CS and computer vision

or im just too scared to take calc 3 lol

[wide spear]

I think that the content of an intro numerical analysis is also generally pretty straightforward to learn on your own

[clever siren]

ok tysm ❤️

[late kettle]

By the way If you want a good introductory book to learn it on your own I highly recommend the one by James F. Epperson

[brave crypt]

guys how would you solve this question? Ik you need to get the formula of Sn = n/2[2a1 + (n-1)d] but how do you do ahead of that?

[wide spear]

#precalculus

[brave crypt]

ohk sorry

[molten knot]

splines will probably help with graphics and optimization stuff for cv

[wide folio]

Hi

I'm running a simulation rn, and I'm having a lot of trouble generating convergence plots, especially temporal convergence.

Do any of you have an suggestion/recommendation on a book or something that I can refer?

[wide spear]

What sort of simulation

[wide folio]

It's a FEM + SDC implementation of the convection - diffusion of a gaussian pulse

[wide spear]

SDC = spectral deferred correction?

What happens when you try to do a convergence plot

Does not converged at desired order?

[wide folio]

Yes

[wide folio]

Yea 🙁

it's a pretty large deviation 😔

[wide spear]

Both spatially and temporally?

[wide folio]

spatial convergence seems ok

Temporal convergence is where I'm running into problems

[wide spear]

Have you tried another time stepping method

[wide folio]

Well, I've tried implicit Euler

No problems there

[wide spear]

Does it converge correctly without convection

[wide folio]

Uhhh so it's like this

I'm also using the semi lagrangian method, to take care of the convection term.

I'm not super sure if it's relevant here

But I do have a code for just diffusion, but I haven't run any convergence tests on that

[wide spear]

Oh goodness

There is so much going on

[wide folio]

Yea 😔

[wide spear]

Ok with a semi-lagrangian advection, as you decrease the time step, you remesh more, which increases the amount of interpolation error you are accumulating which can adversely affect the time convergence

So I would test it without the advection

If there are still issues, then you SDC has problems

[wide folio]

Oh

Thanks a lot for the inputs though

I'll try running it without advection and see if that changes anything

[brave crypt]

is this the channel for Num Lin Alg ?

[wide spear]

Yes

[cursive heron]

pi_n(x) goes to 0 as n -> inf for both equispaced nodes and Chebyshev nodes, right? So it's sufficient to prove that M_n is asymptotically smaller than n! in order to prove that interpolation converges, correct?

[wide spear]

Does pi_n go to 0 for equispaced nodes?

[cursive heron]

oh wait, I need some assumptions on the range of x 🤔 it definitely goes to 0 for x in [-1, 1]. And for x in [a, b] it's atleast bounded, so the convergence is still determined by the asymptotics of M_n vs n! I think?

[wide spear]

Are you sure it's bounded

[cursive heron]

huh, I think I just understood the utility of chebyshev nodes a bit better by thinking about this. If you want to prevent pi_n(x) from converging to 0, you would pick an x near the endpoints of the interval [a, b], let's say x near b, because then there could be a term (x - x_0) with x_0 near a, forcing pi_n(x) to be big. So with chebyshev nodes you increase the density of nodes near the endpoints so that you get lots of terms (x - x_i) close to 0. Hmm, it made more sense in my head than written down tbh 😅

[cursive heron]

hmm, not really... so when you partition the interval [a, b] into n equispaced subintervals, for a given x you increase the number of nodes close to x, so you get more factors (x - x_i) close to 0, but you also increase the number of nodes far away from x, so you also get factors (x - x_i) that are much bigger than 0. I'm not sure whether these effects cancel each other out, let me think

[wide spear]

For equidistant nodes, || on [-1,1], pi_n is bounded by n!*h^(n+1) where h=2/n ||

[cursive heron]

I see I meant it should go to 0 for x in [0, 1], or generally for intervals with length <= 1. But I see now how pi_n grows pretty fast on larger intervals

I just did some experiments in python, and pi_20(x) on [0, 10] attains a maximum of about 111397184318

thanks for the help btw, you're pretty much the only person keeping this channel alive

[wide spear]

Also you should think about how this applies to a function like 1/(1+25x^2) on [-1,1]

[late kettle]

I think of the utilization guaranteed by the Chebyshev zeros through the use of residues or complex remainder. If you define a family of curves C(p) that are smooth concentric curves about the midpoint of [a, b] you reduce by using Chebyshev nodes by making contours become true ellipses having foci at some x value (say positive negative 5) and that should make it easier to choose a contour which completely encloses the interval without hitting the singularities. A discussion on that could be summarized by the theorem in this included image.

You could also avoid that complex-typed justification but it yields only partial answers. Both justifications can be found inside the notes named On the Runge Example by Epperson.

You can also find some useful treatment on this file starting with the function Ange mentioned

[cursive heron]

Interesting, thanks! is it correct that using Chebyshev points give the smallest ||f(x) - p_n(x)||_inf for any function f? ie. choosing chebyshev nodes always give the best interpolation?

[wide spear]

No, the chebyshev points are not optimal

Trefethen has a paper about this

[cursive heron]

okay, I just saw this in the lecture notes above, but I think I'm misunderstanding what they mean

[wide spear]

This is a common misconception

[cursive heron]

so that statement is just wrong?

they prove that 2^{-n} * T_(n+1)(x) is the monic polynomial of degree n + 1 with the smallest infinity norm, but then they jump to this conclusion without any elaboration

[late kettle]

They are near-optimal for wide classes of smooth functions, but they do not guarantee that for every possible f. Just a very good generic choice but not universally optimal

The file still useful though I didn't mention it particularly for this conclusion. The Epperson notes is more elaborate and kind of authentic material to follow from of course

[cursive heron]

I see, thanks

[wide spear]

The Remez algorithm is probably the most popular method for constructing L_inf optimal interpolants

Closely related is the search for optimal Lebesgue constants for various sets of nodes

In the sense of Lebesgue constant, the Chebyshev nodes are asymptotically optimal but can be slightly shifted to improve them

[sturdy spire]

https://www.youtube.com/watch?v=8E2858kWVEw&list=PLAVG7GMBpcYArR9QLXm3DVvqYhRdF6Tsj

https://www.youtube.com/results?search_query=numerical+linear+algebra+series

[cursive heron]

If A is a matrix then presumably ||Ax|| is a vector norm. Then the argmax is the vector x with ||x|| = 1 that maximizes ||Ax||

I think it can be any norm, I don't think it's common to assume the 2-norm is used when not specified

[cursive heron]

Quick question: I'm interpolating a function for which I know interpolation converges. The spikes in error for large n here are likely due to floating point errors, right?

[wide spear]

Please plot this on a log-log plot

[cursive heron]

okay, give me a sec

[wide spear]

Do those seem like floating point errors to you

[cursive heron]

hmm... not sure what floating point errors look like tbh

[wide spear]

Floating point errors should be O(10^{-16}*number of operations)

I suspect that there are issues with the interpolation

[cursive heron]

so if the number of operations increase linearly with n, then the floating point errors should also increase linearly

I can send the interpolation code if you're interested, it's a bit ugly though

[wide spear]

I do not want to see your code

I don't even want to look at my own code

[cursive heron]

lol

tbh, it's not at all clear to me why it's more likely to be a problem with the interpolation rather than floating point error

[wide spear]

Those are too large for floating point errors

You are not doing 10^16 floating point operations with degree 60 interpolation

For interpolation a lot of the time the floating point error also does not grow too large

How are you doing the interpolation

[cursive heron]

just with numpy, close your eyes if you don't want to see:

# computes the k-th Lagrange basis function, and evaluates at x
def lagrange_basis_function(k, xs, x):
    return math.prod( (x - xs[i])/(xs[k] - xs[i]) for i in range(len(xs)) if k != i )

# computes the Lagrange interpolation polynomial for points (xs, ys) and evals at x
def lagrange_interpolate(xs, ys, x):
    return sum( lagrange_basis_function(k, xs, x)*ys[k] for k in range(len(xs)) )

[wide spear]

I suspect there is a cancellation issue where x is close to xs

Where are the Barycentric Lagrange weights?

Oh

You aren’t using the Barycentric form

[cursive heron]

Nowhere, because we haven't learned about that 🙈

[wide spear]

The non-Barycentric form of Lagrange interpolation is prone to numerical instability

[cursive heron]

I see

[wide spear]

Does the basis correctly evaluate to 1 when x=xs[i]

[cursive heron]

yep, even with n = 100

I gotta admit something, but you have to promise not to get mad
-# I used equispaced nodes

with chebyshev nodes there's almost no increase in error for large n

[wide spear]

Well well well

[cursive heron]

I mean, this is just an assignment, the goal is just to learn stuff. And I learned that equispaced nodes suck

but thanks for your help 🙏 I'll look into barycentric interpolation when I need to interpolate something important

[sturdy spire]

I have a question about this diagram

I don't see which one gives sparse solutions (it's supposed to be L1)

Since the circle (L2) and the square diamond (L1) both hit the x and y axes at some points, wouldn't you argue they both contain sparse solutions?

[wide spear]

Ok

L^1 norm does give sparse solutions and L^2 norm does not

But

Horrible picture for this

https://stats.stackexchange.com/questions/45643/why-l1-norm-for-sparse-models

[sturdy spire]

reading it rn

[sturdy spire]

Ok I understand the article, which is basically saying: "L2 mostly penalizes larger entry and not the smaller ones, whereas L1 treats them equally".

But how does that relate to this diagram (we can use this better one)?

[wide spear]

Is this diagram better? I'm not sure

https://miro.medium.com/v2/resize:fit:1600/format:webp/1*_e8BLNA749W_7yxi7hz-DA.gif

[sturdy spire]

Ok there's one with a convex function

looks like an ellipse that rotates and shifts randomly

And it just so happens that most of those rotations leave the vertex as the constrained min

[wide spear]

So the ellipses are level sets of a distance function from the target point

[sturdy spire]

yes

[wide spear]

For most rotations/variants of the distance function and most locations of the target point, one of the L^1 vertices is the closest point

[sturdy spire]

That looks geometrically intuitive

[wide spear]

Another way to think about this is if you pick some region of the plane and randomly pick points, on the circle, the closest point will be uniformly distributed

Whereas for the square, the closest point will more commonly be one of the corner points

[sturdy spire]

Thanks for that, it helps me understand

How would you get that distribution

How often will the square's vertex be the solution?

[wide spear]

Does this make sense

You can try to compute the distribution yourself

[sturdy spire]

How did you choose the circle's radius

Or it doesn't matter does it

[wide spear]

In the limit as R->infinity it shouldn't matter

[sturdy spire]

As R -> infinity the square becomes a little point and the bands become negligible

so 0 right

or wait nvm, that's the constraint not the loss

The radius is the distribution of possible loss functions

Wow, this way of understanding loss is a lot harder (I thought it would be easier)

https://tenor.com/view/facepalm-gif-9004202298607570582

[late kettle]

Are you referring to the idea that setting a very large regularization parameter (approaching infinity), the penalty term will completely dominate the objective function. So you start thinking the only way to minimize this massive penalty is to set all coefficients to zero? Because it's not clear to me what you're thinking of exactly here.

The image Ange posted was indeed clear (emphasizing symmetric loss functions) and in the more general case, loss functions can be asymmetric and at an angle which results in more zeros for L1 (both are shown in these following comparable images)

You can generate random loss functions with scaling and compute the distribution to check that. Otherwise what's the issue?

[cursive heron]

Can someone help me interpret what exactly I'm supposed to do here? Should I create 10 different plots, one for each n in {1, ..., 10}, with each plotting k against the error rate?

[cursive heron]

I can't see any other way to interpret it, but it's just weird to have to make 10 different plots... not sure if I should try to fit all 10 graphs in the same plot, or if I should just have 10 plots, one after another

[sturdy spire]

Ok I want to simulate this. What would be a good objective function? |Ax-B|?

Wouldn't it behave differently for every problem (so I might as well focus on my favorite type...)

actually it wouldn't be that hard or inefficient to program several different objectives (logistic/normal regression, proximity, classification any kind of convex objective) and see the distribution for myself once I've coded the main simulation in 2d

[sturdy spire]

would it be approach to start by trying various |Ax-b| setups with different b vectors and see the result (should hopefully verify the first diagram)?

[wide spear]

You should be able to do all 10 on the same plot

[late kettle]

I can give you a simple hint but in the inverse form. Plot something like this log-log plot but inverse the idea of h -> 0 and make it as a function of the number of discritzation as stated

[late kettle]

You can use this expression to generate random loss functions (And you can obviously combine a logistic loss for example with the regularization in other realistic cases)

For this expression a~U(0, 10) should scale the bowl in the beta_1 direction, and b~U(0, 10) the same for the beta_2 direction, c~U(-2, 2) should control the amount of tilt/angle that's away from vertical or horizontal, and (c_1~U(-10, 10), c_2~U(-10, 10)) is the position in coefficient space of the minimum loss function value. U(j, k) means uniform random variable between j and k.

You should better use soft constraints as it would be more intuitive for the implementation

[cursive heron]

Thanks for the hint I ended up with this; it's a bit ugly but I think it's okay

[wide spear]

It isn't an error rate

A rate is a change with respect to something

It's just the error

[cursive heron]

oh, good point, thanks

[wide spear]

This reasoning makes sense right

[tall solar]

I know a univalent polynomial p(x) in Pn of degree n is uniquely determined by n distinct points.

Is it also true that a polynomial is uniquely determined by specifying n many initial conditions?

For example, is there a unique poly f in P5 such that f(0)=0, f’(0)=1, f’’(0)=2, f(1)=3, and f(2)=4?

[wide spear]

You need one more condition

Polynomials of degree 5 have six coefficients to specify

[tall solar]

Thanks do you happen to know why 6 conditions uniquely specify a polynomial?

I think to show 6 distinct points specify a polynomial uses the vandermonte determinant because a unique solution to the interpolation problem corresponds to a Linear system whose determinant is the vandermonte determinant.

Idk how I would prove this similar result. I think it’s a determinantal condition. This also feels like a result about odes

[wide spear]

You can just set up the linear system for the coefficients

The matrix has full rank as long as you aren't specifying multiple conditions for identical function values

You can search for results about confluent vandermonde matrices if you want

[late kettle]

You can establish a simple connective justification through the Rank–nullity theorem, but as Ange said such things will follow immediately from well-posed conditions and linear independence.

[late kettle]

I haven't worked extensively with quadrilateralized spherical cubes, but I can see the logic. Have you tested it with a posteriori data points?

[wide spear]

Someone else confirmed it for me

[late kettle]

Good!

[wide spear]

It simplifies to 6n^2d^2+2

[sturdy spire]

Yeah, about that... the paper I'm interested in uses hard constraints, and so does the diagram above. Thankfully, I should be able to jack up the constraint coefficient to effectively infinity to remedy that (what do you think, Lagrange Multiplier or Indicator function)?

[sturdy spire]

(this diagram)

I appreciate the formula, will lyk if I make use of it

[late kettle]

Lagrange multipliers will suit this type of setting within the previously mentioned formula

[bronze gazelle]

Hi, is there any way to get mathematica for free?

[wide spear]

Through your university

[fathom rain]

You can install it for free in jupyter notebooks/just use the scripting language.

[bronze gazelle]

Yes I just did this! Tysm

[tall solar]

I think your right that the exact thing I’m looking to understand is: why the confluent vandermonte determinant is nonzero on any point of the closure of the open simplex defined by a<t0<…<tn<b ?
thanks

[quiet wren]

@final pine

[wide spear]

[stoic thicket]

porting question #advanced-lounge message

hello everyone

might be a long shot, but does anyone know any papers which discuss using multiple species in graphical CFD? I am working on a project and I was wondering if anyone has described discretization methods for admixtures of gasses.

[wide spear]

Are the gasses very different? How large are the concentrations?

Are there sources and sinks for the various species?

[stoic thicket]

I intend to use some simple gasses, but also might implement some fun ones with outlandish properties. at basic, its meant to simulate air, CO2, hydrogen.

there are not intended to be static sources and sinks but created ones, like a gas canister or an air siphon.

I was thinking about also maybe handling a fluid sim in the same system, just with some added parameters for gravity etc, but I think I should keep the system simple in principle for now before I worry about that.

I'd like to use a eulerized system with voxel gridding for densities.

[wide spear]

You can just keep a mass percentage for each gas you want to track

[stoic thicket]

correction- yes, sinks and sources.

[stoic thicket]

could I just ease between these percentages during the diffusion step? idk if that would result in normalish looking behavior.

[wide spear]

If they are well mixed then you just apply a Laplacian diffusion to each component

[stoic thicket]

hm

[wide spear]

But for gasses, viscosity is near 0 anyways so you can omit viscosity

[stoic thicket]

yeah the environment will not be working with gas densities where viscocity would matter most of the time

for liquids that's a different story

[wide spear]

Water is as inviscid as air

[stoic thicket]

for water like liquids yeah I probably don't need to care either

and I could maybe cheat for viscous things

doesn't need to be perfectly accurate to real life, just convincing.

the only thing that I think might complicate the system is explosive decompression which I would like to attempt to model, but I may need to cheat there as well.

[wide spear]

That should be fine as long as you are careful with your numerical discretization and time stepping

[stoic thicket]

Now that I think about it I don't think I need to do anything too fancy for liquids... its essentially just an incompressible gas.

Just needs a gravity parameter.

Blegh, documenting this in code is gonna be a nightmare. I might need to just pseudocode it up in an .md

[wide spear]

[stoic thicket]

Think a gauss-siedel would work just fine?

[wide spear]

Why are you solving a linear system

[stoic thicket]

all of this is new to me so I'm going off of the papers I've already read on the subject, from what I've looked into when you're trying to work in real time trying to solve densities backwards as a system of linear equations is leagues faster computationally than just marching forward

[wide spear]

It depends

If you share the papers you are looking at I can give more detailed explanations

[stoic thicket]

sure, that would be wonderful

one moment.

mainly jumping off from these papers.

http://graphics.cs.cmu.edu/nsp/course/15-464/Spring11/papers/StamFluidforGames.pdf

https://cg.informatik.uni-freiburg.de/intern/seminar/gridFluids_fluid-EulerParticle.pdf

https://www.cs.cmu.edu/~kmcrane/Projects/GPUFluid/paper.pdf

[late kettle]

Oh Stams papers are for graphics purposes I don't think you're required to do that much of physical accuracy

[wide spear]

Are you trying to create a solver to embed in a video game?

[stoic thicket]

indeed

[wide spear]

Ok

How accurate do you need the model to be?

[stoic thicket]

Only needs to be accurate up to a point where a player with basic engineering knowledge would understand what acting on these gasses would entail, and is convincing to the eye.

in terms of physical spaces, I believe cells will be relatively large. somewhere on the order of a quarter meter maybe, not sure yet.

[wide spear]

How large is the whole domain

[stoic thicket]

I'd like a level will encompass an area somewhere around~ lets call it 200m^3.

[late kettle]

If that's for testing the qualification for entering the general area of Numerical analysis research, then you're expected to be able to supply any requested details, such as proofs and derivations, etc. Check the topics if you have an examination syllabus!

There are numerous facts and techniques that should be known by every PhD graduate (for example) who claims a research interest in numerical analysis. A 90-minute examination is inadequate to test all of these, but a reasonable subset is usually to focus on the material typically covered in the courses that are forming a basis for your examination syllabus.

For each topic (concept, theorem, problem, algorithm, etc.) listed or covered, standard questions that might be asked include its definition, existence, uniqueness, characterization, derivation, proof, applicability, sensitivity, stability, accuracy, convergence, computational complexity, etc., as may be relevant.

So these are well expected.

[stoic thicket]

I intend to design the game such that two levels can coexist and interact, but no more than two.

200 is probably on the conservative side but I'll see how much I can bite off.

[wide spear]

I think maybe lattice Boltzmann would suit your requirements

[stoic thicket]

I'll look into it, on the surface it looks like it is modelling the same thing that paper no.2 mentions where vectors are discretized along cardinal directions

marker and cell

[late kettle]

Then they're typically short, focused vivae testing your understanding of key concepts, algorithms, error estimates and simple proofs rather than original research.

You can still be asked to justify or prove a norm inequality you used, especially if it underpins an error bound, but questions stay shallower and more applied than phd qualifying orals.

You don't need that deep theoretical generality. But as I understand from your provided image above on your grading scales. Understanding some theory is required for Pass with distinction which is above the standard pass grade.

So, spend some time on the theory if you want to rock it

[stoic thicket]

thank you for the suggestion, this looks much closer to what I'm going for

[late kettle]

Feel free to ask any questions regarding that -- if it looks vague to you.

[stoic thicket]

if I were to make levels connect at arbitrary locations (like "portals"), would it make sense to just have a small additional dimension for the domain and just set a boundary condition that may have holes poked through it?

I don't know if that makes sense.

maybe there's a smarter way to do it.

[wide spear]

What do you mean by levels

[stoic thicket]

two separate 3d volumes

[wide spear]

How are they connected?

[stoic thicket]

portals interface both volumes in a noneuclidian fashion

at arbitrary locations

might just have to be a separate process

[wide spear]

But like

What is happening physically

[stoic thicket]

uh like in the game?

two "universes" connecting together

they're overlayed physically in 3d space

[wide spear]

You can just connect the grid points then can't you

Or "identify" two of them or however many you need

[stoic thicket]

yeah, I think it can just be treated as another spatial dimension

in principle I don't think it should cause anything to blow up?

since they're closed systems

aside from gravity maybe

but that's fun

could be clamped

[stoic thicket]

Oh, wait I see what you mean

like just directly tying them?

hm

that's probably way simpler than what I was imagining

[wide spear]

Just treating it as a single point

[stoic thicket]

yeah OK that's very doable then

[mystic wadi]

What is the fastest way to solve transcendental and linear equations

[late kettle]

For nonlinear equations the general approach is Newton's method but that will be drawn differently according to the type of convergence you want (local or global convergence) and efficiency ofc. Such cases you'd use different variations of the previous method or generally another method (might include another variation) for robustness for such different situations.

The same idea goes for linear equations -- like is it a dense problem or a large sparse system? For both ways you need to add some additional property you want to achieve and be more precise.

[stoic thicket]

is there a parameter you can look at to decide whether a system is relatively dense or sparse?

[wide spear]

Number of nonzero entries in a matrix

[stoic thicket]

ah

[late kettle]

It's also (to be fair) worth to mention in this case that no universally agreed-upon exact percentage of zeros that qualifies a matrix as sparse exists for a particular efficient computational operation. It's all heuristic, so you have to get a sense of when it becomes computationally beneficial to use special data structures and algorithms that only store and process the non-zero elements inside your sparse entities.

Be a little bit wise of the needs of your algorithms and storage. The threshold for using sparse property is a little bit context dependent for the different choices of algorithms, or operations, etc.

And asside from this idea of "Sparse measure" -- just stick with the mathematical strict idea of the less number of nonzero entries and don't confuse yourself.

[mystic wadi]

Thanks

[cursive heron]

Is the Lagrange interpolation polynomial the best polynomial approximation if you pick the optimal set of nodes to interpolate on? I mean, if you select nodes x_0, ..., x_n that minimizes the 2-norm ||f - p_n|| where p_n is the Lagrange interpolation on x_i, is p_n the polynomial that minimizes ||f - p_n|| among all degree n polynomials?

[wide spear]

There is a single polynomial that interpolates f the best

(barycentric) Lagrange interpolation is one way of constructing a polynomial that is close to optimal using pointwise data at a set of interpolation points

Another set of polynomials, the Bernstein polynomials, have theoretically desirable properties

Like

You can use the bernstein polynomials to give a constructive proof of weierstrass approximation theorem

However, they aren't a nodal basis and are less attractive to compute with

[cursive heron]

I see There's a question on my assignment that asks what the relation is between the best approximation polynomial of degree n, and the Lagrange interpolation polynomial on nodes that minimize the 2-norm. I have no idea what to answer - like, they're not the same, and you can't use one to compute the other. Any ideas what they might be looking for?

They also ask for a way to compute the optimal interpolation nodes without using gradient descent. We have covered orthogonal polynomials, so it might be related to that. Based on how the question is written, I feel like they're implying that the Lagrange interpolation polynomial actually is the best polynomial approximation

[wide spear]

[cursive heron]

I asked the TA, and he just said the best polynomial approximation isn't necessarily a Lagrange polynomial, but we're trying to approximate it by Lagrange interpolation on optimal nodes. I still have no idea what to write, so I'm just gonna write something vague and generic

[wide spear]

Lol

[cursive heron]

Another tricky question they asked was to discuss the computational cost of Lagrange interpolation vs splines. I have no idea where to begin - do they want me to talk about asymptotic performance? Or the performance of evaluating an n degree polynomial vs indexing into a list of k-degree polynomials and evaluating it? I feel like it depends so much on the specifics that it's impossible to give a useful generic answer

[wide spear]

Yes, asymptotics as a function of the number of nodes

[cursive heron]

Hmm... so evaluating an n-degree polynomial requires O(n) multiplications (or maybe more? You can calculate $x^n$ in O(log(n)), so an n-degree polynomial would be on the order of $\sum_{k=i}^n log(k)$ I think, not sure if that is superlinear). But using splines, if you have n nodes divided into k subintervals, then it would be O(n/k)

[pine jetty]

sheddow

[late kettle]

I think it's simpler than that so we avoid such misconceptions.

Splines can't evaluate in O(n/k) unless you mean something (Horner's method)-related which makes part of total the evaluation O(1) by interpreting n, k for something equal.

Splines then in the typical "standard manner" evaluate by O(log(n)), and you can interpret that as a binary search time for interval finding. Having constant time won't change that.

Also I don't see a problem in providing an asymptotic total complexity and not just evaluation. Because it won't make a notable difference between the two methods as splines would still be better.

So you can say (assuming Barycentric for precomputing) the total effort for Lagrange is O(n^2) + O(Nn) under the same order. Highlighting which is which in that order

[late kettle]

You asked about this too.

Orthogonal polynomial roots would imply the use of the Golub-Welsch Algorithm as far as I can telI. I don't usually think of gauss nodes exactly in this way but maybe you're targeting something L2 related optimality

[pine jetty]

mq

[cursive heron]

I think I get it by keeping k proportional to n, we get that evaluation on each subinterval is constant, but looking up the subinterval is log(n)? What is N in O(Nn)? I'm guessing evaluating an n-degree polynomial is O(n^2)?

[cursive heron]

I just looked up Horner's method, seems I was a bit optimistic in thinking polynomials could be evaluated in O(n)

[late kettle]

Well, this is wrong.

The total evaluation is O(log(n) so adding O(1) won't change that. Your interpretation assumes O(1) is itself the total effort for "splines" which isn't.

You can also think of the spline evaluation as two steps each with an effort.

N typically denotes the number of evaluations so evaluate once means O(n) -- this is now general I'm not talking about splines at this point just clarifying what N means.

O(n^2) denotes the weighted precomputaion done once in the Barycentric form

[wide spear]

Well

It all depends on how the polynomial is written

Surely if you are doing barycentric lagrange then you are also using chebyshev points so there is no precomputation needed

Usually barycentric Lagrange interpolation is considered O(n)

[cursive heron]

Your interpretation assumes O(1) is itself the total effort
What do you mean? I'm saying that evaluating a polynomial of constant degree takes O(1). If you have a list of polynomials of degree c, then evaluating one of them takes O(c^2), which is constant time if c doesn't change

[cursive heron]

oops, I misread Horner's method btw, I thought it said O(n^2), but it uses O(n) multiplications and additions

[wide spear]

When using horner's method you have to be careful about numerical stability issues

[late kettle]

@cursive heron This is again if you're not doing precomputaion. I just wanted to clarify that in case of the standard definition you can write the precomputaion part (if it's there) along the subsequent evaluation O(n) if you wanted

[fathom otter]

Hey guys.

Does anybody know of any books or articles that discuss solving a least-squares optimization problem when the Jacobian matrix is both ill-conditioned and rank-deficient? I'd like to know how these properties affect the optimization algorithms such as Gauss-Newton or Levenberg-Marquadt.

What I know until now is that when the Jacobian is rank-defficient, then there are several solutions to the minimization problem. In that case, which one do the algorithms naturally choose?

Regarding the condition number of the Jacobian, how does that affect the convergence of the algorithms, and are there any algorithms best suited for this kind of situation?

[wide spear]

Use a SVD

[fathom otter]

My problem is large-scale. Is that a good way of approaching it?

[wide spear]

Well

It's one way

Alternatives being what

Sobbing and crying

[fathom otter]

I am using Levenberg-Marquadt's method for now, and it works fine, even if the Jacobian is ill-conditioned and rank-defficient. The advantage being that it works for large-scale problems as well. But I wanted to better understand to what solution is the method converging and how the condition number affects the convergence.

[wide spear]

For levenberg-marquadt the choice of damping parameter should help with ill-conditioning and rank-defficiency

[fathom otter]

Is there an optimal way of choosing that parameter based on the ill-conditioning?

[wide spear]

The damping parameter lambda should be chosen dynamically based on the statement of the optimization

If the update is smaller than desired, increase lambda

[fathom otter]

Do you have any references in mind for that suggestion? That'd one way of improving what I have already.

What about other types of regularization, other than Levenberg Marquadt's choice? Where can I learn more about them? Do you know any good book on numerical least squares problems that discusses this topics?

[late kettle]

Maybe Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion by Per Christian Hansen for such ill-posed problems

That should give some insight even if you're tackling nonlinear equations

[late kettle]

You can follow something like this. It's heuristic in practice

The minimization problem is then close to a standard linear least-squares problem

You'd just increase it for stronger regularization or as Ange told you. The rest of of the logic follows

[fathom otter]

Thank you very much, I'll give that book a look.

[fluid cliff]

Hi, I'm currently investigating some stuff linked to the smoothness and stability of 2 dimensional polynomials using a tensor degree structure i.e.

P(x,y) = p(x)*p(y)

this introduces higher order cross terms, say i use p of order k

id now get terms of the form $x^i y^j$ such that $i+j \geq k$

I'm searching for some references about the influence of these higher order cross terms on the total polynomial, as I'm doing some specific discontinuous galerkin FEM related stabilization research for which this information would be of great importance, only i can't find anything, as basically all mathematicians just use a total degree polynomial such that they remove all $i+j \geq k$ modes...

also, i do not have a clueeee as to what channel this question belongs to so feel free to tell me to gtfo and kindly provide me another more suited channel 🙂

[pine jetty]

thijs2725

[wide spear]

Can you quantify what sort of information you want?

For the higher order FEM (spectral elements) I am familiar with, people do not ignore the higher order terms

But people also don't write their polynomials in the monomial basis

If you're working on a rectangle, which I assume you are because of the tensor product polynomials, then working with all the terms is the natural thing to do

Because you also probably are working with a tensor product grid

If you are working on triangles, you have to be more careful and ignoring the higher order terms is the correct thing to do

This is the so called "triangular truncation"

[fluid cliff]

Yes, the thing is, im working on shock capturing methodologies, for which we indicate troubled elements of the DGFEM solution, and to which we apply a limiter, this limiter im researching is applied to a hierarchy of polynomials within the total polynomial which in our case is a tensor product form

so i see 2 possible solutions, 1: i define this hierarchy as just the subsets of lower degree tensor product spaces, or 2: i define the hierarchy as the total degree spaces within the tensor product space. so where option 1 limits basically the entire layer of polynomials such that either i or j is equal to degree k, the option 2 limits the cross terms first... either of these of course works perfectly fine as at the end of the day its just limiting to ensure stability, the idea is that it is desirable to not limit more then necessary, so im looking for information surrounding the introduction of instabilities and if these will first introduce themselves in the highest order cross term ie x^k y^k or if this will first manifest in just the entire layer x^i y ^j where either i or j =k...

do you follow??

and ofcourse my hypothesis is that the cross term will first introduce the instabilities as this is the highest total degree, however my proof for this is basically well uhhh i feel like thats the case, and thatsss why i need references

[wide spear]

Why are you writing the polynomials in a monomial basis

[fluid cliff]

because thats a lot easier then writing orthonormal legendre polynomials

[wide spear]

No but you want to be using an orthonormal basis right

[fluid cliff]

yes

[wide spear]

Do you know what points you are working with in each rectangle

[fluid cliff]

how do you mean what points?

[wide spear]

The best thing in my opinion would be to write the polynomials in terms of a nodal basis

[fluid cliff]

modal basis is one of they key aspects of indicating instabilities though as they introduce themselves in highest order modes first

which is kindof the most important part

[wide spear]

What

Nodal not modal

[fluid cliff]

no? modal? as in the orthogonal basis functions?

[late kettle]

Hesthaven & Warburton "Nodal discontinuous galerkin methods"

Might be useful

You can combine both

With some projection I guess

[wide spear]

A nodal basis is one like the Lagrange basis

[stable trellis]

hi this is a somewhat stupid question but i'm trying to do a nonlinear least squares regression and the way i've written this is as follows:

import jax.numpy as jnp
import jax, jaxlib

def nls(fxns, guess, tol = 1e-4):
    """
    Gauss-Newton nonlinear least squares solver.

    Args:
        fxns:   List of residual functions, y - f, where f has parameter dependence and y
                the numerical evaluation. We want to minimize residuals squared.
        guess:  Input parameters that we then iteratively update.
        tol:    Tolerance for convergence of parameters.
    
    Returns:
        Optimal parameters that minimize the sum of residuals squared.
    """

    def vec_fxn(params):
        """Vector-valued function containing all of the parameters."""
        return jnp.array([f(params) for f in fxns])
    
    jacobian_matrix = jax.jacfwd(vec_fxn)

    error = 10.0
    beta = guess

    while error > tol:
        residuals = vec_fxn(beta)
        J = jacobian_matrix(beta)

        JTr = J.T @ residuals
        JTJ = J.T @ J

        JTJ_reg = JTJ + 1e-8 * jnp.eye(JTJ.shape[0]) # regularizer

        delta = jnp.linalg.solve(JTJ_reg, JTr)
        beta_s1 = beta - delta

        error = jnp.linalg.norm(beta_s1 - beta)
        print("Current iteration error: ", error)
        beta = beta_s1
        print("Current parameters: ", beta)
    
    return beta

what often ends up happening though is that i get stuck in a loop where the error doesn't decrease, but i don't understand why

However, when I change the JTJ_reg to instead include

JTJ_reg = JTJ + 1e-6 * jnp.diag(JTJ)

this now converges, and i don't really understand like, why

[fluid cliff]

Hi, can anyone here point me in some directions to look for computing lower degree polynomial projections of high degree polynomials?

i.e. if i have a 5th degree polynomial, i need to find the best fit solutions for 0th 1st 2nd 3rd and 4th order polynomials. The most straightforward way (and sadly the only way i know) is an L^2 projection.. only this becomes EXPENSIVE in certain cases due to curved spaces, high dimensionality, and high orders... therefore if anyone has some other ideas on how to attain lower degree polynomials at a lower cost id very much like to hear those

[wide spear]

Hand compute some updates where it gets stuck

[wide spear]

Write your polynomials using an orthonormal basis

[fluid cliff]

curved space, so this does not work as elegantly

[wide spear]

What do you mean by curved space

Like

A sphere?

You just need to be careful about your choice of basis

[fluid cliff]

no mate, when using orthonormal bases, on a regular interval -1 to 1 these bases would be orthonormal to each other which results in a diagonal mass matrix which lets one just remove higher order degrees to attain a L2 projection into a lower degree polynomial space

In CURVED space however, this is NOT possible, one can only define an orthonormal basis on the reference element, only the polynomials are not defined in reference space but in phyiscal space, and due to the curved nature of the domain, this mapping is NONCONSTANT which results in a non diagonal mass matrix, meaning the truncating of higher order modes is no longer valid as the coefficients would be wrongly scaled

[wide spear]

Can you define what you mean by curved space

[fluid cliff]

i can NOT define an orthonormal basis on (a), maybe i COULD, but it would not be feasible for millions of elements at the same time and therefore it is never done

therefore the basisis orthonormal on (b) and a nonlinear jacobian maps between these

[wide spear]

Ok

What are your computational requirements?

How were you computing the L2 projection?

[fluid cliff]

compute a mass matrix for the lower degree polynomial which is a square matrix, compute a mass matrix for the original polynomial, which is rectangular, and then compute the solution for q of:

M_q q = M_u u

in this case i made the assumption that the mass matrix is a precomputable constant matrix... of the size nDOF(lower degree) x nDOF(lower degree) and nDOF(lower degree) x nDOF(original degree)

however i also have certain computational cases in which the mass matrix is NOT a consant and has to be recomputed every single interation... which makes inverting it every iteration pretty much unfeasible... so conjugate gradient methods or something would have to be used... but thats also sheit and way too expensive in contrast to just truncating modes as is possible for rectilinear elements

[wide spear]

Can you write the mass matrix change as a rank-n update?

[fluid cliff]

nooo, veeeeryyy nonlinear change

[wide spear]

What original DOF do you use?

[halcyon thunder]

Do you get nicer behavior choosing points near chebyshev points? Or do you have to be exactly on them?

[late kettle]

Deviating from the exact roots certainly opens a slight room for reintroducing the less near-optimality behavior especially with bad regularity and increasing degrees. You have to be precise about how nice you want the interpolation

[wide spear]

I think redpoplar was asking thijs2725

Maybe

Unclear

I don't know

[late kettle]

[wide spear]

[fluid cliff]

Is variable but a DOF of around 125 seems to be the norm im aiming for

[wide spear]

5x5x5 in each cube/deformed cube?

[fluid cliff]

Yes

[wide spear]

Gauss Legendre points?

[fluid cliff]

For numerical integration yes

[wide spear]

Why don't you want to use say 3x3x3 in each cell?

Or lower order

[fluid cliff]

Because more = better and im more focussed on the mathematical procedure itself then a precise number, but the target audience indeed often chooses anything between 3x3x3 and 5x5x5

Lower than that it is not accurate enough and higher than that is often too expensive

[wide spear]

Well for lower order you balance it with smaller cells

[fluid cliff]

Well that is also more expensive

Key is to find a optimum

[wide spear]

Hardy please see #latex-help

[shy talon]

alr

thanks

[fluid cliff]

But for now the issue is the L^2 projection is pretty damn expensive however it is meant to be used to go to lower orders anyway, therefore it might not matter that much to have a slightly worse approximation

So if there are other methods which can relatively cheaply approximate a k’th order polynomial using a lower order polynomial that would be worth trying out

[wide spear]

The cells are not changing right

Can you precompute an orthonormal basis in each cell using G-S and then reuse that at all later iterations

[fluid cliff]

The cells are not, however the mass matrixes are, due to uhmm well a lot of math which is used to attain stability, making the mass matrixes dependent on both the basis functions as well as the solution itself

[wide spear]

The mass matrices are not what

They are changing right

[fluid cliff]

Oh yes brain fair yes the mass matrices are fhanging

[vestal roost]

Hi, i have received a file with x and f(x) values, but i dont know how to import the file into my program. Does anyone know how to do that? I am doing a problem that wants me to find the distance a rocket travels within a certain time, given the graph of the velocity function (but not the function itself - i only have the x and y values).

[vestal roost]

this question might be more related to computing-software?

[late kettle]

I assume the file is csv or whitespace-delimited with time and v columns?

[vestal roost]

this is what the file looks like:

the row on the left hand side would be the time, while the row on the right hand side would be the velocity

[late kettle]

read it with pandas.read_csv or numpy.loadtxt into arrays, sort by time if needed

then approximate distance = numpy.trapz(v, t) (or use cumulative trapezoid via numpy.cumsum of 0.5*(v[:-1]+v[1:])*diff(t)) to get displacement

[vestal roost]

so i need to write this in my code before i import the file?

[wide spear]

#computing-software

[wide spear]

Yes you were correct

[late kettle]

At this point yes :+)

I'll continue with you in #computing-software

[radiant salmon]

Spps. $\alpha >1$, construct a sequence ${p_n}_{n=0}^{\infty}$ that converges to 0 with order $\alpha$

Can someone help me w/ this?
I understand what the question wants but i just can’t seem to think of an example, if someone could lead me in the right direction that would be great

[pine jetty]

KySquared

[radiant salmon]

Wait nvm i think i found one
Looks like i was overthinking this

[fathom otter]

Hey. Suppose I have a linear inverse problem:

y = J x + n,

Where J is constant, but ill-conditioned and rank defficient and y are given measurements with noise n.

In my case, x are the constant changes in conductivity values (change compared to some background) on the elements of a FEM mesh, and the values themselves are not so important, but the image they form is.

What are the conventional ways of correcting for the effect of noise in distorting the image formed by x?

I was looking a bit into JP Kaipio's "Statiscal and computational inverse for problems", and if x and n are drawn gaussian random variables, then y is also one, and Kaipio gives expressions for the posterior mean and covariance.

Apparently, for the case of white noise and white prior, the posterior mean is the solution of the classical tykhonov regularized problem for y = Jx. I looked into some ways of using the poaterior covariance to correct for noise, but none seem to visually improve the reconstructed image quality, so I was wondering if you know of amything that'd possibly work.

[wanton dawn]

what is the "image formed by x"? Is part of the problem deciding what that is?

[wanton dawn]

gemini gave a good answer with several things to try: 1. Total Variation (TV) Regularization. 2. Sparsity-Promoting (L1) Regularization. Advanced: Hierarchical Bayesian Models, Learned Regularization (Deep Learning).

[fathom otter]

Hey @wanton dawn , thanks for the reply.

The aim of the inverse problem is to image the inside of a domain, for example, image the lungs in the chest cavity. In this case, it uses the fact that when you breed the conductivity in your lungs changes with the air inside them. For this application, it's not so important to actually measure the conductivity of the lungs, but to observe their change is. Reconstructing the conductivity change, in this case, would get you an image of the lungs working. It's not part of the problem deciding what the reconstructed image is. Hope I clarified what I mean by the image formed by the conductivity. The inverse problem i'm talking about is electrical impedance tomography, even though I'm working on a variant of that.

I think, of your suggestions, the simplest thing to try is going to be TV regularization. Indeed, different regularizations seem to be more or less resistant to noise, I tried Tykhonov and another one.

However, I was looking for suggestions of an a posteriori correction to the reconstructed x. In this PhD thesis I was looking into, they compute the standard deviation of each element of x ( I think they assume the distribution of x to be gaussian, but it is very unclear how they compute the std), and then employ this:

x_new = x/std(x)

As a simple correction, but I am not able to reproduce their results, and was wondering if anybody is aware of this a standard technique to correct for noise?

[teal eagle]

Hello great numerical analysis fellows, sorry for spamming. I'd would have some help regarding a problem I have in numerical analysis. the question could be found here: #help-47 message (help 47).
Thank you for your time and consideration.

[indigo cape]

Hello, does anybody here know how to use the method of momenta on antenna design?

[alpine mirage]

never worked with antenna design but if you mean method of weighted residuals then it's usually just galerkin method
so the main challenge is in choosing the basis functions and the algorithm is just solving a linear equation

[crisp lake]

Hmm. Is there a recurrence or semi-closed form to generate two-point Lagrange-Hermite interpolation basis polynomials with the endpoint knots of arbitrary but equal orders? Hopefully also with some sort of way to represent constraints to be of lesser degree.

[crisp lake]

I found https://iiif.library.cmu.edu/file/Traub_box00027_fld00028_bdl0001_doc0001/Traub_box00027_fld00028_bdl0001_doc0001.pdf but I can't entirely make sense of it.

[tall solar]

Hey sorry I don’t know numerical analysis but I skimmed that paper. So are you trying to find a recurrence relation just like there is supposedly a recurrence relation for the so called bell polynomials?

I don’t quite understand the exact polynomials you are looking to describe. I know it’s the case hermite interpolation works if you specify enough points. I’d be interested to read more about these recurrences your interested in

[woven haven]

anyone here has experience working numerically with the Kadomtsev–Petviashvili equation? or good literature references?

[wide spear]

What problem are you trying to solve

Periodic domain? Free space? Boundary values?

Do you have a method you want to use?

How much accuracy do you need

[woven haven]

I was thinking about something like integrating factor pseudo spectral for space variables in the x direction, not sure about the y direction. Not sure how to deal with the nonlocal term too. Time-stepping something explicit like Runge Kutta 4th would be preffered if possible. Homogeneous neumann boundary conditions on a rectangle. Not sure about accuracy

so basically the strategy I'd use for the KdV plus something smart for the x integral of u_yy

[wide spear]

Ok it sounds like you have a rough idea of what you want

Concerns I have

Ummmm

Well ok

Have you implemented this

[woven haven]

I have implemented the KdV with pseudospectral scheme (+ integrating factor). Its a 1+1 dimension problem, the method works fine against the analytic traveling solitons, very simple. I want to see if I can get a similar scheme for the KP which is 1+2 dimensions. Essentially its an initial boundary value problem, I'd like to plug in the know line-soliton solutions and see if they evolve correctly, just like the KdV. I want a finite difference based method for space and explicit time-stepping scheme for evolution in time.

one can think of the KP as the KdV + a term which is:

$\partial_x^{-1} u_{yy}$

[pine jetty]

Sydd

[woven haven]

I dont know how to treat this at all

[wide spear]

Well, you have a psuedospectral method which to me means that you chopped up your domain and have a local medium to high order polynomial basis in each element

Can you differentiate and integrate these basis functions?

[woven haven]

yes, more specifically I am using a Fourier method so the basis functions are just the exponentials with different frequencies

so I guess I can diff and integrate by just multiplying and dividing by these frequencies, respectively

I think I'll try just the spectral basis on the x directions and use a simple second order finite difference operator for the double y derivative

the inverse derivative here translate to multiplying by -i/k in fourier space, I'll try this out with fft and see how it goes

so actually its gonna be periodic in x

and neumann in y

the boundary conditions

[wide spear]

If it's a fourier basis then why is it pseudospectral

[woven haven]

this is a bit of semantics, but okay, let's call it a fourier method then

[wide spear]

At any rate, differentiating and integrating in the fourier basis should be pretty straightforward I think

[crisp lake]

I was trying to specialise Traub's more general equation for n points of order p to the case n=2 & get a computationally useful simplification for that case. What I've found is that Bell polynomials have properties enabling their reductions to Stirling numbers of the first kind, which then reduce to a falling factorial, which then cancel a factorial & it's simplifying relatively well.
The basic idea is to form a sort of basis for polynomials of a higher degree, say, 2p-1 via the derivatives at interval endpoints up to the p-th order & impose constraints so they're of lower degree. When the smoke clears, I basically want to solve for C^d piecewise rationals & cook up a basis that makes it significantly easier.

[cursive heron]

I'm struggling to understand how to compute a multistep method in practice. On the left side here you get a sum of x_{n-4}, ..., x_{n+4} and the right side also relies on x_{n-1}, ..., x_{n+1}. Am I supposed to solve this linear system for each iteration? Or do I take an average or something?

[wide spear]

You are given x_{n-4} through x_{n+3} and find x_{n+4}

[cursive heron]

Ah yeah, that's pretty simple dunno why I was confused, I think maybe because wikipedia mentions implicit and explicit methods wrt to multistep methods

[paper hollow]

is this where i can come to ask questions about fixed point iterations?

[wide spear]

Sure

[crisp lake]

Once I get the basis I think will simplify the calculations worked out, I think it's going to end up as a QCQP w/o many if any useful properties. It may be that changing the basis yet again could help get some helpful properties. I think Farouki-Rajan Bernstein bases could be a plausible fallback if the Lagrange-Hermite affairs fall flat. It's a long road to producing improvements on the error bounds & runtimes given in Gelfgren (1975 & 1978).

[crisp lake]

I think objectives could be arc length or variance, which are still pretty tricky. Maybe tanh or tanh o log of the obvious objectives could compensate for their blowing up to infinity when denominator roots fall within the interval. Someone somewhere has likely commented on ways to compensate for the singularities.

[wide spear]

If you have a singularity of degree d at a point c, you can multiply the function f by (x-c)^d which removes the singularity

Then work with (x-c)^d*f

If you get an approximant for (x-c)^d*f, then divide through by (x-c)^d

[crisp lake]

It wouldn't be obvious that a rational function's denominator has a root in the interval & the singularities would move from the POV of that.

[crisp lake]

It's a C^d piecewise rational interpolant, so the variables are the coefficients of the rational functions within the subintervals.

[crisp lake]

I might want to learn more about when nonsingular/bounded C^d piecewise rational interpolants can & can't exist. Involving multiple intervals b/c the derivatives aren't part of the problem makes it complex.

[cursive heron]

We're supposed to solve this differential equation using various methods up to time T = 10^5 with a step size of h = 0.1. But even solve_ivp with RK45 from scipy takes almost 10 minutes. Am I doing something wrong? Or have they set Tmax way too high? If every method we implement takes as long as solve_ivp, this notebook will likely take more than an hour to execute

I guess it's hard to say if I'm doing something wrong if I haven't shared my code, but my question is rather, is it common/expected that it should take this long?

[wide spear]

Let's see

10^5 end time with a time step of 0.1 means you're doing 10^6 time steps

So you're doing 10^5 time steps per minute

I don't think this is slow

If you want something faster then just write C

[cursive heron]

Yeah, that's what I feared, thanks for confirming

[minor vigil]

Standard Gradient descent with constant step size can be cast as an ODE

In this new setting, it is possisble to prove with stable manifold theorem that it avoids strict saddle points almost surely

I wonder if its possible to extend this result to the case where we have armijo line search

What would be a natural ODE in this case ?

[wide spear]

Instead of thinking of it as a different ODE, maybe you can think of it as an ODE discretized with an adaptive time stepping

[pine jetty]

a

a

[left silo]

ok

so gradient descent with a constant step size corresponds to sampling time with a fixed delta

and armijo would lead to different delta_k's

depending on the previous x(t_0+k*delta_k-1)

but, doing this you are no more analyzing the iterative method as a dynamical system, right? but as a discrete update rule

[wide spear]

I don't know what jeca tatu is trying to do

[vagrant whale]

need help trying to troubleshoot my MATLAB code: it's attempting to estimate the largest eigenvalue of the given matrix A using inverse iteration to an accuracy of 10^-6 but keeps converging to the wrong one

x = [1;0;0;0];
A=[[1,2,3,4];[4,5,6,7];[2,1,5,0];[4,2,1,0]];
lambdas = [0,0];
a = 15;
B=A-a*eye;
[L,U]=lu(B);
for i = 1:2
    y = U\(L\x);
    x=y/norm(y);
    lambdas(i)=x'*A*x;
end
while abs(lambdas(i)-lambdas(i-1))>10^(-6)
    y = U\(L\x);
    x=y/norm(y);
    lambdas = [lambdas,x'*A*x];
    i = i + 1;
end

should go to 11.1055 but instead terminates at -0.0397

[vagrant whale]

oh i may be a dumbass lmao

[vagrant whale]

huh im not even sure what im doing wrong it still converges incorrectly

[fathom rain]

set a=0 and you will get an eigenvalue

[vagrant whale]

it gives me the smallest magnitude one as expected in this case yea

i should've made it clearer that the hw wants me to find the largest one with a being the smaller of the 1 or infty norms of A

and x being e1

full problem statement (holy yapping)

[vagrant jackal]

since you didn't specify a size of eye, it defaulted to 1 (scalar). this didn't raise an error since matlab assumes you want to subtract the matrix with all entries of that scalar

[vagrant whale]

op

😭

lemme see if that fixes it

wow this is what my dumb ass gets for not reading the documentation closely

AH

that fixes it

thank you so much 😭🙏

[idle warren]

Which university program is this for?

[wide spear]

For future reference debugging questions should go in #computing-software

[cursive heron]

I'm struggling to understand how to implement multistep methods in a way that reduces roundoff errors. I'm currently reading Solving Ordinary Differential Equations, and it starts off easy enough, but when it comes to the (numerically accurate) implementation, there's 3 different matrices, and you need to split the characteristic polynomial for some reason (how?) and there's a shift operator etc... Is there an easier explanation somewhere? Or maybe an actual implementation I can look at?

[sage river]

Hey guys, a question in signal processing: given an open, irregular path drawn on the surface of a sphere (which can be represented in spherical coordinates; imagine drawing a freehand line on a ball), I want to perform FFT and iFFT on it to achieve low-pass filtering and uniform amplitude scaling. What is the appropriate way to do this?

[wide spear]

Ok

Let's see

So open path = not a loop right?

Usually spherical coordinates are a bit unpleasant to work with because the coordinates themselves have a singularity

One thing you can do is consider the path in (x,y,z) coordinates and take the FFT/iFFT component wise

Parameterizing it as (x(t), y(t), z(t))

It depends on what you want to do with the FFT though

[sage river]

Do you mean perferoming fft to x(t), y(t), z(t) respectively? My concern is that, since I plan to do low-pass filtering (remove high fequecnies from the original spectrum), and use ifft to reconstruct the signal, (new_x(t), new_y(t), new_z(t)) may not stay on the surface of the original sphere.

[wide spear]

The easy thing to do would just be to project the new_x,new_y, new_z to the surface of the sphere

If you want to avoid this then you can use spherical coordinates

Apply a rotation so that the curve avoids the poles

If you want something that takes into account the surface of the sphere is slightly widen the curve so that it's a 2d patch on the sphere, then take the spherical harmonic transform of the indicator function of this 2d area

[sage river]

I really appreciate your help.
For your spherical coordinates suggestion, I don't quite understand how "apply a rotation" is performed;
For your spherical harmonic suggestion: this is suppopsed to be an unit sphere and I only want to process this curve without considering the surface of the sphere. Given this situation, would you still recommend doing it this way? If so, could you provide a sample or a reference link for how to handle it?

[wide spear]

For rotation: https://en.wikipedia.org/wiki/Euler_angles

What do you mean you don't want to consider the surface of the sphere

[sage river]

Sorry I may have stated it incorrectly. In fact, I am analyzing the motion trajectory of the endpoint of a limb segment (for example, the forearm). The starting joint (such as the elbow) is fixed and the limb length does not change, therefore this can be simplified as analyzing the trajectory of a point moving on a unit sphere. Can spherical harmonic expansion be used to analyze the motion trajectory of this point?

[wide spear]

@minor osprey maybe you have better ideas about this

[minor osprey]

Relative to the fixed joint, yes. It might be better to take the Euler angles page Ange sent and view the trajectory as relative to the joint (consider that the origin) to acquire the necessary 6d Euler angles from the odometry data you are trying to analyze

[sage river]

I understand, thx!

[next hollow]

Hey guys am confused about something, for the bisection method is the order of convergence 1 or 1/2,our professor says it's 1/2 while in some other sources it's 1

[wide spear]

You should be careful about the terminology

Order of convergence vs rate of convergence

[terse harness]

Can anyone help me? I'm trying to perform an error analysis on the following algorithm but i just cannot seem to find any way to start this

[wide spear]

Floating point error?

Have you seen a floating point error analysis done before

[steep stone]

i think u can compare error bound with Horner's method and perform backward error analysis

[terse harness]

only on the naive summation and product, this was like the third exercise i attempted

i don't understand anything about this course, while PDE feels like a breeze

i'm so confused

[stoic thicket]

anyone know which techniques are commonly used for fluid solvers with large domains?

[fathom rain]

finite volumes or finite differences

[wide spear]

This is a complicated question

It depends not just on the domain size but also the flow you are trying to capture

[stoic thicket]

I suppose finite volumes, in particular i had in mind eulerized cells.

so, high cell-counts

[wide spear]

Also the grid spacing matters

Different methods are better for LES vs DNS

[fathom rain]

I don't know what "Eulerized cells" means. If you want high accuracy (and domain is appropriate) you might want to use DG or spectral methods.

[fathom rain]

DNS is not really applicable for "large domains".

[wide spear]

Well also don’t know how large of a domain they want

Large domain means different things to different people

[normal sparrow]

Ime DNS (and plain RANS) are only used to test and validate before scaling up (decreasing grid size, time step, duration etc). That is it is very rarely used in practice for an actual end-to-end simulation of a novel structure

You'd either want a modified RANS (e.g. PANS) or the many variations of LES

I'm happy to be corrected

[wide spear]

We don't know what sort of fluid sim they want to do at all

So I don't think we can really say one way or another what they should be doing

[stoic thicket]

just so the specifications are a bit more clear, this is intended for a subsystem for a game. by "large domain" I meant something on the order of 2^4096 cells.

the simulation isn't intended to be too extremely accurate, only convincing.

the cells themselves are meant to be roughly half meter cubes.

[wide spear]

Wait didn't I suggest lattice boltzmann for you before

Also are you sure about 2^4096

4096^2?

[stoic thicket]

side length of 65336 cells

[wide spear]

2^4096 is more than the number of atoms in the observable universe

[stoic thicket]

I think I fucked up the math then hold on

[wide spear]

65536 squared or cubed?

[stoic thicket]

cubed.

[wide spear]

(2^16)^3 is 2^48

[stoic thicket]

my bad

[wide spear]

Ok at any rate

[stoic thicket]

but yeah LBM is probably the approach I'm going to use but it feels like I'd either have to 1. vastly reduce this domain size

or else find some way to compress the represented data

[wide spear]

Well this is for a video game right

[stoic thicket]

correct

[wide spear]

Do you need the entire domain always?

Or can you do chunking and only evolve the fluid in the portion currently relevant to the gameplay

[stoic thicket]

no, and really the size was only ambitious. I was just curious if this was somehow a known-solved problem

[wide spear]

Ok so the way you solve the problem of very large domains doesn't really vary too much by the numerical method you use

You just parallelize your fluid solver

Each parallel rank gets a subdomain of the entire thing

[stoic thicket]

so rank in this case would be like

octation?

[wide spear]

?

No parallel rank meaning like

An mpi process

With an independent memory, perhaps running on a full computational node

[stoic thicket]

I see

[potent river]

hey do you people have to use numerical methods for FFT, or do you just use inbuilt function?

Also 🥲

[wide spear]

The inbuilt functions are numerical methods

I'm not quite sure what you're asking

[normal sparrow]

do you mean to ask whether you should implement fft yourself or use a library's implementation?

[potent river]

yes

[potent river]

and we have different kinds of those methods as well

[cold linden]

never realized how cracked game devs were

ran across a blog post of RK4 visualization from a game dev

[normal sparrow]

woahh

can you send?

[normal sparrow]

well we would need a lot more context and domain knowledge to help with that but some rules of thumb:

the library implementations are often very good and have settings that can be tuned

unless you're doing research level work or deploying on an extremely resource constrained device

likely library implementations are good enough

for your question of the different libraries

they will have their tradeoffs

e.g. it says here [1] that scipy fft is sometimes more performant on larger data but in general the visible differences are likely negligible for everyday use cases

[1] https://stackoverflow.com/a/6363572

There are variants within libraries too such as FFTW but again you would just need to read up on them and see if it would be best for your case

[wide spear]

Ok also so

If you need a FFT, you should start with a simple preexisting implementation

If it’s too slow, then you further investigate more advanced implementations

[crisp lake]

I'm not sure those two things are mutually exclusive.

[cold linden]

https://www.haroldserrano.com/blog/visualizing-the-runge-kutta-method

[wide spear]

Doesn’t Wikipedia have a similar picture

[cold linden]

does it

ngl i dont even think the visual helped my intuition that much

[wide spear]

But yes the picture doesn’t help that much

[cold linden]

i just tried to think abt it like

a more complicated version of eulers method

my professor is making us memorize the butcher table for RK6

😭

[wide spear]

[normal sparrow]

[normal sparrow]

ah kinematics

interesting to see it on a game dec post fs!

[cursive heron]

I'm struggling to remember all the theorems about error estimates of stuff. For example this one about the interpolation error of Lagrange polynomials. Is there any intuition for this, or do I just have to memorize the formula?

[wide spear]

I feel like this is pretty intuitive

[cursive heron]

It is? What's the intuition?

[wide spear]

Well the pi_n+1 is from the use of the lagrange basis right

It's exact at the interpolation points

And degree n polynomial interpolation always has error that depends on the (n+1)th derivative

And you divide by (n+1)! because it's like a taylor series

[cursive heron]

hmm yeah, part of that makes sense to me. Why do degree n interpolation have error that depends on the (n+1)th derivative?

[wide spear]

A degree n polynomial interpolation exactly interpolates polynomials of up to degree n, so the values of the n+1 th derivative is how higher degree behavior you have

[cursive heron]

I think that makes sense, thanks it's not completely obvious to me that the formula should be exactly that, but I think I see why those parts are involved at least

The xi variable is a bit weird to me, but I think that follows from Rolle's theorem in some way

[wide spear]

Yeah

[wanton dawn]

One method of proof is to apply $E = I - P_n$ to the integral form of the remainder from Taylor's theorem. The claim is that $E$ applied to the remainder is the RHS

[pine jetty]

L

[wanton dawn]

You can also look into divided differences, which I believe generalizes your result

[cursive heron]

Yeah, I found this: https://en.wikipedia.org/wiki/Mean_value_theorem_(divided_differences)
The funny thing is the proof seems to be almost exactly the same as for the theorem above

[fleet igloo]

How can I compute the operation count of Classical Gram-Schmidt（CGS) method?
I tried a detailed derivation but it seems doesn’t work basically.

The result should be something like O(mn^2), but my derivation shows like it is the stuff at the right bottom of my notebook

[fleet igloo]

<@&286206848099549185> I'll also try google online

[mystic jacinth]

any numerical enjoyers understand what's causing a noise like signal between boundaries when handling this first order wave PDE with leapfrog?

[wide spear]

Unsurprising with leapfrog

Try adding a filter

Like the robert-asselin filter

[mystic jacinth]

would that supress the reflected noise but preserve the second pulse or would it turn the noise into something meaningful?

Ive not covered filters

[wide spear]

Is this for a class

[mystic jacinth]

it is

[wide spear]

Ok no filters then

What have you tried to fix it so far

[mystic jacinth]

just trying to understand the behaviour at boundaries

[wide spear]

What boundary condition are you using

[mystic jacinth]

not necessarily anything , but do want to understand what causes it to arise

[mystic jacinth]

def staggered_leapfrog(t0, t1, dt, x0, x1, dx, A, v, d, x_beg):
    times = np.arange(t0,t1+dt,dt) 
    x = np.arange(x0, x1+dx, dx)
    
    print(f'Running to time {t1} with Courant number {v*dt/dx}')
    

    # f_nm1 == f -> u - 1 , f_np1 -> u + 1...
    f_nm1 = analytical_solution(x, 0, A, v, d, x_beg)
    
    #f_nm1[0]  = 0
    #f_nm1[-1] = 0
    
    f_hist = [f_nm1]
    
    f_n = f_nm1.copy()
    #f_n[1:-1] = f_nm1[1:-1] - dt*v*(f_nm1[2:] - f_nm1[:-2])/(2*dx)
    f_n[1:] = f_n[1:] - dt*v*(f_n[1:]-f_n[:-1])/dx # upwind alternative (actually stable ?)

    f_hist.append(f_n)

    for _ in times[2:]:   
        
        f_np1 = f_nm1.copy() 
        
        f_np1[1:-1] = f_nm1[1:-1] - dt*v*(f_n[2:] - f_n[:-2])/dx
        
        #f_np1[0]  = 0
        #f_np1[-1] = 0
        
        f_nm1 = f_n # u-1 <- u   
        f_n = f_np1 # u   <- u+1

        
        f_hist.append(f_n)

    return(np.array(f_hist),x,times) 

[wide spear]

[mystic jacinth]

there was no difference in clamping the edges to 0 and I sure hope this implementation is correct lmao

[wide spear]

I will not be looking at your code

[mystic jacinth]

yeah no worries it just has the boundaries in there

[wide spear]

Can you describe the boundaries

[mystic jacinth]

set as this initially

that gets bootstrapped with upwind

and then the rest are leapfrog

[wide spear]

That's the initial condition

[mystic jacinth]

the edges at 0 and arr[-1] are clamped to 0

I know but just full context

[wide spear]

Your x looks like it's going from 0 to 60?

As a general note, I think that the type of 3d plot you've displayed is very unhelpful to look at

I would recommend making gifs if possible

[mystic jacinth]

yeah, I dont think that would change the behaviour tho right? I mean there is a spatial boundary always there but it wouldnt change how it traverses

[mystic jacinth]

it was just a quick thing to provide some context

ill make a gif rn

@wide speargif do be giffing

[wide spear]

What fascinating behavior

How are you computing the spatial derivative

[mystic jacinth]

the exact code in the method?

- dt*v*(f_n[2:] - f_n[:-2])/dx

just this, not sure if thats what you meant mb

[wide spear]

Try using a centered difference

[mystic jacinth]

I dont think I explicitly calculate the centred difference, but it is a simplification of this

[wide spear]

Are you saying this is a centered difference?

Why isn't it divided by 2*dx then

[mystic jacinth]

it would just be multiplied by 2 afterwards anyways

the full step calculation in code looks like

f_n[1:] = f_n[1:] - dt*2*v*(f_n[1:]-f_n[:-1])/(2*dx)

so I just dont bother with the 2's

[wide spear]

Ah ok

[mystic jacinth]

but yeah even with different ways of bootstrapping it, it just has the same behaviour

if I enforce f[1] = 0 it even reflects by pi , behaving like a physical boundary

but yeah not really sure what causes it to happen

[wide spear]

Ok I do believe that this is a result of leapfrogging

Try swapping the time stepping

[mystic jacinth]

will do, might take a second

    for k, _ in enumerate(times[2:], start=2):

        if k % 2 == 0:
            
            f_np1 = f_nm1.copy()
            f_np1[1:-1] = f_nm1[1:-1] - dt*v*(f_n[2:] - f_n[:-2])/dx
        else:
            
            f_np1 = f_n.copy()
            f_np1[1:] = f_np1[1:] - dt*v*(f_np1[1:] - f_np1[:-1])/dx

        f_nm1 = f_n
        f_n   = f_np1
        f_hist.append(f_n)
```This is just from alternating methods

SOOO WEIRd

I have a tutorial tomorrow with the prof so I might just ask him to see if he knows anything too

[wide spear]

Ummmm are you sure that is implemented correctly

Why does the wave disappear

[mystic jacinth]

I think(????) that should happen if its travelling off the boundary

im not entirely sure to be honest , ill give him a shout tomorrow , and see if anything comes from it

thanks for the convo <3

[wide spear]

In my experience a dirichlet boundary condition is usually reflecting but that might only be for the second order wave equation

[mystic jacinth]

I hadnt enforced f[0] = f[-1] = 0 , Im just gonna attribute it to leapfrog being non-dissipative

[mystic jacinth]

true, depends on the PDE

[wide spear]

What you can get with leap frogging is alternating time states splitting apart

So like

u(x,t_0), u(x,t_2), and so on will form a solution

And u(x,t_1), u(x,t_3), and so on will form a solution

[mystic jacinth]

so Odd–even decoupling ?

[wide spear]

Yeah

[mystic jacinth]

hmmmmmmmmmmmmmmmm

[wide spear]

This is what filters fix

If you're interested you can search for the robert asselin filter

[mystic jacinth]

im gonna take a look

there seems to be some sort of 'computational node'

that is refenced a lot

and I remember reading odd even decoupling along with it

[wide spear]

A computational mode is a wave that is not physical but gets propogated/enhanced by the numerical scheme

Lots of numerical methods, if you're not careful, will have a computational mode with wavelength 2 delta x

[mystic jacinth]

@wide spear thought id ping just in case your interested

taking the absolute value of the backwards propagation seems to restore the original pulse, would the amplification factors complex parts become more dominant on the first bounce leaving whatever mess that is on the left?

[wide spear]

What are you taking the absolute value of

Which backwards propagation?

[mystic jacinth]

the value of the numerical solution at all t's

[mystic jacinth]

i didnt know how else to describe it

but the noise travelling with ~-v

[errant crystal]

Hi, it's the first time I've written in a channel, I'm passionate about number theory, any advice on how to study it? (I'm Italian so I'm using Google Translator sorry for the translations)

[wide spear]

Do you have a motivation for doing this

[wide spear]

Number theory or numerical analysis? Wikipedia tells me that number theory is teoria dei numeri and numerical analysis is analisi numerica

If you do want number theory, then you are in the wrong place and should go to #elementary-number-theory

If you do want numerical analysis, I'd be happy to provide a few book recommendations

[errant crystal]

then I got the wrong channel

[wide spear]

[mystic jacinth]

other than the amplification factor being complex valued no

I figured that maybe the boundary causes the computational node to become dominant, and as long as the phase information hadnt been polluted with errors its magnitude should be equal

I think the phase errors would just arise from having a grid spacing that doesnt resolve with high enough resolution for a given v and dt

I did find it really coold that it did restore the original profile tho

[cursive heron]

Are Newton divided differences polynomials the same as Lagrange polynomials, just calculated differently?

[wide spear]

Newton divided differences and lagrange basis polynomials are two ways of constructing interpolating polynomials

Given a set of d+1 points and corresponding function values, there is a single degree d polynomial that interpolates the data

Newton and lagrange are two ways of getting to that polynomial

[cursive heron]

I see, thanks

[wide spear]

<@&268886789983436800>

Posted in every channel

And appears to be an exam

[wary oriole]

Ha I just reread it today. Two different ways. Lagrange interpolation construct a identity matrix for Ax=b. Where newton construct a lower triangular matrix for the system. Divided difference is another way to get to the system (better conditioning even). Due to special form of A both takes O(n^2) to solve instead of full matrix constructed by plain monomial basis

It’s about two different yet both smart way of choosing basis for the polynomial resulting in special form of basis matrix. Instead of just using 1,x,x^2,…resulting in full matrix

[cursive heron]

thanks for the explanation, didn't know there was a linear algebra way to view this

[spice valley]

Hi, can someone help me with this? I am really confused and don't quite know how to start. So far we learned about Newton-Cotes, Gauss-Christoffel, Gauss-Legendre, Euler-Maclaurin, and Romberg quadratures/formulas. But tbh I have no idea how to use any of them and far less how to use it for a generic polynomial...

[wide spear]

You aren't using a rule you already know

You are deriving a new quadrature rule

Did you learn about how quadrature rules are derived

[spice valley]

Not that I know of tbh

let me take a look at my notes again

The numerical approximation of a given integral is called numerical quadrature. The basic approach of numerical quadrature is to first approximate the function in the given integral by a function that can be integrated exactly. This is often a polynomial. Then a rule is developed, a so-called quadrature formula, which integrates the approximating function, for example the polynomial, exactly. If the interval is long, it can be divided into subintervals beforehand, the quadrature formula can be applied to each subinterval, and finally the integral over the entire interval can be approximated by summing the results on the subintervals. This approach is called a summed quadrature formula.

That's about everything I could find. Everything else in the script is about the different rules there are.

[spice valley]

I don't think I have unless I am missing something

[wide spear]

Ok

So what you do is

You have 4 unknowns

lambda 0, 1, 2, 3

You plug in f(x) = polynomial of degree 0, polynomial of degree 1, polynomial of degree 2, polynomial of degree 3

You get a system of linear equations

Then you solve for lambda 0 through lambda 3

[spice valley]

for this is it enough to put f(x) = c, f(x) = x, f(x) = x^2 and f(x) = x^3? or does it have to be the entire generic polynomial like ax^3+bx^2+cx+d?

[wide spear]

Yes the monomials are enough

Integrals are linear so you can add monomials

[spice valley]

amazing! thank you so much! I'll give that a try

[spice valley]

it worked :D feels almost like magic

thank you again!

[wintry elk]

First quarter of Numerical PDEs is almost over. Two more quarters to go

[fathom osprey]

gg

[wintry elk]

Really just felt like Taylor series, linear algebra, and some coding

[grave spoke]

Which methods are covered in the course?

[wintry elk]

We did Finite Difference, Finite Element, Stability, Consistency, Crank-Nicolson, and some other stuff

[fathom osprey]

gg @wintry elk fun stuff coming up soon 🙂

[wintry elk]

Yeah, we have Stochastic Numerical PDEs next quarter

Spring will be finite volume. I've never really done numerical/applied stuff before so it's been a whole new ball game for me

[fathom osprey]

I took pde znd it was super fun

[grave spoke]

so, with a focus on second-order linear elliptic and parabolic PDEs?

[wintry elk]

Basically, yeah. The pre-requisite class is graduate Numerical Analysis and the graduate sequence on PDEs

[fathom osprey]

like in my phd, new stuff fr me coming frm a ag and cat backgrnd in undergrad

like i did take pde classes

but that was cz theory and fourier stuff n strichz estimates

but the numerics was super new and super fun; wish i could d again

[wintry elk]

I'm enjoying Numerical, but god damn does HW in third year suck

[stable trellis]

I gave up on taking numerical next term bc it's too much debugging

Idk how academics in numerical analysis do this sort of thing

[lime moat]

Lmk what textbook you guys use and what is covered

[wintry elk]

Textbook = Lecture Notes

[lime moat]

Lol

[wintry elk]

It's tough, because I had this guy for numerical. He's very good if you're extremely familiar with the subject. But if you're unfamiliar, it's very hard to follow

So he'll be stating a theorem, and then go on a tangent about how he read a recent paper that actually improves the convergence with slightly stronger/weaker assumptions

And talk about the applications of that paper for like 2-3 minutes

Then go back to doing the proof, but while doing the proof he'll continuously compare the method that we're learning to a method he recently learned

That none of us have ever seen

[lime moat]

🤣

My NLA professor is exactly like this

He will go into very large and interesting rants about the deep relation between Krylov subspaces, spectral calculus, and the Kalman rank condition for operator theory, or the Frechet derivative, or about how Schur's triangularization trivially resembles Jordan normal form or how Ritz-Galerkin is a poor man's FEM, when half the class hasn't even taken PDEs

I kind of prefer it though, I'd rather have more information than I could ever need than not enough

[hollow lily]

Hello. I just joined and i dont know if this is the correct channel but i have a question. Im a cs studnet and im writing a program that approximates the Poisson equation on a round plane, but i dont understand why are we doing this. We already have a function g. Why are we approximating it? I really dont understna dthe sense of executing a program just to find a worse answer than if we plugged points in g

[wintry elk]

You can treat it as a test case: if you do have an exact solution, then how well does the approximate solution work? What is the error given by plugging in an exact solution? There are times in the future where you don't have the exact solution, only the approximate solution

Sure, it will be worse, but how much worse? Is it an acceptable error for your application? If so, why or why not?

Is the approximate solution easier to compute than the exact solution? etc.

[hollow lily]

So im just benchmarking? Is there any other use of for this? Or is this just a sort of test to see how well does an approximation methode work on a known solution?

[wintry elk]

I mean, what is the use of an approximate solution if not to approximate a solution?

If it's known, then that's great. But many times it won't be known

[wide spear]

What do you mean a round plane

Are you solving the poisson equation?

[wintry elk]

I took it to mean the unit disk

[wide spear]

If you're solving the poisson equation then what do you mean by plugging points in g

[wintry elk]

I also took $g$ to be the solution to the (somewhere) given PDE

[pine jetty]

MoonBears-C-

[wintry elk]

Oh D is an annulus

[hollow lily]

Yes sorry did not know that word

[wide spear]

Ok so you know the function g

But you are trying to find the function u

[hollow lily]

yes my question is why

why search for a function u when we know g

[wide spear]

Solving the poisson equation is important in a number of physics problems

[wintry elk]

Ohhh, $u$ is the unknown

[pine jetty]

MoonBears-C-

[wide spear]

Electromagnetism, gravity, fluid dynamics

https://en.wikipedia.org/wiki/Poisson's_equation#Applications_in_physics_and_engineering

[hollow lily]

I have read that but i dont understand why

[wide spear]

Why what

[wintry elk]

It's like saying I know x^3 + 3x^2 + 2x + 1 = 0. Why solve for x when it's just zero on the right hand side?

[hollow lily]

what is the sense of all this if we have g why we are searching for an approximation

[wide spear]

Frequently you don't have g in a closed form

[hollow lily]

hmm

[wide spear]

If you have g given as a data field from some measuring device for example

[wintry elk]

Here's a a toy case, take $\Delta u = 0$, can you find all such $u?$

[pine jetty]

MoonBears-C-

[hollow lily]

does this mean that I should consider delta u as a solution for g?

[wintry elk]

You want to solve $\Delta u = g$. g here is given, u here represents the uknown function

[pine jetty]

MoonBears-C-

[wintry elk]

Here's a more apt example, say $f'(x) = e^x$. Can you say what $f$ is?

[pine jetty]

MoonBears-C-

[wintry elk]

Or if $f''(x) = x^2$. It's in the same vain as these problems

[pine jetty]

MoonBears-C-

[wintry elk]

The goal is to always find the function

[hollow lily]

In the exercise im asked to find, given the boundery conditions at r=1 and r=2 and g, all the values of the points on the annulus that are created when we subdivide the area using an n. For example let n=4 then i subdivide the area in 4 angles 0,pi/2,pi,3pi/2 and 6 radii 1,1.2,1.4,..2. My question is why do we use approximation with the Poiss eqaution and not just plug the points into g so g(0,1),g(pi/2,1)... is this just becouse thhis is the exercise so to see how approximation works or is there another reason? I know that its stupid but i would really like to understand.

[wide spear]

You already know g

You want to find u

[wintry elk]

What is the function u?

[hollow lily]

im sorry i just undersooud that i wasted your time. My initial question was given wrongly but now trying to explain myself i understood the sense. I was not really given g i was given g in a programming sense as an input i dont have a known g but the question was posed as if I have a known g and Im searching for a u.Like u``=x then u=integral(integral(x)). Im sorry for wasting your time.

My question was in the sense why should we aproximate the points that could be g when we have g but i just understood that the sense was the approximation which yall tried to tell me from the start.

Thx

[wintry elk]

Glad we could help

[dawn gazelle]

Hi! I am almost done with my bachelor specialized in numerics and stochastics (probability). I was thinking about leaning more on the numerical anaylsis side in my master and maybe contuniue with a PhD programm e.g. optimazation

Does anyone have experience or has any tips in advance?

[lime moat]

Try to do research as an undergrad, look at 50-100 programs you think might have interesting research in your area, look at what professor research is done there, do your GRE and GRE Math, start submitting applications. Good luck

Look at 50-100 as in look at the programs and see if the research area fits your needs and desires

Also important to look at things like if the location, etc match your desires

[wintry elk]

Try not to get too married to a subject early on. What you find interesting now most likely isn't what you found interesting three years ago. Keep an open mind, but push towards your goals. A MS program is usually a good place to get your feet wet with either a MS thesis, or take a broad range of classes to see what you wanna do

I switched what I wanted to do from CC, to University to MS program, to research between PhD, and now in PhD I do something different from everything before

Just follow your interests as they come, and work on math when you can. Bonus points if you can find some research mentor

[dawn gazelle]

I see, yeah makes sense, thanks you both

[brave crypt]

Who here is a numerical analyst?

[normal sparrow]

It's generally better to ask your question (whether a specific one about a numerical problem you have, or advice about career etc) than to ask for specific *people*

[halcyon thunder]

What classes do you need to take to do research, other than the standard 1 year of numerical analysis?

[wide spear]

Depends on what sort of research you want to do

If you are doing numerical analysis applied to some subject area, then you should take classes in that area

If you want to do numerical linear algebra, then functional analysis is useful

[haughty ermine]

Hey guys

Im trying to read up on numerical methods on an intro level

Any recommendations for texts?

Books you liked, I'll mainly be using Python for tools but MATLAB is OK as well

[late kettle]

"An introduction to numerical methods and analysis" by james f epperson is a lovely introductory level text. I'd also go for the first course lecture notes by Ralph Hiptmair at ETH Zurich (||C++|| )

[plain pendant]

I took a lot of numerics courses and have never found a textbook I actually liked

[wide spear]

Iserles is nice