numerical-analysis

[worn night]

yea

but you also have to note the property that the each b_i represent the SUM of the deviation values of each row essentially

i.e. for first row, in the case of 1 coin, there is guaranteed NO coins that are counterfeit, since b_1=0

but for two coins, I suppose you could have one positive one negative (which would have to be x_1=-.2 and x_7=.2

and for the case of still no counterfeit coins in row 1, you basically get down to the case where x_6=0 and x_2=.2 and x_4=.2

am I thinking about this wrong

[worn night]

ohhh wait

I can combine these to just get one system, basically make everything that isn't x_1, x_2, x_4 or x_7 0

then show x_1=-.2, x_7=.2 x_2=0, x_4=0 a solution, then x_2=.2, x_4=.2, x_1=0, x_7=0 a solution

This doesn't obey addition property of subspaces right

because the entries in which they have zeros can vary and what not and it becomes less sparse?

[worn night]

and this isn't a norm because it doesn't obey the scalar properties of a norm?

[proper valve]

Yes you can show it

[drowsy roost]

here are the quadratics $p$ that minimize $|e^{-x^2} - p|_p$ for $p=2$ and $p=\infty$. is it expected that the infinity norm would suck as hard as it does compared to the 2-norm (least squares)

(on [-1,1])

[pine jetty]

winston gergill

[wide spear]

You use p to mean two different things in that expression

[drowsy roost]

you right

[wide spear]

On what domain are you minimizing

(-1,1)?

[drowsy roost]

[-1,1]

yeah

[wide spear]

Your minimizer for the 2 norm error clearly has a smaller max norm error than the purported minimizer of the max norm error though

[drowsy roost]

indeed

but the only thing changed in the calculations was normalizing a set of orthogonal polynomials according to either the 2 norm or max norm

i guess a little context: used gram schmidt to get those ortho polynomials ${\phi_n}$, then normalized those wrt a norm. call those polynomials ${\psi_n}$, then compute the minimizing quadratic as $\sum_{j=0}^n \beta_j\psi_j$, where $\beta_j = \langle f, \psi_j\rangle = \int_{-1}^1 f(x)\psi(x)dx$ (here f = exp(-x^2))

[pine jetty]

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

[drowsy roost]

well that's the tex i wanted lol

clearly does fine for the 2-norm, not so well for the sup norm

[wide spear]

Have you tried it for the 1 norm

[wide spear]

Looks wrong

[drowsy roost]

yes

eh it's probably just late and there's a symbol wrong somwhere

[sleek sundial]

I’ve got an ODE that I’m trying to solve numerically:
$$w''=2tw'-10w, w(0)=0, w'(0)=120$$
I know that the solution to this is $$32t^5-160t^3+120t$$. The issue is that I need to solve it for t from 0 to 10, and for Runge-Kutta methods to be stable there, I need stupidly small step sizes ($10^{-7}$ at least, and my laptop simply doesn’t want to calculate that), otherwise it blows up. What change of variables could help me here?

[pine jetty]

BIGfoot496

[sleek sundial]

I tried doing $t=\ln{(x+1}$, which results in the equation
$$(x+1)w''+w'=2\ln{(x+1)}(x+1)w'-10w(\ln{x+1})$$
but the solver says that it cannot solve a delay differential equation that depend directly on time variable.

[pine jetty]

BIGfoot496

[hoary light]

your problem is linear. If you cast it as a first order 2D system, have you checked the region of stability?

[sleek sundial]

It depends on the step size, no?

(I’m bad at this, sorry, if I’m being dumb)

[hoary light]

Maybe more to the point, check the eigenvalues of the 2x2 matrix (it depends on t) and from that you can hopefully choose an appropriate time stepping method

[sleek sundial]

Well, the eigenvalues are weird. They are complex until t becomes sqrt(10), and real afterwards

I don’t really know what that entails

[hoary light]

it might entail using a different method for different time intervals

often all that matters is the real part, however

you're familiar with checking the region of stability for Euler, backward Euler, etc?

[sleek sundial]

If I understand correctly, the main issue with this particular problem is that the solution grows really quickly, and the step size for RK to be stable needs to be significantly smaller than 1/max(w) ~ 10^-6

[hoary light]

I see. I guess I cannot comment on that without actually working something out on pen and paper myself. But basically what I'm getting at is that you should be thinking about how to use a different method, not trying to change the equation so that you can use a particular method

Some clever transformation of the ODE won't lead to any general principle that can be used in other problems.

[sleek sundial]

The problem specifically said to use Euler, Implicit Euler and RK of orders 2-4

[wide spear]

Can you use implicit RK

That should avoid the small time step issues

[sleek sundial]

Eh, I tried. Still unstable

It does a decent job for t<7.5, but after that it blows up to infinity

As in the error becomes many orders of magnitude larger than the solution

So, basically, the only way to solve the problem with the methods they ask is to reformulate the equation somehow

[delicate elk]

#linear-algebra

[idle quartz]

Do you use petsc or Trillinos

?

[grave spoke]

dedicated numerical analysis channel

[minor osprey]

We also have #modeling and #optimization

[grave spoke]

yeah, I saw them

[merry ginkgo]

hasn't this channel already existed but was just renamed?

[minor osprey]

It used to be #applied-computational-math

[merry ginkgo]

ah right

is the name change just to make the focus of the channel more specific?

[minor osprey]

As a part of making the server more friendly to Applied people, we separated channels into narrower categories

Auto correct moment

[sleek elbow]

casserole math

[merry ginkgo]

ah okay so modeling and optimization used to go here?

[minor osprey]

Yeah, this channel used to encompass all of applied math

[merry ginkgo]

I see, good change imo

[minor osprey]

Which is like saying all algebra goes in #groups-rings-fields

[exotic pewter]

does anyone know how numpy samples gaussians, i.e. what is the procedure python uses to approximate a gaussian?

[uneven thicket]

I'm not sure exactly how numpy does it; however, I strongly suspect it uses "inverse transform sampling"

Certainly for univariate distributions with continuous support this is the most straightforward (and often fastest) technique

In the case of the multivariate normal distribution you can use an analogous technique

I can't remember off the top of my head the exact details

[brittle atlas]

For multivariate would they transform from standard multivariate?

[uneven thicket]

Yes

[brittle atlas]

What about general multivariate distributions? A while ago I was wondering how to generate CUE eigenvalues

[uneven thicket]

Sample from the standard multivariate then transform using the Cholesky decomposition of the covariance matrix of the target distribution

iirc

For general multivariate distributions you will usually need something like Metropolis-Hastings

[brave crypt]

MCMC methods basically

[uneven thicket]

Exactly

[brittle atlas]

Hmm ok I will look into that thanks

Right now no one's using the uni servers since winter break so I just have them making unitary matrices but it's of course very slow

[wide spear]

Didn't I suggest that you parallelize it

Or put it on gpu

[brittle atlas]

yes but it's still going to be way slower if I have to make a giant unitary matrix and then solve the eigenvalues

also I have not figured out how to run on gpu

can't install the right driver on linux and the easy cupy package doesn't have non-hermitian solver

but I did finally get access to my university's server that has a gpu so I might experiment with it over winter break

(numpy does parallelize np.linalg.eigvals at least)

[wide spear]

The lapack routine you want is cgeev

[brittle atlas]

I don't know if I want to know how to directly call lapack

[wide spear]

I promise it isn't that scary

[brittle atlas]

what if the only language I can sort of code in is numpy

and bibtex

[wide spear]

washingbear learns c++ over winter break

[brittle atlas]

too much power

[random canopy]

lapack good and eigen good

[brave crypt]

I suggest julia for if u want to do big computations

[exotic pewter]

this is cool, but this then begs the question, how do they "uniformly" sample from the interval [0,1]. They only have limited precision, i.e. number of bits, so how does that impact the precision of the inverse transform sampling.

[dense echo]

I'd expect most libraries generate uniform integers in [0,2^53) and divide them down. It just seems too much of an invitation for subtle biases if you try to exploit the higher precision in the lower end of [0,1).

[uneven thicket]

Great question that I am unfortunately unqualified to answer in sufficient depth. (Pseudo)random number generation is a huge field with a lot of depth. As far as I know, the most common general purpose pseudorandom number generator is the Mersenne Twister, though.

As Troposphere alludes to, I expect a fair bit of thought has to be put into dealing with the fact that floating point numbers are themselves not a "uniformly distributed" subset of the reals.

[vale willow]

Good books on numerical analysis? Specially interested in numerical integration

I found the one by Stoer and Bulirsch, which seems pretty nice

[fathom rain]

not a book but this post might be of interest https://fredrikj.net/blog/2017/11/new-rigorous-numerical-integration-in-arb/

[feral lark]

"An introduction to numerical methods and analysis" by james f epperson is a lovely book.

[brave crypt]

By numerical integration, are you more into integration in general or solutions to differential equations?

Any decent numerical analysis text should cover numerical integration until quadratures

But not all of them cover monte carlo integration

[stiff dawn]

Is it typical to know matrix decomposition algorithms?

[plucky kayak]

it's probably important to have a rough idea how you can efficiently compute decompositions, but I would say you can spend your entire life studying more and more efficient versions and discovering new algorithms

it also highlights the pros and cons of using some decomposition

[stiff dawn]

I dont really understand this statement

Ive had impressions that matrix decompositions are just all known

so what do you mean by spend entire life studying more efficient algorithms?

sorry for my ignorance

[plucky kayak]

well you said "matrix decomposition algorithms", so I assumed you meant computing various decompositions

because usually the standard algebraic way to compute some decomposition that gets taught is far from being fast or stable in practice

[stiff dawn]

oh i didnt know

i also never learned any decompositions

so im trying to find resources on the internet to learn stuff like cholesky,lu,qr

[plucky kayak]

in that case I'd just look into the general idea behind computing them and focus more on the applications of these decompositions and their existence, stability
in general there's only a fairly short list of decompositions to know - SVD, Spectral, Cholesky, LU(P), QR and perhaps something else I'm forgetting, maybe complex numbers people enjoy using polar decomposition, but I never looked into it

the actual computation of each of these is it's own kind of science, let alone computing sparse decompositions or low rank

[stiff dawn]

cool

excuse me about the questions but i have had interest in the algebraic side of things and never took the time to branch out how I feel proper

thanks for educating me

[hoary light]

It's worthwhile to know something about the existence and mathematical meaning behind the basic decompositions (LU, QR, SVD, etc). Learning algorithms is a much more in depth endeavor and depends on how serious you are about diving into applications. In practice, decomposition algorithms are area-specific because many matrices are structured rather than arbitrary.

[stiff dawn]

i dont really know what QR is for

but i know LU is helpful for solving linear systems because it reduces the situation to forward substitution and backward substitution.
SVD is helpful for visuallizing how a matrix transforms vectors? i can only assume QR which I know to be orthogonal something has to deal with eigenvalues but idk

because ik main reasons for factorizing is to either make operations like exponentiation easier in case of diagoanlizing, solving systems of equations, or finding eigenvalues

however i cant compute them and im not sure about this numerical stability thing

[plucky kayak]

QR-algorithm is based on QR decomposition and is used to compute all eigenvalues of the matrix and it's actually quite nice algorithm when you need all eigenvalues, you can also use QR algorithm to solve least-squares normal equations without introducing squared condition number - A^T A x = A^T b turns into Rx=Q^T b where R is upper triangular

tbh it's not a big deal if you can't compute any decomposition by hand, it's not really worthwhile spending time on this

SVD is mostly interesting as a low-rank approximation to the original matrix because by Eckart-Young theorem we know that first k left and right singular vectors and singular values (that is cut SVD) is best rank-k approximation to the original matrix and it provides a direct way to estimate the error

you can also generally tell the approximate rank of the matrix by looking at singular values in log plot and how rapidly they decay - rapid decay in log scale is something that decreases quicker than straight line, low rank matrices can have rank far less than their dimension therefore you can replace them by low rank SVD and loose no information

[stiff dawn]

i was curious about this especially because i kind of felt insecure about this

at the same time i have never been required to decompose a matrix

low rank approximation means?

[plucky kayak]

well if you have a matrix say 1000x900 and it's rank is only 200, then you can pick 200 linearly independent vectors and that'll be enough to recover all the information in the original matrix

while the dimension is high, the rank of the matrix is comparatively low

[stiff dawn]

yeah so what is dimension supposed to mean in case of a matrix? a matrix isnt a vectorspace, so are we saying the dimension of the vectorspace the matrix is in?

[plucky kayak]

here dimension is just how many rows and columns a matrix has

[stiff dawn]

okay lol

so svd gives a way to approximate the original matrix as if it were a lower rank

[plucky kayak]

yeah, though a lot of matrices, especially really huge ones, are low rank in fact

[stiff dawn]

thats weird

id assume that sparse matricies are high rank

ig ur just talk about matricies with elements chosen randomly between some range

[plucky kayak]

it might sound weird, but matrices filled with gaussian noise are far from low rank

[stiff dawn]

wtf

so like

i know another cool thing in numanalysis is sampling algorithms

i want to learn a bit soon too

but if you sample you elements of your matrix from gaussian then they are typically high rank?

is that what you mean by gaussian noise?

[plucky kayak]

no I mean if you were to generate some matrix with random elements which are gaussian noise, then that matrix won't be low rank

there are of course examples of very sparse yet full rank matrices such as diagonal, bi-diagonal, tri-diagonal, etc - those are nasty in that sense (though it's very nice to solve systems with those)

[plucky kayak]

no those are different concepts, sparse means there's lots of zeros in the matrix and low rank doesn't make assumptions about which elements are filled in and which aren't

[stiff dawn]

oh ok

is low rank a precise term?

[plucky kayak]

in the general sense no, there's more precise terms such as numerical rank that uses epsilons to measure the rank that you would find using a computer (because in real life we can't compute rank exactly anyway because it's very sensitive to small perturbations)

[stiff dawn]

so i tried looking up gaussian noise

but it seems to be represented by random vectors or stochastic processes

[plucky kayak]

yeah you just take random gaussian vectors

gaussian noise just means you're sampling from gaussian distribution if that's what you're confused about

[stiff dawn]

ok yeah i thought this

but i phrased it poorl

ty for the education

[plucky kayak]

here just generated random 1000x500 matrix and all singular values are well above 10, while we would hope most of them would be nearly 0 (for low rank)

[long lake]

this is called the eckart-young thm if you're interested (or eckart-young-mirsky)

[stiff dawn]

thanks

[long lake]

gilbert strang has a good video on it

[upper cape]

Hello I am currently working on the network simplex algorithm where the following question came to my mind.
Is it reasonable asking how "stable" the network simplex is, i.e lets say we have a network/graph and the network simplex yields a optimal solution, what happens if we add a node with some arcs to our network?
I wasnt able to find anything about that nor do I know if it is even reasonable asking such a question because:
How would you measure how the solution changes? like something reminiscent of a norm for a directed weighted graph.
If anybody has some idea or encountered this it could help me alot.

[plucky kayak]

perhaps what you're interested in is sensitivity analysis via lagrange multipliers, once optimal solution is computed they tell you how restrictive each constraint is or how much more progress in decreasing the objective you could've made if not for the constraint

[dense glen]

oh this channel exists?

Currently trying to implement the Remez algorithm but having some issues lol

[dense glen]

my implementation is not working

[inner briar]

Okay i'm having a little brainlag here i think. Say i have a fixed point problem $\Phi(x)=x$ and say i know that a fixed point exists, call that point x^. Is $|\Phi'(x^)|<1$ enough to say that the fixed point iteration $x_{k+1}=\Phi(x_k)$ converges to said fixed point for every starting value $x_0$?

[pine jetty]

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

[inner briar]

whatever

|\Phi'(x)|

i am being asked to prove or disprove that f(x) = e^-x has a unique fixed point and the fixed point iteration converges for every starting value to said fixed point for context

the lecture notes are being very vague here

[random hornet]

Every starting point x_0 sounds morally wrong, at best I would like to think this could work in some neighbourhood of x_0

Because like

The function could explode outside some neighbourhood

[inner briar]

yeah i mean it does that's the thing

but the task is to comment on said convergence for every possible real starting value

[plucky kayak]

yeah doesn't that reduce to saying that newton method is a special case of fixed point iteration and that famously diverges for quite a lot of starting points

[inner briar]

yeah that's the thing, my lecture notes say that the abs value of the derivative at the fixed point decides about the convergence or divergence

but i think they have to choose the starting value in some "good" neighbourhood for that

but i'm unsure how to show that for the given example

the fixed point is at like 0.567 and i think choosing a starting value smaller than 0 would give me issues since the function blows up there

i'll just try around a bit for now i guess

[idle quartz]

Wasn't there an applied computational mathematics channel?

[plucky kayak]

it got split into #numerical-analysis , #modeling and #optimization , the history of that channel is in #numerical-analysis

[dense echo]

For a negative starting value the iteration blows up at first, but right after that, a large positive value will always map to a a nice small output.
After two iterations you're sure to be in the interval (0,1) no matter where you started.
After three iterations you're in (e^-1, 1).

And once you're in (e^-1,1) the Banach fixed-point theorem applies.

[vale willow]

,w definite integral from 0 to 1 of 1/x*cos(log(x)/x)

is this accurate?

[pine jetty]

Results provided by WolframAlpha

Click here to refine your query online
Upgrade to WolframAlpha Pro!

Requested by Croqueta#3385

[dense glen]

I’ve seen this integral before

[vale willow]

[dense glen]

Is it right? It should be

Yeah it is

[vale willow]

How do you know?

I was skeptical at first because for entirely non-mathematical reasons I suspected that wolframalpha would not be able to handle it

[dense glen]

See https://mathworld.wolfram.com/Hundred-DollarHundred-DigitChallengeProblems.html problem 1

[dense glen]

Various other methods to compute the integral are given

[vale willow]

nice, thanks

How do they obtain that series?

That sum below seems to diverge

I dont understand what thats supposed to be

[dense glen]

I think it does diverge but you can get a value out of it - if that makes sense

Similar to how 1 + 2 + 3 + … = -1/12 by the zeta function

[vale willow]

coool

[dense glen]

I’m curious as to how you found this integral lol

[vale willow]

Professor sent this problem

[dense glen]

Which professor?

huh

interesting

[vale willow]

I think they didnt know wolframalpha had it as folklore haha

[dense glen]

Wolfram alpha has gotten a lot better since these problems have been posed (2002 I think)

[vale willow]

Its not of any class, they post problems and if you solve them you get points, its like a tiny competition thing. Pretty fun

[vale willow]

tbh, I have no idea how wolframalpha works

[dense glen]

For a numerical analysis class?

[vale willow]

no, they post problems of any topic

And any student can solve them

[dense glen]

Oh interesting

Strange

[vale willow]

and when one gets solved, they post more

[dense glen]

Is it a class you’re in? I wanna join now haha

[vale willow]

its like a tiny competition where every student can participate. So they post three problems (in a webpage), when a student sends a solution, he gets a certain number of points which is the number of days that have past since the problem was posted, and then another problem gets posted. But we have zero interaction with the people who post these problems lol, its not something serious

[dense glen]

Oh fair fair

[inner wharf]

Hi guys! I am interested in massaging a function to make its second derivative closer to a certain other function! Probably need to go in increments since the target ψ'' also changes WRT the function, although not as much.The main purpose is in 3D, but I think this applies to 1D as well.

2nd deriv is curvature. So, I figure if the target ψ'' is too low, we can lower ψ at that point a bit, which raises its curvature there. Perhaps by a constant multiplied by the difference between ψ'' target, and ψ'' as measured from ψ using finite-difference etc.

The good news: This actually, on initial attempt, works really well. In 3D, and 1D. (1d to help troubleshoot this; works fine I think) However, if the nudge amount is too big, the function goes crazy. I can post a pic if that would help. And using small nudges can be too computationally intense since you have to do a lot of them to get a good result.

Are there any approaches you recommend over this proportional nudge? It was the first thing I thought of, and provides great result, but is prone to chaos, and is computationally-intense. Thank you!
(Yes, I am solving the Schrodinger equation numerically for arbitrary potentials)

(Re-post from the calc channel; I think it's probably more relevant here)

[deft citrus]

Dear @inner wharf ,

Since you've already tried some attempts,
I don't think I can get a better idea about it without actual running.

Nevertheless, let me try.
Let F be a primitive function of a primitive function of the target second derivative.
We need to make (psi - F) be flat without too much distortion in psi.

Now the motivation from the physics told me to define,
"energy" E(G) := integral | (psi + G - F)' |^2 + (constant A) * integral | G |^2
and then minimize the energy(by SGD or any optimizer you want).
The only parameter here is A, of course.

If you gonna try it, please let me know whether it works.

[inner wharf]

Thank you very much for the detailed approach! I will try this, and post back here with results

Btw. Pre-nudge (Sum of 2 analytic fns). 2d slice of 3d:

[deft citrus]

Nice visualization... how did you draw it?

[inner wharf]

Post nudge: This appears (I think) to actually be a valid wavefunction, since it matches the schrodingern eq

I'm using a custom graphics lib, using the WGPU library in Rust

(Started with matplotlib, but it was too slow, and I couldn't make it interactive)

[deft citrus]

Cool. Thank you

[inner wharf]

And with too much:

[deft citrus]

lol

[dense glen]

i have no idea what is wrong with my implementation

have some code trying to implement the Remez algorithm and I've looked at the wikipedia as well as some papers, but I am not getting any kind of convergence

idk if anyone wants to try to check my code, too niche

[nimble garden]

Any good resources for solving PDEs given a set of boundary conditions through finite difference numerical analysis?

[dense glen]

https://lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/ ?

[silk axle]

is there any way to numerically verify that a function is bijective?

say i just have some black box f(x) and the domain/range it is valid over -- is there any way, up to i suppose floating point error/precision, verify that it is bijective?

i suppose one can just iterate over a range of numbers but that seems a little unscientific

[dense glen]

is the function continuous?

[silk axle]

hmm in my case no

piecewise continuous

[dense glen]

hmm not sure

[deft citrus]

Since it is not regular in any extent (linear, differentiable or even continuous),
I think the only way is a brutal, discrete method to find-out the collision,
which is known as "birthday attack" in the cryptography.

T be any data structure for easy sort & insert. (like a red-black tree)
for iter = 1, ... , N
pick a random input a_i
compute an output f(a_i)
cut off f(a_i) up to some precision, so make it an integer b_i
check if b_i in T
if there is, f is not bijective, break
if not, insert b_i to T

This is very brutal, but quite good because:
if f is not bijective and you iterate N times,
the probability to find out output-collision pair is roughly proportional to N^2.
(cf. birthday paradox)

[silk axle]

Hmmm yes I suppose there’s no better method. Thanks for the info!

Just curious, if it was known to be continuous, would you suggest something else instead?

[deft citrus]

If it was continuous, it is bijective iff monotone,
so I think it would be more efficient to check whether it increases, or decreases on random increasing inputs without any exception.

[warm otter]

how much code is it?

[faint patrol]

Is there something like continuous least square method for 2 variable functions?

[idle quartz]

I have been following dealii step 15 https://www.dealii.org/current/doxygen/deal.II/step_15.html tutorial for solving a non-linear scalar PDE where linearization is dealt with using a newton rahpson technique.
I want to solve a vector PDE, is it too different? Does the directional derivative become tensorial when dealing with vector valued PDEs?

[hoary light]

The value of a directional derivative lives in the same space as the output values of the original function

Because the directional derivative tells you how much the function changes by, if you perturb the input in a particular direction

[dense glen]

uploading the code i have - the context is i am trying to solve https://en.wikipedia.org/wiki/Hundred-dollar,_Hundred-digit_Challenge_problems problem 5

i must be doing something wrong

[brave crypt]

https://cdn.discordapp.com/attachments/791350441408987146/1058895718192001164/toyota.mp4

[dense glen]

toyota

[minor osprey]

Why is this here.

[brave crypt]

toyota

[dense echo]

Keep random spam in #chill, please.

[dense glen]

finally got the 10 digits lol

[hoary light]

those are great problems. maybe I'd enjoy solving them too

[dense glen]

5 is easily the hardest

[dense glen]

Well actually maybe not - but I’ve implemented solutions of 9/10 of them, good overview of numerical analysis

[hoary light]

I think first I'd have to prove that the reciprocal of the Gamma function is holomorphic

[dense glen]

It’s at least meromorphic

[hoary light]

I want holomorphy on the unit disk so that I can use the maximum principle

[dense glen]

This is a well known fact since the gamma function has no zeros, the reciprocal of the gamma function is therefore entire

“well known” = on the wiki page

[hoary light]

haha, yes I guessed as much, but I was trying to figure it out from the definition

[dense glen]

Oh yeah yeah

I looked up the proof because I hate analysis

[dense glen]

Next question is on random walks and I cannot read my notes from years ago when I first was working on this

[wraith wing]

I think this fits here, please correct me if I am wrong.
I started reading a chapter on Finite Element Method and the author gave this equation but I can't wrap my head around how this came. Can anyone please explain that? Thanks

The book is 'Numerical Techniques in Electromagnetics with MATLAB' by M.N.O Sadiku

[grave spoke]

Read the sentence before the underlined equation. It's the most common choice for the approximation space (linear triangular finite element), where the solution is approximated by a linear polynomial within a triangular element.

A basis for linear polynomials in 2D is given by the constant 1 and the monomials x and y, so any linear polynomial can be written as a linear combination of the basis functions.

Although, usually, you use nodal basis functions instead that are transformed from a reference triangle.

[wraith wing]

hmm okay, I kinda find FEM to be more difficult to understand than FDM but I will keep reading and that might help me understand it better.

[grave spoke]

It requires some familiarity with weak formulations (and Sobolev spaces), but linear FEM is not that hard to grasp, good luck!

[wraith wing]

Thank you!

[hollow wraith]

can someone explain me the bolzano weierstrass theorem, or am i in the wrong channel

[grave spoke]

#calculus or perhaps #real-complex-analysis

[burnt pendant]

Hello , could someone help me understand how I'd be able to use Gram-Schmidt method to produce a polynomial in order to approach a function using the Least Square Method in specific points ? I am kind of in a mess after tons of hours of studying 😓

[drowsy idol]

not sure if this is the right channel but

like a year ago I had this as a homework

and I still can't figure it out

Let $f \in C^2([a, b], \mathbb{R})$, $a = x_0 < \ldots < x_n = b$ with $n \in \mathbb{N}$ being a subdivision of $[a, b]$, and $s$ an interpolated cubic spline to $f$ and the knots $x_0, \ldots, x_n$.

Prove the following inequality
[ \int_a^b (s''(x))^2dx \leq \int_a^b (f''(x))^2dx ]

for the cases \

1. $s$ is a natural spline; $s''(a) = s''(b) = 0$ \
2. $s$ is a complete spline; $s'(a) = f'(a)$ and $s'(b) = f'(b)$ \
3. $s$ is a periodic spline; $s^{(m)}(a) = s^{(m)}(b)$ for $m = 0, 1, 2$

[drowsy idol]

no idea how I would even begin proving this

[grave spoke]

variational calculus?

[drowsy idol]

?

[grave spoke]

Nvm I was thinking of another related proof. For your problem, you can write $$\begin{align*}0&\leq\int_a^b(f''(x)-s''(x))^2\dd{x}\&=\int_a^b f''(x)^2\dd{x}-2\int_a^b(f''(x)-s''(x))s''(x)\dd{x}-\int_a^b s''(x)^2\dd{x}\end{align*}$$ and try to show that the middle term vanishes.

[pine jetty]

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

[warm otter]

have you tried showing it for n = 1? Then s is just an interpolation polynomial of degree 3

[round spade]

Hi, the last number is supposed to be a composite number but in my code, it returns true. (i basically said if a number within 2 and n-1 divides into n, then it is not prime but clearly this doesn't work???)

[still forum]

just try is_prime(15) and see if anything has gone wrong there. You may be able to figure out what's going wrong from that.

[fathom rain]

none of the input arguments are divisible by two , so it returns True at the first iteration in the for-loop

[round spade]

what does ((n^1/2)+1) have anything to do with prime numbers?

[dense echo]

There's no computer ever built that has yet had time to count from 2 to 2^67-2 since it was put together. Since your program has already given you a result, it is not doing what you think it is.

[round spade]

it's grand now, i fixed that code

[fathom rain]

no need to check for more than that (but there are much better algorithms than this brute force exhaustive search)

[wanton plaza]

Hi, I'm trying to do a plot to show convergence for my algorithm that solves the following bvp using shooting method with derivatives represented using finite differences. But I don't have the exact solution. Atm I'm varying the grid spacing and plotting that against the 1-norm of the difference in solution between adjacent grid spacings (plotted) but that's probably wrong. I also considered replacing the exact solution with the solution using the smallest grid spacing.

[wide spear]

Usually what you do when you don't have the exact solution is just to compare against something that has a smaller grid spacing than everything else

So it looks like you go down to 2e-4 grid spacing

So you could compare everything against a grid spacing of say 1e-4 or 1e-5

Also usually you take the 2-norm difference

[wanton plaza]

Ah I see. Thanks a lot.

[dense glen]

this may be dumb question but is there an equivalent to power iteration for matrices over a finite field?

Actually nah lol

[gray tundra]

Guys do u have a recommendation for learning about fft

[tiny lake]

I'm studying De Boor's algorithm for finding the coefficients of a B-cubic spline, but I don't quite understand what e.g. x_{-3} would be defined as

We have n + 1 arguments as part of the input: x_0, ..., x_n (+ the values of the function corresponding to them ofc)

Every resource I've found disregards that issue completely, so maybe I don't see something

[grave spoke]

I assume that you have an open knot sequence, so $$x_{-k}=\dots=x_{-1}=x_0,\quad x_n=x_{n+1}=\dots=x_{n+k}.$$

[pine jetty]

zan #annoyedbirdemote

[tiny lake]

Wouldn't that result in the denominator being zero for $$V_{-3}^3$$

[pine jetty]

Thingoln

[tiny lake]

When k = 3 ofc

[grave spoke]

Yes, however, the convention here is to treat 0/0 as 0.

[swift mirage]

I am asked to calculate global error for x_n=0.16 n=2^k for (0,2,3...15) and h=x_n/2^k where h is the step size. However i get that for small h the Euler and Adam Bashforth method gives me a better result than runge kutta is that normal?

also i am asked to draw the log(h) against log(En) graph and i need to comment on the accuracy of the result and comment on the theoretical accuracy. doesnt theoretical accuracy depend on local error not global ones?

[swift mirage]

nevermind I'm dumb lol

[tulip bane]

Can someone explain me what is that algorithm? it is like Symplex Algorithm from Operations Research but you multiply the pivot column with -1, in my numerical analysis professor's slide it says Gauss-Jordan Elimination but i took linear algebra and im pretty much sure this is not Gauss-Jordan elimination

[tulip bane]

even it is gauss-jordan elimination i couldnt find anything online about it

this is the algorithm i think, but im just looking for the name of it

i posted here because the course that i take is numerical analysis, thx for any help

[tulip bane]

when u do the same operations to the last matrix u get this matrix

and it says the roots are these

i didnt post all slide because it is in a different language

[plucky kayak]

gauss-jordan is just a fancy way to say gaussian elimination to confuse students

[tulip bane]

I'm aware but this is not gaussian elimination right?

[plucky kayak]

doesn't look like it, no

[tulip bane]

i feel pretty dumb I'm trying to figure that out for 3 days

[hoary light]

this is indeed Gaussian elimination, but it's just storing the L factor in the lower triangular region simultaneously

that's why you're able to compute A inverse and x at the same time

I can explain further later, if this is not clear

[tulip bane]

ooh i get it now thank you, but i need a little more to connect everything together in my mind, to not to give any trouble to you is there a specific source that i can learn the connection

[hoary light]

you can look at a textbook such as Trefethen's "Numerical Linear Algebra" for the process of performing LU factorization. Basically, try doing Gaussian elimination on the problem that you posted, in the normal way. And then try doing LU factorization on that matrix A. You can see that you can actually do both at the same time, and the numbers at each step will match the numbers in the picture

[vale willow]

In general, when you want to test whether a series diverges or not, naive numerical experimentation won't really be helpful right?

[warm otter]

you can use a shanks transform to accelerate convergence if it converges, see the quote on the wikipedia page https://en.m.wikipedia.org/wiki/Shanks_transformation

[tulip bane]

thank you, you helped me a lot. I understood very well

[eternal flame]

Anyone that have pdf of this book:
Applied numerical methods with MATLAB for engineers and scientists fourth edition?

[fathom rain]

http://library.lol/main/10C3CF6DA33B96782A3264CC5FE3EC18

[eternal flame]

Tack! 😀

[lean cypress]

this is from page 4 of this paper: https://journals.sagepub.com/doi/pdf/10.1177/1756829317695564

how did they get the update rule and how is (3) the normal equation?

[neon cedar]

What is a Julia set? I've been told it's the complementary in a Riemann sphere of a Fatou set but I don't understand what that means.

[hoary light]

This isn't really numerical analysis, but did you already read the description on Wikipedia? The Fatou set is like the "stable" set, hence the Julia set is like the "unstable" complement

[neon cedar]

Thanks, and what is this for then?

[red jackal]

Mostly making pretty fractals

[neon cedar]

oh I meant the chat hahaha

[red jackal]

Numerical methods for analysis...scroll up

[wide spear]

Numerical analysis is not numerical methods for analysis

[ornate sierra]

Hey guys, any help with this one please?

[raven dew]

What is x_{k+1}? It isn't defined anywhere

[ornate sierra]

Yeah sorry, 1 sec

I was trying to play with it around but no luck so far

[raven dew]

I would just expand f(x_{k+1}) - f(x_k), since we know we should get lots of calculations and wind up with just one term left

[ornate sierra]

Yeah was trying that earlier but no luck, let me write that on the computer show u what I got

[ornate sierra]

Had a typo

[ornate sierra]

never mind 🙂 figured it out

[broken halo]

Let N be a numerical integrator (for example using RK45), $t_0$ a start time, $t_e$ an end time.
Let there be a ODE $y^{\prime} = f(t,y)$ that is linear but has no constant coefficient the solution would be over the complex domain.
Let there be $t_1, t_2$ such that $t_0 < t_1<t_2 <t_e$. And $N(f(t,y), t_0, t_e) = y(t_e)$ therefore applying the integrator results in obtaining the value of the function at $y(t_e)$.
Is this true?:
$$N(f(t,y),t_0,t_e) = N(f(t,y),t_1,t_e) = N(f(t,y),t_2,t_e)$

So therefore the history of the integration would be irrelevant? I do not think so but i havent taken any formal education on numerics and just thought about it when thinking about parallelization of an algorithm

[pine jetty]

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

[hoary light]

It's not true, because numerical integrators are not exact except for certain types of function f

In other words, N(...) does not equal y(t_e) exactly, it is only an approximation in general

As a consequence, the starting time does matter. For a general function f, the further away the start and end times are, the more error will accrue.

[daring shuttle]

Hey all, I'm still new so I hope this is a right place to ask about this.
I'm currently working on a C programming assignment, where I got caught by surprise by the necessity to implement Laguerre polynomials in a program (this is my first semester and I haven't worked with numerical analysis yet). Specifically, I need to write functions that calculate a Laguerre polynomial of nth degree for x, as well as the first, second and third derivative of one.

I've found recursive formulas for both the base polynomial and the first derivative on Wikipedia and I'm definitely going to stick with them, I just need help with the formulas for the second and third derivative.

[upper cape]

So I have to do a multidimensional quadrature over the unit triangle $T$ using the Duffy-Transformation:
$$[0,1]^2\rightarrow T , \ (x,y)\mapsto (x,(1-x)y)$$

I already have the gauss legendre weights and points for integration over [0,1]:
$$\int_0^1f(x)dx=\sum \frac{1}{2}w_i f(\frac{1}{2}x_i+\frac{1}{2})$$
Since I never did multidimensonal quadrature I was wondering if it is just:
$$\int_Tf(x,y)d(x,y)=\int_0^1 \int_0^1 f(x,(1-x)y)(1-x)dxdy$$
$$=\sum_i \sum_j \frac{1}{2} w_i \frac{1}{2} w_j f\left(\frac{1}{2}x_i+\frac{1}{2},(\frac{1}{2}-\frac{1}{2}x_i)(\frac{1}{2}y_j+\frac{1}{2})\right)(\frac{1}{2}-\frac{1}{2}x_i)$$

[pine jetty]

Enoo58

[hoary light]

What you're doing looks correct; this is called a product quadrature, where you treat the double integral as iterating two integrals, integrating in each of the independent variables x and y, in which case you can apply the 1D quadrature in each direction. This only works if your final domain is a square, or more generally, the set product of two 1D domains.

[wide spear]

They're integrating on a triangle though

[hoary light]

They were doing a change of variables via the map they gave above but I admit I only glanced at it, didn't check if they did the Jacobian part correctly

[crimson hamlet]

Any resources anyone can recommend on numerical integration wrt the haar measure for some specific Lie groups? I'm particularly interested in locally compact affine groups, i.e. semi-direct products like R^2 ⋉ SL(2,R) and R^2 ⋉ GL^+(2,R).

[wide spear]

Weyl integration formula may be helpful

[spark summit]

5x10

pls

/help

[crimson hamlet]

doesn't this apply only to compact groups? unfortunately I'm working only with local compactness

[wide spear]

Oh unfortunate

[warm otter]

im not sure about how to handle the semidirect products, but i think for the lie groups SL(2, R) and GL^+(2, R) you can just use integration with respect to the lebesgue measure on R^2 because it seems like the haar measures of SL(2, R) and GL^+(2, R) are absolutely continuous wrt to the lebesgue measure

[crimson hamlet]

I am not a mathematician so I'm sorry if my reply doesn't make sense. I assume you mean breaking up the integral into something like this? Do I not lose left/right invariance of the haar measure in this case? I should have clarified from the start this is required.

[bold yew]

"dd" is terrifying notation here lol

need some \mathrm{d} in there

[warm otter]

i mean the haar measure is the unique measure that has this invariance property and the right hand side you obtain by a change of variable giving you d\mu_{GL_2(\R)} in terms of the usual lebesgue measure

[upper cape]

I am learning about FEM and I got this excercise which shoudlnt be to hard but I am stuck:

I have a non-degenerative Simplex $T\subset \mathbb{R}^d$ with corner points $a_0,...,a_d\in \mathbb{R}^d$. Show that there exists constants $C \geq c > 0$ independent of T, such that for
$$u\in\mathcal{P}1:={ p:\mathbb{R}^d\rightarrow \mathbb{R} \mid p(x)=\sum{|\alpha|\leq 1}\lambda_\alpha x^\alpha , \ c_\alpha \in \mathbb{R}} $$

$$c \varrho_T\sum_{i=0}^d|u(a_i)|^2\leq | | u| |{L^2(T)}\leq C h_T \sum{i=0}^d|u(a_i)|^2$$

$$| |\nabla u||{L^2(T;\mathbb{R})}\leq C \varrho^{-1}| |u| |{L^2(T)}$$

$\varrho$ is largest radius of inner circle of T
$h$ is diameter of T

[pine jetty]

Enoo58

Enoo58

[hoary light]

Maybe to start, can you prove it for a standard simplex?

[upper cape]

I tried that but not that successful. My strategy was Proofing for the standard simplex and them maybe extending it with the Map to any other simplex. My main point of trouble is getting from the inequality for the determinant to the one with the L2 norm

[hoary light]

So you can at least prove the first step?

[upper cape]

not really 😦

[hoary light]

How about this, can you prove it for d=1?

[upper cape]

wouldnt that be:
$$\rho_T(|u(a)|^2+|u(b)|^2)\leq \int_a^bu(x)^2dx\leq h_T(|u(a)|^2+|u(b)|^2)$$
and $\rho_T=h_T=b-a$

[pine jetty]

Enoo58

[upper cape]

i think $c=1$ works and for C I need something larger

[pine jetty]

Enoo58

[upper cape]

oh and since we are in 1-d $u=\lambda_1 x+\lambda$
and then $$\int_a^b (\lambda_1 x+\lambda)^2dx=\frac{(ab+b)^3-(a^2+b)^3}{3a}$$ with that I should be able to find a C I dont think i need to be more specific

[pine jetty]

Enoo58

[hoary light]

I apologize for being blunt but this calculation is complete nonsense

Furthermore you'll need to justify why you think c=1 works (I don't think it does)

Anyways my general suggestion is to try to unravel the 1D case and see how you might generalize that argument. Focus on the standard simplex whose vertices are the origin and unit coordinate vectors.

[upper cape]

wait can you explain whats wrong? Did i integrate incorrectly or is it even wrong that i think u=\lambda_1 x + \lambda. I thought since we are in 1-d that would be correct

[hoary light]

Your final integrated value has no lambdas in it

[upper cape]

lamo yeah I didnt even realize sorry, I just plugged it in wolframalpha, cause i was to lazy. How would I generalize that to d-dimensions tho

[hoary light]

You're really concerned about the generalization step but in my opinion you haven't done the d=1 step correctly. Once you do that I think it'll be slightly clearer. In particular you might want to use a Cauchy-Schwarz type inequality

[upper cape]

so if i do this on the [0,1] I get with $u=\lambda_1x+\lambda_0$:
$$\lambda_0^2+(\lambda_1+\lambda_0)^2\leq \frac{\lambda_1^2}{3}+\lambda_0^2+\lambda_1\lambda_0\leq \lambda_0+(\lambda_1+\lambda_0)^2$$
because $h_T=\varrho_T=1$ With this I just have to put in c,C and then rearrange

[pine jetty]

Enoo58

[upper cape]

I think I went to deep into the 1-d case. I didnt even need to use Cauchy-Schwartz

[hoary light]

What you wrote above fails to hold for lambda_0 = 0, lambda_1 = 1, so there's a mistake somewhere

[upper cape]

Yeah its missing the c,C constants sorry I should make that more clear

[hoary light]

Honestly I'm still not convinced you understand the 1D case. I'm sorry maybe I did not indicate a strong enough direction to work towards here, let me write out what I meant

$\int_0^1 (u_0 + u_1 x - u_0 x)^2 \leq 3 \int_0^1 u_0^2 + u_1^2 x^2 + u_0^2 x^2 \leq 4 (u_0^2 + u_1^2)$

[pine jetty]

TheMipchunk

[upper cape]

Sorry I am not sure I am really following you here.

[hoary light]

In 1D the integral on the left is the one you care about. I compute an upper bound using an inequality similar to (a+b)^2 \leq 2(a^2 + b^2), and this allows me to provide an upper bound in terms of the sum of squares of u at the endpoints, which was what was asked in your problem

This is not the only way to do it, but you need to express your integral in terms of the endpoints values of the function u and in every calculation you showed above so far, you didn't do that

Furthermore, your calculation should, at least in principle, be independent of the interval of integration

[vapid plume]

Suppose that $y = y_h + \mathcal{O} (\beta(h))$ and $z = z_h + \mathcal{O} (\beta(h))$, for h sufficiently small. Does it follow that $y - z = y_h - z_h$ (for h sufficiently small)?

[pine jetty]

Sup?

[vapid plume]

I'm thinking that the errors are constant multiples of B(h) but not necessarily equal and so it is possible that y - z = y_h - y_z + O(B(h))

Another one is that show $\sum_{k=0}^{n} r^{k} = \frac{1}{1-r} + \mathcal{O} (r^{n+1})$.

[pine jetty]

Sup?

[vapid plume]

We know $\sum_{k=0}^{n} r^{k} = \frac{1-r^{n+1}}{1-r}$ and we can rearrange terms so that $\abs{\sum_{k=0}^{n} r^{k} - \frac{1}{1-r}} = \abs{\frac{-r^{n+1}}{1-r}}$

[pine jetty]

Sup?

[vapid plume]

Having trouble in finding an upper bound for |1/(1-r)| though

[wide spear]

Seems reasonable

[wide spear]

You need conditions on r

[junior swift]

I am stuck on 2 very 'basic' exercises. I am definitely overthinking them. I know I am supposed to start with f(x)=x and I gotta solve for x^2+c=x....

[cursive kite]

are you familiar with the quadratic formula?

[junior swift]

Yeah, I am... is that all I need to use? Cause I can't believe it is that simple

[cursive kite]

…

if you prefer to solve x^2 + c = x another way, by all means do so.

[dense echo]

I don't think finding the fixpoint(s) is all there is to the task -- especially part 1. It is a beginning.

[cursive kite]

but it was the only specific thing they mentioned.

[dense echo]

Hmm, right.

[quartz monolith]

Can anyone help me with this, I'm getting weird values. I have pretty much no idea about difference equations, I was just following the process my teacher told me to follow. I can't seem to get the answer

[hoary light]

This is an inhomogeneous system. Do you know the procedure for how to deal with the inhomogeneous term?

[quartz monolith]

I have to find particular solution is that right?

[hoary light]

You do, yes. Because your inhomogeneous term is 2^n (that is, some type of exponential of n), there is a very specific way to find this particular solution, did you learn it?

[quartz monolith]

I know the process, there's no problem with that, but the values for t I'm getting doesn't match the answer given in the book. Did I do something wrong?

[hoary light]

I see nothing wrong so far. But you should work out your full solution and then check if it satisfies the equation.

[dusty garnet]

Anyone know metric spaces?? is d(x,y) = 0 a metric space??

is it not a metric space because it should be d(x,x) or d(y,y) as for d(x,y) = 0, x =y

[inner briar]

Are you asking if setting d(x,y)=0 for all x and y defines a metric?

if so then no

[dusty garnet]

x,y are elements of X

[inner briar]

but either way this isn't the right channel

[dusty garnet]

there aint no metric spaces channel

sorry boss

[inner briar]

#point-set-topology

[dusty garnet]

shhhhit went right over my head

thnx man

[tepid river]

Anyone here ?

[wide spear]

If you have a question, just ask

[brave crypt]

is this the right channel for this?

[wide spear]

Sure

Are you intended to code this up?

It looks like you are

[brave crypt]

we can either code it up or give a simplified math version of it

im doing it by hand

[wide spear]

Ok do you know how fast fixed point iterations converge?

[brave crypt]

nope

i should show you the entire quesiton because

i did the first part

[covert ivy]

Hi! Does anyone have experience with the bilinear transform as a numerical method? I am familiar with classical numerical differential equations and quadratures, but have not encountered much on control systems.

In particular, I am looking at a the equation

$$x(t) = Ax(t) + Bu(t)$$
(coming from a state-space model) and the way that I would handle this is by multiplying by the matrix exponential and integrating to get:

$$x(t) = e^{At}x(0) + \int_0^t e^{A(t - \tau)}Bu(\tau) ,dtau $$

and then handling the convolution integral on the right hand side using some standard quadrature. However, the authors of a paper I'm looking at use the bilinear transform.

My understanding of the bilinear transform is that it computes the transfer function of the trapezoid rule and the transfer function of the continuous integral and equates them to give a transformation to apply in "Laplace transform space" to go from continuous to discrete time. Is this correct?

Now, If I take the bilinear transform as a given and apply a Laplace transform to my system, I understand that after a bit of algebra I will arrive at the discrete Laplace transform (or z-transform?) of the system:
$$x(t + h) = (I - h/2 A)^{-1}(I + h/2 A)x(t) + h(I - h/2 A)^{-1}Bu(t)$$
which I am not sure how to derive from the original system using a traditional numerical method (I tried just applying trapezoid rule + pade approximation but don't see immediately how to get the $$h(I - h/2 A)^{-1}B$$ portion).

My questions are:

How does this compare to other numerical methods? In particular, what is the order of this approach and what would a be a comparable classical numerical method? What are the benefits/drawbacks of this as a numerical method? Is there a way to derive this as a special case of a classical quadrature using the convolution integral as above over each timestep?

[pine jetty]

Probably_Jason

[wide spear]

Are you missing some time derivatives

[covert ivy]

Ah, yes! the differential equation should be
$$\dot{x}(t) = Ax(t) + Bu(t)$$

[pine jetty]

Probably_Jason

[wide spear]

At any rate, something that jumps out to me immediately is that you went from a method that requires no linear solves to a method that does

Do you know anything about the matrix I-h/2A?

[covert ivy]

Yes, but really the linear solves are super minor. Even with an adaptive timegrid (with some extra work) there are ways to precompute $$(I - h/2A)^{-1}$$ and cache it for many applications

[pine jetty]

Probably_Jason

[wide spear]

Is it well conditioned though

[covert ivy]

The first bit is the pade approx $$(I - h/2 A)^{-1}(I + h/2 A) \approx e^{Ah}$$ which is understandable, it's that second one in that I do not understand

[pine jetty]

Probably_Jason

[covert ivy]

In my case, yes

for what it's worth, most quadrature methods involving the matrix exponential in the convolution integral will need to compute the same pade approximation

So caching it is the best option

[wide spear]

Yes ok so you can do something with say a super LU

Ok that isn't a concern then

Have you tested this out numerically?

[covert ivy]

yes, it's numerically fine

I forgot to mention a big and important bit!

A is diagonal

[wide spear]

You can numerically check the order then right?

Oh A being diagonal is such a scam

Unfair, even

Why can't my matrices be diagonal

[covert ivy]

Yeah I suppose that's the best way. I believe there are some nice reasons for applying this transform in particular but I just can't find a comparison to the methods I'm more familiar with

[covert ivy]

follow-up question. Is the discretization of the system via the bilinear transform only viable for a fixed grid? Ie. would it possible to apply it to a general partition like the numerical method?

[wide spear]

You mean for a general quadrature grid right

I have no clue

[covert ivy]

yep. Or even just dyadic

[wide spear]

I feel like it should be possible, the update will just be a lot messier

[covert ivy]

hm :/ For the quadrature of course it is dead simple

[warm otter]

id like to solve this problem numerically, is there a standard way to somehow transform the limit condition to a boundary condition?

[wide spear]

So this is already a boundary condition

It’s called a radiation boundary condition

I’m sure there are ways of handling it

But you’ll be working on a finite spatial domain right

You might be able to rewrite this as the derivative of phi at the right endpoint is 0

[warm otter]

i thought about approximating it by applying a 0 dirichlet boundary condition at the right endpoint

is there a justification for applying it as a neumann condition?

[wide spear]

No I don't think a 0 dirichlet boundary condition is the right way to do it

For example, f(x)=x(1-x) satisfies f=0 at x=0, 1, but lim f(x) as x->infinity is certainly not 0 right

Actually maybe the neumann boundary condition isn't good either

What if you did a transform

Like x=arctan(z)

That turns infinity into pi/2 so you do have a dirichlet boundary condition there

[warm otter]

ill try this, but i have the feeling that this will make the equation a lot nastier 😅

[hoary light]

there needs to be some restrictions on a(z). I assume it is always positive? When a(z) is like a positive number, then you can observe the solution will be exponentially decaying, whereas if a(z) is like a negative number, then the solution is like a sinusoid

[warm otter]

i want to solve the equation for some a(z) which look like z^r for r in (-infty, 1) or linear combinations and see how phi'(0) behaves

[shy thistle]

solve the following system of equations;
2X1 + 2X2 − X3 + 3X4 = 10
X1 − 3X2 + 2X3 + 2X4 = 20
−X1 + 4X2 + 10X3 − 2X4 = −5
3X1 + 4X2 − X3 + 10X4 = 12
Using

1. Naive Gaussian Elimination
2. Partial Pivoting
3. LU decomposition

[wide spear]

So what do you need help with?

What have you tried?

[autumn mica]

can someone help me with this. In solving the problem

[wide spear]

#calculus

[vagrant tide]

Is numerical analysis and numerical methods the same thing? I'm in a numerical methods class and was just curious

[wide spear]

The first non-parenthetical words of the channel description are "Numerical methods"

lol

[vagrant tide]

Yeah I probably should've read that lol I see it now thanks

This may be a dumb question but I'm still new to numerical methods. But if you have a set of 4 data points how would you go about finding how many polynomials of degree 2 pass though all the given points?

[wide spear]

Set up a linear system, find the solutions

[vagrant tide]

Alright, that's kind of what I was thinking because someone else in my class did it that way but I wasn't sure if that was the correct way

Thanks

[simple vigil]

hmm I am wondering if something like what I am looking for already exists

essentially, I am looking to have a 3d "bowl" with decently sharp walls

Actually NVM i think i can do this w some easy manipulations to a gaussian

[brave crypt]

Hello, could someone please elaborate on what the "K" in the definition means? It is from the book "Numerical Mathematics" by Quateroni. I believe that it wasn't defined and yet it is being used.

[plucky kayak]

I'd guess some function that depends on eta and data? in other words the change in the solution dx is no bigger than some number K times the change in data

and eta in that sense would have an interpretation of some sort of perturbation quantity since the change in data is smaller than eta

[brave crypt]

Thanks a lot. This makes sense!

[unreal mural]

it is just an upper bound that depends on eta and d

[stiff dawn]

true or false?

[wide spear]

Ah my original statement may have been a bit dramatic

There are no good introductory numerical analysis books

This is because introductory numerical analysis is more of a collection of separate techniques

[twin sonnet]

Bro was out for blood for a message from a year ago 💀

[wide spear]

Lol

No one has disagreed yet

[stiff dawn]

i just lurk.

[wide spear]

https://arxiv.org/abs/2302.11474

May be of interest to some people

[fossil badger]

someone wants to go through the book "Numerical Analysis" of Burden and Fairs together?

[quartz monolith]

Can anyone tell me how to find particular solution for this kind of problems?

[hoary light]

Do you know how to obtain the particular solution if the RHS were just one of the two terms?

[leaden parrot]

im trying to figure out how to get to the RHS. here $t \in [0,1]$, $l>0$ and if for some a,b such that $||a-b|| \leq l$ then $(||a-b||-l)^2$ vanishes (=0)

[pine jetty]

Eso

[leaden parrot]

basically trying to show its convex

[deft venture]

Is it fine to ask a question about the taylor remainder of a taylor polynomial here?

[dense echo]

If you want to compute/estimate it in practice, yes.

[deft venture]

Yes I want to use it later on in combination with the error function and MATLAB

Since it's about the function e^(-x²)

[wide spear]

Well what’s the question

[deft venture]

How would I go about determining the Taylor remainder of the nth degree Taylor polynomial of the function above?

[wide spear]

You can find a formula on Wikipedia for the remainder

[deft venture]

Actually?

I did try googling first but wasn't really able to find anything on the remainder of e^(-x²). Perhaps I missed it haha

Mind linking me perhaps?

[wide spear]

There will be one for general functions

You can then plug in your specific function

[deft venture]

Right but how would I know what the nth derivative of e^(-x²) is?

Since it does not seem to be something simple

[wide spear]

You can find formulas for it online as well

Alternatively you could try to derive a formula for it

[deft venture]

I actually did try looking for it but did not find it. I did try deriving it myself too and got something, but was unable to figure out the coefficients for each term sadly.

[wide spear]

You may look at the wikipedia article on hermite polynomials

[deft venture]

Oki I will look there!

Awesome I think I kinda get it from there. Thank you!

[obtuse furnace]

Hi,
What are some of the big problems/challenges that people in numerical analysis are currently trying to solve?
This is a broad and kind of vague question, so I'd like to hear about anything that comes to mind for anyone.

[fathom rain]

I do efficient methods for high precision computations of eigenvalues and eigenvectors of structured matrices (especially important for non-Hermitian matrices where typically normal methods fail)

[wide spear]

Lots of stuff on scaling various algorithms to supercomputers effectively

Lots of stuff in various domain areas

Randomized numerical linear algebra has been hot recently

Reduced/mixed precision algorithms

[obtuse furnace]

Are you referring to problems that arise specifically with numerical methods in certain applications?

[wide spear]

Yes, improving techniques for specific problems

[carmine swan]

Is there any standard way to find some polynomial having some conditions for eg. I want to find a polynomial p(x) passing through origin such that p(x)-1 has m positive zeroes and p(x)+1 has n positive zeroes for arbitrary m,n ??

[wide spear]

DId you read my response

[carmine swan]

Oh sorry.

[wide spear]

Did my response make sense

[carmine swan]

I am sorry coz was unable to reply .Yeah your response make sense . I am still trying that problem: a polynomial p(x) passing through origin such that p(x)-1 has m positive zeroes and p(x)+1 has n positive zeroes for arbitrary m,n ??

Can you please help me to construct such .

[mossy prawn]

can someone explain the QR-Decomposition?

[wide spear]

Pick m points x_i so that p(x_i)=1 and n points y_i so that p(y_i)=-1

You then have m+n equations to solve for a degree m+n polynomial

[wide spear]

What questions do you have?

[carmine swan]

But that can be problematic. Because they may have repeated root at one point.

[wide spear]

This polynomial will not have any roots, a degree d polynomial is uniquely determined by d+1 points

https://link.springer.com/book/10.1007/978-3-031-08121-7

Also new book

May be of interest to some people

[carmine swan]

Yeah but it may happen that there is another z rather than x_i and y_i's such that f(z)=1 or f(z)=-1

That would be a problem right?

[wide spear]

Ok new idea

Built the polynomial based on how many times it crosses 0

Oh no nevermind

Hmmm

[carmine swan]

Actually I think my claim is false for both m,n odd / even.

[hazy jolt]

[hazy jolt]

Hi guys, could you give me an idea for the solution of this exercise?

On a binary machine, what is the unit rounding error for 48-bit mantissae?

[ionic karma]

heyy...any book on numerical methods ?..any journal or paper would be better...

[wide spear]

There are many books

I linked a book if you scroll up a bit

There are also many journals

And many many papers

[ionic karma]

ahh this one nice...thanksss

lol sometime i fell 'numerical analysis' or 'numerical method ' is a weird name ..i mean word 'numerical' has a broad meaning..

[hazy jolt]

Guys could someone give me an idea with this exercise please

[dense echo]

What is your problem with it? Do you know what the exercise means by "unit rounding error"?

[hazy jolt]

I give up, I closed all the pages, I need help, can I clarify my doubts and answer my questions with you ?

[hazy jolt]

can help me ? :c

[dense echo]

Not without knowing what "unit rounding error" means ...

[hazy jolt]

look

Maybe it is not the same as what you know but that is what I have, I have been reading the book that my professor uses and it gives as a result this procedure then we have this

More reduced this 😄

[hazy jolt]

[inner briar]

the rounding error depends on the value you put into fl

are you looking for a bound on that error?

[hazy jolt]

On a binary machine, what is the unit rounding error for 48-bit mantissae?

I was talking to a friend of mine and we came to the same conclusion.

If a MARC-32 does not correctly round numbers but simply removes excess bits, what would be the rounding unit?

[hazy jolt]

Then I have this question that is complicated to do I have the solution to this exercise but I have no idea.

I don't know where that 2^1 and fl(x+E)>1 come from.

[hazy jolt]

[minor osprey]

2^1 is because (0.100…)_2= 1/2 in base 10 and fl(1+ eps)>1 is the floating point result to the input 1

[ionic karma]

hey..is is possible to find all the roots of given function using numerical method?

[wide spear]

If a function has infinitely many roots, no

If a function is discontinuous, no

If your function is continuous, you can use a bisection method

[ionic karma]

so what's the use of numerical methods when there are more good and effective methods?

[wide spear]

What

The bisection method is a numerical method

[ionic karma]

nevermind my bad..i realised i asked bad question..lol

[hazy jolt]

I see that no one likes numerical analysis.

[wide spear]

?

[raven dew]

Except for those that do

[amber slate]

So I need a help with the numerical solution of PDE(s). The differential equation (two-dimensional) that I'm trying to solve is [\nabla^2\phi=1] I know that I could use the Finite Difference method for this, and I did, but the problem is that the boundary condition is the Neumann boundary condition instead of the Dirichlet boundary condition, where in ALL bounds, [\pdv{\phi}{x}=0] and [\pdv{\phi}{y}=0] I tried to calculate this using a 4x4 two-dimensional grid, and built the coefficient matrix to solve using Gauss Elimination. But the problem is, I can't seem to get the pattern of this matrix, and I want to know how (the problem isn't actually 4x4, it's actually 201x201, so I'm trying to find the pattern so when coding it's easier for me to construct the coefficient matrix). I'll attach the matrix below.

[fathom rain]

create the toeplitz with -1 2 -1 (2 on the main diagonal) but change the topleft and bottomright corners to 1. call matrix L then do kron(L,In)+kron(In,L) assuming a box and n is same in both dimensions and In is identity of size n

the changing to 1 makes it Neumann

[amber slate]

will do

thanks

[fathom rain]

and you should instead use 1 -2 1 matrix for your problem i guess

[amber slate]

the a b c comes from this:

$$\pdv[2]{\phi}{x}+\pdv[2]{\phi}{y}=1$$
Using the finite difference approach, where,
$$\pdv[2]{\phi_{i,j}}{x}=\frac{\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}}{h_x^2}$$ and $$\pdv[2]{\phi_{i,j}}{y}=\frac{\phi_{i,j+1}-2\phi_{i,j}+\phi_{i,j-1}}{h_y^2}$$
Then,
$$\frac{\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}}{h_x^2}+\frac{\phi_{i,j+1}-2\phi_{i,j}+\phi_{i,j-1}}{h_y^2}=1$$
Or,
$$\frac{\phi_{i+1,j}}{h_x^2}+\frac{\phi_{i-1,j}}{h_x^2}+\left(\frac{-2}{h_x^2}+\frac{-2}{h_y^2}\right)\phi_{i,j}+\frac{\phi_{i,j+1}}{h_y^2}+\frac{\phi_{i,j-1}}{h_y^2}=1$$
\begin{equation*}
\begin{split}
a&=\frac{1}{h_x^2}\
b&=\left(\frac{-2}{h_x^2}+\frac{-2}{h_y^2}\right)\
c&=\frac{1}{h_y^2}
\end{split}
\end{equation*}
Hence,
$$a\phi_{i+1,j}+a\phi_{i-1,j}+b\phi_{i,j}+c\phi_{i,j+1}+c\phi_{i,j-1}=1$$

[pine jetty]

calicals

[fathom rain]

a = 1 b=-2 c =1 (and put the 1/h^2 in front of the created L

literallt create the 1d problem and then just use kron

everybody makes this so complicated 😛

[amber slate]

that's true

[fathom rain]

i can do the code for you in a few lines in julia

[amber slate]

that would be helpful!

[fathom rain]

w8

[amber slate]

alrighty, danke

[fathom rain]

using LinearAlgebra

function laplace_1d_dirichlet(n)
    return diagm(0=>-2*ones(Int64,n),1=>ones(Int64,n-1),-1=>ones(Int64,n-1))
end
function laplace_1d_neumann(n)
    L=laplace_1d_dirichlet(n)
    L[1,1]=-1
    L[n,n]=-1
    return L
end
function laplace_2d_neumann(n1,n2)
    L1=laplace_1d_neumann(n1)
    L1=L1/(n1+1)^2
    L2=laplace_1d_neumann(n2)
    L2=L2/(n2+1)^2
    return kron(L1,I(n2))+kron(I(n1),L2)
end

should be correct

[amber slate]

that's very helpful, thank you!!

[fathom rain]

(there are multiple ways to impose boundary conditions so you might want it slightly different)

[amber slate]

will definitely take that into consideration

but all of this is very helpful

time to work

[vapid plume]

For the Durand-Kerner method to find the roots of a polynomial, once we have the region where we know where the roots are supposed to be, can we take complex numbers randomly to be our initial guesses?

Example from the book: For the polynomial x^4 - x - 1, we get the region |z| <= 2 and the initial guesses for the roots were taken to be: 1, -1, i, -i. My question is, could we have taken the initial guesses randomly from the region |z| <= 2 and still get convergence?

[austere moon]

I am struggling on this exercise. I am dont see the form that takes a symmetric quadrature. Can someone give an hint please ?

[wide spear]

Ok, let's start with a simpler problem

What is an example of a symmetric quadrature rule

[austere moon]

We can choose let's say $Q[f] = 2f(\frac{2}{3}) +2f(\frac{1}{3})$

[pine jetty]

Julien J

[wide spear]

What polynomials is this scheme exact for

[ionic karma]

heyy..for this equation x^5 + x + 1 = 0 is there way to find a root ..by the method of how we solve quadratic or cubic equation ??

i mean not just by trial or error method or any iterations

can we solve this symbolically without using numbers ???

[wide spear]

No, and this would be better suited for #groups-rings-fields

Without using numbers = not numerical analysis

[ionic karma]

fine thanks..

[austere moon]

I don't understand the question

[wide spear]

Have you seen, for example, gaussian quadrature before?

Or you might have seen it called gauss-legendre quadrature

[austere moon]

Just seen it very recently

So I am not very familiar with it

[wide spear]

Ok you know how gauss legendre quadrature is exact for polynomials up to a certain degree right

[fathom rain]

You can see that it is same as $(x^2+x+1)(x^3-x^2+1)=0$ and solve those two sub problems. Again, wrong channel

[pine jetty]

Sven-Erik

[austere moon]

2n+2 right ?

[ionic karma]

heyyy

how close the number should be in order to be approximate of exact value

??

[amber slate]

btw got a 9/10 on that assigment thanks

also got a new fancy thing to say: toeplitz

never heard that before tbh

[wide spear]

Toeplitz

Circulant

Hankel

[amber slate]

damn

am i that dumb in math

[wide spear]

Just different types of matrices

[amber slate]

but also im a physics undergrad

didnt take math because i thought i could never do it

turns out i have to dive deep down into linear algebra for QM

[fathom rain]

Toeplitz matrices are very common in physics. All PDE discretizations (on structured grids) are Toeplitz

(or Toeplitz-like)

[random canopy]

Dare I say vandermonde

[hazy jolt]

guys a question

[wide spear]

Well what's your question?

[hazy jolt]

Show that 4/5 cannot be represented exactly in MARC-32.
What is the closest machine number?
What is the relative rounding error that occurs when this number is stored in MARC-32?

It is this exercise

4/5 cannot be represented on the MARC-32 machine.
since its binary number base two gives 0.110011001100110011001100....

now for the second point tell me which is the closest number to the machine, the closest would be 0.11001100110011001100110011001100110011001100 base 2 which would be the closest number, so I am truncating the number 4/5 with respect to the MARC-32 machine.
0.110011001100110011001100 =0.79999995231628
is this correct?

[hazy jolt]

[hazy jolt]

help me pls

[dense echo]

I don't know "MARC-32", but wouln't you get a closer approximation it you round up at the least significant bit instead of truncating?

[hazy jolt]

Could someone guide me in this exercise?

[grave spoke]

#computing-software

[modest saffron]

I have a problem with terminology in finite volume method (though it applies also to other PDE discretizations). How to call the numerical error from the difference between the initial value u0 and interpretation of numerical state after u0 is converted to the discrete state? In my case that'd be the difference between u0 and piecewise-constant averages of u0 in cells. I do not think this is a round-off error as it does not have anything to do with floating point accuracy. But it also does not seem right to include this into the global truncation error as that one is an aggregation of local truncation errors.

[hoary light]

Perhaps spatial discretization error, as opposed to error accrued from temporal discretization.

[wide spear]

I agree with spatial discretization error

[modest saffron]

is there any article or book that mentions it? (i need to use very formal terms in my work)

[wide spear]

Yes lots of articles and books mention it

This is terminology I would use in a paper

[modest saffron]

thank you.

[leaden parrot]

I need help with showing $(|x|-k)^2$ is non-convex for some fixed k

[pine jetty]

Eso

[wide spear]

Well what have you tried

[leaden parrot]

well i tried so far starting from assuming its convex and then showing the supposed greater than or equal to side is less

i tried fixing lambda to 1/2

[wide spear]

Sorry what

[leaden parrot]

like

$\lambda (|x|-k)^2 + (1-\lambda)(|y|-k)^2 < (|\lambda x + (1-\lambda) y |-k)^2$

[pine jetty]

Eso

[wide spear]

Ah ok

[leaden parrot]

if this holds for for some lambda in (0,1) then we have proved its not convex

[wide spear]

Yes ok

So you do this

Are you successful?

[leaden parrot]

well apparentl not

thats why im asking for help

[wide spear]

Does this simplify

Also you can fix x and y as well

[leaden parrot]

i know but idk how to proceed

ive ended up with $\leq 1/2(|x|^2+|y|^2)-k|x+y|+k^2$

[pine jetty]

Eso

[wide spear]

Have you looked at a plot of the function

[leaden parrot]

the function is any dimension

[wide spear]

So?

[leaden parrot]

is there some algebraic manipulation im missing

i dont see how looking at a graph helps me

[wide spear]

Do you understand what that definition of convexity means

[leaden parrot]

oh you probably mean

give an example

[wide spear]

Graphically

[leaden parrot]

i hoped i could show this more "rigorously"

[wide spear]

You can, after developing some intuition about the function

[leaden parrot]

hmm i kinda see what you mean. at |x| = k its non convex am i right?

or rather its not differentiable

[wide spear]

It is differentiable there, but not convex

It is not differentiable at x=0

[leaden parrot]

right right. my bad

[wide spear]

And once you have this mental picture of the function, you then have a better idea of how to precede right

[leaden parrot]

right right

thx

[leaden parrot]

ok so assuming we fix some point x \in R^n. Then we can I think find a y such that ||y|| = k. however this y must be on some other side of the graph in a sense. This way there exists some lambda1 < lambda2 in (0,1) such that the line between these two lambdas are all below the graph

its simple in the 1 dimensional case but how do to project this y to assure its on the opposite side

[leaden parrot]

nevermind i figured this doesnt hold if you have only one free variable y

doesnt hold for any fixed x

[wide spear]

Indeed

[leaden parrot]

but ofc having both of them free and setting x = -y where ||x|| = k we have that for all lambdas between 0 and 1 that its non convex

[wide spear]

Yes

[marble temple]

My teacher has made a homework problem that looks very similar to one from a book (the picture) (Iterative Methods for Sparse Linear Systems, Second Edition). I'm struggling to find an answer to (ii) (although I have proven that it's converging through manually calculating the iterating matrix $M_{GS}$ for Gauss-Seidel) & (iii). Any ideas how I should go about with this one?

[wide spear]

Do you know anything about the convergence properties of Gauss Seidel? How did you do (i)?

[marble temple]

The Jacobi iterative matrix is pic, so I found out that similar matrices (Toeplitz tridiagonal matrices) have the eigenvalues defined as this pic

This is the Gauss Seidel iterative matrix. I know that if A is strictly diagonally dominant, irreducibly diagonally dominant then the spectral of spectral_rad(M_{GS}) < 1. Also if A is consistently ordered (haven't really understood that yet), then spectral_rad(M_{Jac}) = spectral_rad(M_{GS})^2

this is everything I know ab the convergence properties of Gauss-Seidel so far.

[wide spear]

Ok there is one more convergence criteria for gauss seidel

[marble temple]

hmmmm

interesting

[wide spear]

It is on wikipedia

[marble temple]

oo

are you referring to if A is SPD?

[wide spear]

Yes

[marble temple]

okay fair, I forgot to write that

[wide spear]

Your A is spd if alpha=2 isn't it

[marble temple]

but still my primary question was with what convergence rate it would be

yes exactly

[wide spear]

Hmmm

The convergence rate probably depends on the smallest eigenvalue

Intuitively, the idea is that if there is a 0 eigenvalue, you don't have guaranteed convergence anymore right

So the closer the smallest eigenvalue is to 0, the slower it converges

[marble temple]

From my understanding it's the opposite? Considering that the spectral radius checks for the largest eigenvalue

i.e. if the largest eigenvalue < 1 then it converges

[wide spear]

Oh wait yes nevermind

My bad

The internet seems to indicate that the convergence rate is proportional to the spectral radius

[marble temple]

yeah to be honest all of this is kinda confusing to me still

but yes I think the smaller maximum eigenvalue you have the better convergence

but do you know how I could extract the largest or smallest eigenvalue from a matrix like this?

[wide spear]

You have all the eigenvalues

[marble temple]

But only from the matrix A right? Or do they follow??

[wide spear]

Oh whoops

For that other matrix

[marble temple]

aah okay

[wide spear]

Is the upper triangle all 0

[marble temple]

yes exactly

[marble temple]

yes

this is the definition of the iteration matrix

D & L are just the diagonal and lower part of A

[wide spear]

Oh hmmm

So you can find the eigenvalues of (D+L)^-1 right

Well that doesn't help

Hmmm

[marble temple]

I guess but how would I do that?

okay

[wide spear]

The eigenvalues of a triangular matrix are simply the entries along the diagonal

The eigenvalues of A^-1 are the inverses of the eigenvalues of A

[marble temple]

oh what

okay I have completely missed that 😆

but does that say anything about the iteration matrix?

[wide spear]

No lol

[marble temple]

assuming I know the maximum eigenvalue of (D+L)^-1 and A, does that follow into M

aah fk

[fathom rain]

FYI the eigenvector matrix in that example is the discrete sine transform (typically denoted dst-1)

[marble temple]

Aah fellow scandinavian 🙂 that's cool didn't know that

[pine jetty]

Faputa

[warm otter]

does anyone here have experience switching from a pure mathematics master to a numerical analysis phd?

[fathom rain]

no, but I know many who have.

[warm otter]

so would you say thats quite common?

[fathom rain]

at least not uncommon

[pine jetty]

Roan_22

[faint flume]

does anyone know of a name for when using a higher order polynomial to approximate a dataset, lead to oscilations along the approximated polynomial, something like Runge's phenomenon

cuz it seems liek runges phenomenon is applicable when the oscilations are at the edges?

[wide spear]

Runge's phenomenon is fine

Not necessarily just at the edges

[grave spoke]

I thought #dynamical-systems is used less

[grave spoke]

Before the split into #numerical-analysis #modeling #optimization I guess...

[minor osprey]

Sparse matrices

[grave spoke]

[minor osprey]

Exactly!

[wide spear]

Congrats on mod kirby

[gloomy gorge]

actually only the eigenspace spanned by this eigenkirby and eigenzan is a moderator subspace. the other eigenusers are antimods, a.k.a. trolls.

[dusky lion]

am i a basis of the orthogonal complement of that subspace?

[plucky kayak]

already answered in #calculus

do you know what a jacobian is?

[pine jetty]

Transparent_Elemental

Transparent_Elemental

[pine jetty]

[daddy]Adika

[wide spear]

Try checking the definition of convex

[leaden parrot]

i posted this a while ago

it would likely mean id be done

[wide spear]

If this is about Eso’s question

Lol

[plucky kayak]

when given a function $f(\theta) = \frac{1}{2}||r(\theta)||_2^2 + \frac{\gamma}{2}||\theta||2^2$ would the levenberg-marquardt update formula for $\theta$ look like this? $\ \theta{k+1} = \theta_k - (J^TJ + \gamma I + \lambda_k I)^{-1} J^T r(\theta)$ ? here $J$ is the jacobian of $r(\theta)$ which has rows as gradients of $i$-th component of $r(\theta)$

[pine jetty]

Transparent_Elemental

[plucky kayak]

just want to make sure because everybody uses a different convention or derives the algorithm from a different problem and I can't tell if I'm understanding this correctly

or perhaps the inverse matrix has to be multiplied by the gradient of f? which is J^T * r + gamma * theta instead of just J^T * r

like in newton's method

[hoary prawn]

Suppose we have a function f(x1, x2, x3) = x1 + x2 + x3.

My lecture presentation claims that this is unstable when:
|x1| >= |x1 + x2 + x3|.

I tries this out in a sample C++ program with double fp and I could not find out how is that unstable.
For small changes in an input, I received linear changes in output.

How is that unstable?

[wide spear]

It’s saying that the function evaluation is unstable?

[hoary prawn]

it says that it is bad-conditioned, I believe

[wide spear]

Hmmm

Try playing around with some examples where x1 is much larger than the others while keeping in mind floating point

[hoary prawn]

auto unstable(double x1, double x2, double x3) -> double {
  return x1 + x2 + x3;
}

auto main() -> int {
  auto n1 = unstable(100, -2.499999, -2.4999998);
  auto n2 = unstable(100, -2.499998, -2.4999998);

  fmt::print("n1: {:f}\n", n1);
  fmt::print("n2: {:f}\n", n2);
}

$ g++ ./main6.cpp -std=c++20 -o ./a.out -lfmt && ./a.out 
n1: 95.000001
n2: 95.000002

[wide spear]

Ummmm

Much much larger

[hoary prawn]

Sure, this will perform worse probably in case of numbers that are larger due to loss of precision.
How can I determine, when such numbers are "bad enough" to produce actual errors?

[wide spear]

Are you familiar with the notion of machine epsilon

Lol

[fathom rain]

🙂

(take x1 to be 1 and the others of order 1e-16)

[hoary prawn]

I just learned about it :p.

So - error may be big, but it will be only seen if it is bigger than the machine epsilon?

and apparently in this case it was less than this epsilon

correct me if im wrong

[wide spear]

No not quite

I will let sven erik take this one

Because my flight is about to take off

[fathom rain]

what kind of topic is it about I find the formulation My lecture presentation claims that this is unstable when: |x1| >= |x1 + x2 + x3|.a bit weird

is it a about floating point errors?

(maybe I am missing something)

[hoary prawn]

it is about the propagation of error

[simple vigil]

so... cramming for finite element analysis midterm

one of the details pointed out in the notes that was not expanded on (realizing that professor is speedrunning a ton of aforementioned details) is that the system given by "second-order finite differences" (I recalled there being a difference between centered, left, and right? if someone knows which one? I assumed centered but couldn't get it to work out, though there is a nontrivial chance either side has some algebra mistake somewhere) coincides with the Kc = f using Courant basis for Poisson's equation on (0,1)^2

I can't find anything about this online... is someone familiar with this coincidence?

[wide spear]

Centered usually

What matrix do you get and what matrix do you want

[wide spear]

Ummm

[sacred lotus]

Here lies exaea's gpa 💀

[simple vigil]

flopped on responding but thanks -- talked to another student yesterday just before exam and it just amounted to screwing up the algebra

[hoary light]

There's only four variables but 9 equations. Are you trying to minimize the error in the least squares sense?

[jagged cipher]

Not sure if this is the right place for this. Let A be an n by n matrix and u an n dimensional vector. Let f(x) = (A^x * u)_0/x. I would like to compute the limit of f as x goes to infty, and I have some ideas about how approximations might help here.

Idea 1:
Approximate the entries of A^x * u using a closed form expression (nice enough to include in a limit). This approximation can be terrible as long as the error tends toward 0 with larger x.

Idea 2:
Compute f(x) for a large x and reason about how far away it is from lim x->infty f(x). This is less ideal, but might be good enough for my purposes.

Advice is appreciated

[hoary light]

what does the _0 part of your notation indicate?

nnz?

[jagged cipher]

v_0 denotes the first entry of v, sorry for not stating that. Also realizing f(x) = (A^x * u)_0 / x should instead be f(x) = (A^x * u)_0 / (B^x * v)_0 for a second fixed matrix vector pair, B and v

[river thicket]

Does anyone know why my code doesn't work?

n = size(x)
n = n(1)
vector_of_c_0_linear = zeros(n-1,1);
vector_of_c_1_linear = zeros(n-1,1);

for i = 1:(n-1)
v = [one_linear, x(i:i+1)]
c = v \ y(i:i+1)
vector_of_c_0_linear(i) = c(1) %c_0
vector_of_c_1_linear(i) = c(2) %c_1
end

Trying to find equation

((1,1),(x0,x1)) * (c0, c1)^T = (y1,y2)^T

$c{0} + c{1} x{1} = y{1}$

$c_{0} + c_{1} x_{1} = y_{1}$

[pine jetty]

Jonathan!

[wide spear]

We are not here to debug your code

[jagged cipher]

Maybe I can use the largest eigenvalues of A and B? Is it true that I can rewrite f(x) to be p/q where p = c1L1^x + c2L2^x +...+ ckLk^x, where L are the eigenvalues of A and c are magic constants that I don't yet understand? Similarly for q in terms of the eigenvalues of B. If yes, then the limit will be dominated by the exponential term(s) with the largest base(s).

[river thicket]

okay cool guy

nothing wrong with the code

just misunderstanding of something else haha

if anyone has a code that they are confused with in numerical analysis i am happy to help by the way

[sand heron]

Anyone knows a slightly nicer way to capture the effects of turbulence in boundary layers without letting off some accuracy behind such as in the case of Reynolds-averaged Navier–Stokes decomposition method. I tend to care more about spectral methods but a very accurate NS solution is my aimed goal. Thanks.

[wide spear]

Ah yes

A very common problem

Have you looked into

1. large eddy simulations

2. adaptive mesh refinement

[sand heron]

Uhm. There are some significant reservations here when it comes to Inlet boundary conditions and the predictions of LES. I'll slip more into the two and pick the one that superlatively fit the conditions I am laying. Thanks

[wide spear]

Sure LES is not always the right solution

AMR may not be the right solution either

Given that it can very quickly ruin your CFL condition

[sand heron]

Well, it looks complicated as it gets differently with the endlessly problems that arise in every way as it solves the others in back.

I have been granted a special capacity computer out of some bucket and some special tools and I wanted to try whatever possible. I'll look more to whichever solution is to be right. Otherwise, if not helpful enough I'll seek the coping mechanism since I don't have the capacity myself to approach higher.

It was a sort of project to be presented as a gift for my boyfriend's birthday and I have been finishing a good chunk of it but quite struggling on the last verses of that sonnet. It's either I sacrifice some stuff and present it off or just think about a whole total new project with less problems in the remaining time expectantly.

Thanks for your aid 😊

[fathom rain]

Maybe use nek5000/neko?

[hoary light]

Simulating boundary layers as a birthday gift, that is truly excellent

[tall badger]

Hello! Is Lagrange interpolation for a single point even possible?

[wide spear]

What does that even mean

You interpolate functions on intervals from data points

You cannot interpolate a function given a single data point

[tall badger]

My professor said that for the Neville scheme, having a single data point yields a constant interpolating polynomial. 😳

I then thought, what about Lagrange interpolation?

[wide spear]

Oh I mean sure

Lagrange will give a constant function as well

Any reasonable interpolation scheme should give a constant function when done with only a single data point

[tall badger]

It's because we know that for the single data point (x,y) the lone Lagrange polynomial, which when evaluated at x yields y, is unique, right?

[unreal mural]

I mean...

there are like 3 different ways you could come to that conclusion

lagrange interpolation with k data points will give you k-1 degree or less polynomial that passes thru them

without even knowing that they are unique in general, you know they have to be unique in this case

there is only one degree 0 or less polynomial that passes through some given y value

[wide spear]

https://visualpde.com/

[jagged cipher]

Hi I'm considering implementing the algorithm in this paper in mathematica and I was wondering if it was legit. It's claiming to very quickly compute the Jordan decomposition of large matrices and I'm wondering how issues of numerical stability factor into their algorithm https://link.springer.com/article/10.1007/s44198-023-00108-6

[wide spear]

Indeed, this paper is kind of sloppy

[jagged cipher]

Is it hard to compute the decomposition because you need to be able to distinguish/identify two eigenvalues as being equivalent, which requires symbolic representations for everything? Are they getting around this by exploiting knowledge about the algebraic and geometric multiplicities to maintain the correct structure of the Jordan matrix while still using numerical representations?

This isn't really my forte but I reduced my problem to a matrix problem because I figured it would already be solved

[wide spear]

If you can avoid a jordan decomposition then you should do that

There is the famous example of the matrix [[1, 1], [eps, 1]] where if eps=0 then the matrix is in jordan form and if eps!=0 then the jordan form is [[1+sqrt(eps), 0], [0, 1-sqrt(eps)]]

[hoary light]

the way to resolve that is to recognize that it's not really the jordan form that you want to compute, but the Jordan vectors

(nonetheless the problem is still ill-posed)

[wide spear]

For your problem, could you find another matrix decomposition that could work?

Schur decomp, SVD, etc...?

[jagged cipher]

I need a closed-form expression for matrix powers. For the Schur decomposition A^x = SJ^xS^-1, but the closed-form expression for J^x is too expensive. I don't think SVD allows any inverse matrix cancellations like Jordan, Schur, and diagonalization allow

My matrices are over the integer field if that helps at all with numerical problems. So in the famous example eps=0 or |eps|>=1

[wide spear]

And it isn't necessarily full rank?

Any other properties? Symmetric, hermitian, etc...

[jagged cipher]

It's a square block matrix, [[0,1,0],[0,D,A],[0,0,D]] of size 1 + 2k, where A and D are of size k. D is completely arbitrary. If an entry A_ij in A is non-zero, then A_ij = D_ij

So I don't think we get nice properties

[wide spear]

The first column is all 0s if I am reading that correctly?

[jagged cipher]

That's true

[jagged cipher]

Should step 10 say to add vectors from S that are not in the span of the P_ik?

[next garden]

Is someone here familiar with Hongkai Zhao's fast sweeping algorithm?
https://graphics.stanford.edu/courses/cs468-03-fall/Papers/zhao_fastsweep1.pdf

I am trying to understand what exaclty it does but I am having problems.

[wide spear]

Which part are you having problems with

[next garden]

Right now the update statement is flying over my head

https://en.wikipedia.org/wiki/Eikonal_equation#Numerical_approximation

Even in the 2D case

why do we take the maximum like that?

[wide spear]

So that it's positive