numerical-analysis

[wide spear]

The max is of two approximations of the derivative and 0 so the max of these three things will be >= 0

[next garden]

Sure but that's not really telling me why :p

It tells me that i guarantees something you want but not why it works if that makes any sense.

The argument seems a bit circular if that makes any sense.

Like why is it fundamentally true that the discretization must obey that property?

[hoary light]

Before we consider the sweeping algorithm (or any other algorithm), are you already comfortable with the mathematical properties of the underlying PDE?

This is essentially like upwinding, with the eikonal equation viewed as a steady-state of a Hamilton-Jacobi evolution

[next garden]

Which properties? Like I know that the general formulation is $|\nabla \Phi| = f(x)$ for some $f$ and $\partial \Phi = 0$.

And that this describes the modion of a propagating wave according to some set of rules.

[pine jetty]

Makogan

[next garden]

I know that the gradient is orthogonal to the surface at isocontour 0 i.e. teh boundary

[hoary light]

when f = 1 identically the problem is perhaps easiest to understand, because then you are computing the (signed) distance from the initial boundary

you want to trace along the characteristics and the distance continues to go up (or down) as you go further away

so if you think of this problem as time-dependent, where if you run the simulation forever you get the true solution to the eikonal equation, then this is some sort of hyperbolic PDE

and therefore you want to use upwinding for numerical stability

this explains the discretization

(assuming you want to use finite differences at all)

That is not a complete answer, but is this at least somewhat clear, or is more background information needed?

[next garden]

Are you computing the signed distance? I suspect both the signed and unsigned distance are both solutions to the equation aren;t they? THey both satisfy the gradient and boundary conditions

Let me try rephrasing my question.

Let us say we are given the solution already, i.e. we have phi.

The first order approximation of the gradient is just the vector of the central finite differences across each direction, i.e. teh discretized partial derivatives of order 1.

That's how you would get the gradient if the function was known

On Zhao's paper and wikipedia. if I understand proper ly, what oyu do is, you have a point at which you do not know phi.

You want to use the local information you have to compute phi

Thus you look at the 4 (2D) or 8 (3D) neighbours adjacent to your point across the axial dimensions. This makes sense to me

But then you are saying this:

That is a discretization of the gradient kinda

But it;s slightly different

forward differences would tell you that you maintain the grid ordering

i.e. you would do $u_ij - u_{i-1, j}$ to compute the partial derivative along the x dimension

[pine jetty]

Makogan

[next garden]

But you are not respecting that ordering

you are taking the smaller of the 2 values

which puzzles the heck out of me

first, why is it the correct thing to take the smaller of the 2 neighbours

and ssecond, why do we know this gives something meaningful?

when we initialized the neighbours to who know what?

[hoary light]

Have you learned about upwinding?

[next garden]

Not really, I just recently took a course on nuemrical analysis for PDEs but a lot of it flew over my head

and I don't think we ever covered upwind

[hoary light]

You can see this article: https://en.wikipedia.org/wiki/Upwind_scheme
Essentially, you want to choose a forward or backward difference based on the direction of the flow of information. This is an essential theoretical feature of hyperbolic PDE

[jagged cipher]

A is an n x n matrix of non-negative integers. If I ask mathematica to solve for the eigensystem of A the runtime increase exponentially with n. I suspect this is because mathematica is keeping a symbolic representation of the eigenvalues. If I want to solve for the algebraic and geometric multiplicities of each eigenvalue is this symbolic representation necessary?

[fathom rain]

you can compute with variable precision (first low then higher) to verify the multiplicity up to some tolerance

if you have a reasonable matrix it is not a problem

[amber slate]

uhh so i had a question

it's about finite difference

just confirming, if the value of V44 is already given by the lower boundary dirichlet condition, do i have to consider about the neumann condition for V44 (for the right boundary conditions) ?

or am i doing this wrong

danke:))

[wide spear]

Hopefully the two boundary conditions are compatible

[fathom rain]

Is it poisson second order fd? then make the 1D materices with [-1 2 -1] and set a 1 in one corner of each for the neumann. then A=kron(L1,I)+kron(I,L2) the values of the dirichlet are put in the rhs

[amber slate]

this was the equation

[pine jetty]

calicals

[fathom rain]

then I assume kron(L1,In2)+kron(In1,L2)+(2pi)^2*(kron(diag(xi^2),In2)+kron(In2,diag(yi^2))) for the matrix (and multiply with 1/h^2 for the first two parts)

[amber slate]

thank you

[fathom rain]

In1 is Identity of size n1 (probably n1 and n2 is same for you)

[amber slate]

indeed they're the same

[fathom rain]

(note that this is not how most people generate these matrices, but I find it compact and straight forward way to do it)

[amber slate]

i once asked a similar thing here and found out that your "method" worked just fine! :))

[tight geode]

Hope this is an ok place to ask this (there's a long debate over matrix multiplication order in #linear-algebra). I've run into a problem calculating the SVD of a matrix when the eigenvalues in S are repeated. Full question here: https://math.stackexchange.com/questions/4686323/svd-picking-u-and-v-when-singular-values-are-repeated
TL;DR: When the order of eigenvectors in U and V are indeterminate, how does one choose what they should be?

[hoary light]

Consider the minimal case: you have a 2x2 matrix that is just a * Identity where a is some scalar

From here it is clear that neither SVD or eigendecomposition are unique

and you are allowed to pick any orthogonal set of size 2 to be your basis

I suspect that this does not actually answer your question, so perhaps you can follow up on this to clarify more about your question?

[fathom rain]

if you talk about structured matrices, then for example the eigenvectors can be dst or perturbed dst, and then the eigenvalues can "probably" be ordered by them (this is a research question to order eigenvalues with close or same values, it is very difficult and I have not managed yet for large matrices, there are multiple tests to see that the ordering is good and it always fails somewhere for like n>30. I think a brute force search would be necessary which would be expensive). The ordering of singular values in descending order (for structure matrices) is just "wrong", and is chosen because it is a convenient ordering for most applications, so then it is "right"

[pine jetty]

Optimum

TheMipchunk

[tight geode]

Right, it needs to swap U's rows?

And then, just to confirm: In the case where I have 3 eigenvalues, and 2 are repeated, I can only permute the columns of S that contain the repeated columns - right?

[hoary light]

The need to permute your matrix to make it diagonal has nothing to do with singular values being repeated

[tight geode]

Sorry, in the case where it's nondiagonal and there are repeated eigenvalues

[hoary light]

You can permute any and all columns, regardless of whether they correspond to the same eigenvalue

You can permute whichever way you want to make the matrix diagonal

The reason you can do this to make your SVD is that every permutation matrix is an orthogonal matrix

[pine jetty]

Optimum

[hoary light]

Sure, in that case you don't need to change the third component. But what if that 1 was instead a 3. Then the whole thing should be permuted.

[tight geode]

Ok, thank you so much! I think that answers it

Now to write that out in code

[hoary light]

Great. I don't think I fully ever understood what it was you're trying to do, but I'll accept that you figured it out 🙂

[tight geode]

just make S valid in the confines of someone else's code ¯_(ツ)_/¯

Pseudoinverse relies on SVD, so if the SVD is invalid the pseudoinverse is invalid

[fathom rain]

"invalid"?

you can construct the pseudo inverse without SVD https://en.wikipedia.org/wiki/Moore–Penrose_inverse#Construction

[tight geode]

Yeah, it was just the easiest one for me to understand :/

[hoary light]

The QR method for pseudoinverses the standard method; the QR method is the same as doing Gram-Schmidt orthogonalization on the column space of the matrix.

However, given this clarification, I would be careful in the future about how to ask questions (which you should continue to do!) because your situation is a textbook example of an XY problem.

[tight geode]

In my defense, I was encouraged to pursue the SVD route because matlab and scipy/pytorch calculate the pseudoinverse using it. As for the QR method, I can't find a consistent methodology - do you happen to know a good resource on it? And is it preferred because it has better time complexity or stability? I can't find a good answer to that either ._.

[wide spear]

For how to compute the QR decomposition?

There are two main ways, householder reflections and givens rotations

Have you picked up a book on numerical linear algebra before

Trefethen and Bau is popular, Demmel is good as well

[tight geode]

Not how to compute it, but using the decomposition to get the pseudoinverse

Can't find a consensus on how that's done

[fathom rain]

if you want to see julia's implementation just do ```julia
using LinearAlgebra
@edit pinv(rand(Float64,5,5))

[wise swift]

heya! i have an exam tomorrow for my numerical analysis class. all we have been doing all semester is code but our exam, as our professor said, would be pen and paper. what would the questions be like on the exam because i really do not know what to expect? lol has anyone been on the same boat?

if it helps, the coverage would be on computer arithmetic and error analysis, nonlinear eqns and iterative methods, linear systems, and interpolation

[hoary light]

"find of the rate of convergence of this iterative algorithm"

"compute an upper bound for this interpolation"

Find a forward (or backward) error bound for this algorithm due to floating point arithmetic

Are those types of questions that were covered?

[bitter fiber]

(report from #360643390594875392) hey, there was a CTF with a slightly difficult challenge, it's now over and the numerical solution is quite convoluted, does anyone with expertise in algebra and numerical solvers mind explaining some things to me?

the problem basically boils down to solving
x * M = y
where len(x) = 22, some values in it are known, M is a 22x37 matrix and y is a known vector of len 37

using pinv only gives one solution, which isnt the flag

the solver published after the CTF ended utilizes the kernel of a matrix to iterate through the possible solutions until something that's representable as an ASCII string is found

I get the main idea but there's some core aspects that I don't understand

mainly what is achieved by multiplying two matrix kernels together

[warm otter]

could you provide some more context, maybe the statement of the problem or a link to the solution that you reference?

[bitter fiber]

i can post a txt file and an ipynb

model = nn.Sequential(
    nn.Linear(22, 69),
    nn.Linear(69, 420),
    nn.Linear(420, 800),
    nn.Linear(800, 85),
    nn.Linear(85, 13),
    nn.Linear(13, 37)
)``` 
this is the "model", notice no activation function so it's just matrix multiplications + adding the bias, the goal was to find input X that generates a given output

[midnight python]

the backward euler butcher tableau is given by this

which implies this ^

How does this show that $k_1=f(x_{n+1}, y_{n+1})$?

[pine jetty]

pramana

[cerulean minnow]

Hi folks, I'm Özge. If anyone among you is Turkish and knows statistics, can you write to me privately? Thank you.

[woeful jetty]

im completely lost on this

we never really talked about O notation so i don't really understand what i need to show

[graceful marsh]

Can finite difference methods be used to numerically solve the diffusion equation inside a disk?
Is the finite element method useful in such situations?

[wide spear]

Sure why not

[graceful marsh]

thanks.
sorry for my poor question.
i think , in rectangle domain, i define domain u[i][j] (for 0<=i<=m, 0<=j<=n).
but i cannot come up with a similar method in disk domain.

[wide spear]

Rewrite your problem in polar coordinates

[graceful marsh]

i haven't known analysis with polar coordinate. do u have some pdf or code (wishing c++) writing about it?

[grave spoke]

Or use a Shortley-Weller scheme

[amber slate]

soo i had a question

i tried solving the equation using finite difference

and i got this result

idk if this result "makes sense" because my prof told me that the analytical solution is

[pine jetty]

Who the hell is Cal?

[amber slate]

and well ... they're very different

can anyone help? thank you!

[wide spear]

Without further information all we can say is that your result is not correct

[amber slate]

and that further information is ... ?

[wide spear]

How did you discretize?

What matrix did you end up with?

How do you incorporate the boundary conditions?

[amber slate]

damn yeah a lot of blanks

wait

i discretized the equation using 51x51 nodes both in the x and y directions using the usual second order fd eqs

[amber slate]

honestly, this one i'm still struggling with

[wide spear]

Ok how do you incorporate the dirichlet conditions

Also hopefully someone else comes along and helps because I am now venturing into the egyptian desert

[hazy jolt]

Do you guys know how I can demonstrate this exercise?

It tells me that if A is invertible then prove that....

[wide spear]

What have you tried

[hazy jolt]

[lunar brook]

,rotate -90

[hazy jolt]

[hazy jolt]

look

[amber slate]

i still can't do this in python the solution compared to the analytical solution is still very different

any advices on how to incorporate the boundary conditions?

i still don't know what to do when the boundaries clash at corner points

[warm otter]

your boundary conditions dont seem to be enforced. in the image you sent, your solution should be 0 at (0,0) but there it is around 0.6. Moreover, going from (0,0) along the y- or x-axis your solution should increase, since it is roughly half a sine curve which is also not the case in your image

[wide spear]

The boundary conditions don’t clash though

[hazy jolt]

Hi guys

one question

For this exercise I did it like this

[wide spear]

What is a good upper bound for what

[hazy jolt]

What is a good upper elevation for

[pine jetty]

Awuita Fria

[hazy jolt]

yeh or no?

[amber slate]

i worded it wrong, i meant to say when ... let's say at the upper left grid point, this point is imposed by two boundary conditions, dirichlet on the left side and neumann on the upper side (i hope this is not confusing), what should i do?

[amber slate]

yeah i'm still trying to redo it

[warm otter]

take the one that is easier to implement :^D

especially since the boundary conditions are compatible

so in this case dirichlet would be easier i guess

[amber slate]

why did i have to stumble upon bvps again 😭 i just found out that my solution to an assignment two years ago were wrong

do you have any recommendations on what i should read to understand 2d bvps again with mixed boundary conditions (let's say dirichlet boundaries on the left and lower side, neumann on the upper and left side)?

[warm otter]

i dont have any recommendations but if you use FDM you can just use your calculus intuition

[amber slate]

oh also another question, for 2d grids and neumann boundaries, should i use the center approximation or the forward/backward approximation?

[warm otter]

i think the easiest to implement is to just take the finite difference between the boundary point and the point on the next layer inside the domain

for the other versions you need ghost points

[amber slate]

and that could be a hassle, i reckon?

[warm otter]

i find it a bit annoying to deal with haha

[amber slate]

thank you very much!! @warm otter i'll check my work again

[white oracle]

hii guys, has any tried solving any interpolation question by programming? (divided difference precisely)

[wide spear]

Very common thing to do

[magic spruce]

how do i go about sketching B-splines?

[fathom rain]

that looks like homework.

[hazy jolt]

Guys, could someone help me to answer a question about the finite difference method, I get to the part of the indices to get the matrix, then I was wondering if I take it from J=1 or J=2.

This is the equation he gave me and this is the graph with respect to the nodes the step size is 0.5

[fathom rain]

you can call your nodes whatever you want as long as you are consistent. probably most common is to call boundary nodes 0 and n+1 and a total of n internal nodes, in that case j=1,..,n for you

[plush flare]

Does anyone know how to approach this problem?

Reviewing for an exam and have no idea how to do this

[teal lava]

No

I don't have ideas

[wide spear]

What is the operation count for a LU factorization?

What is the operation count for solving Ax=b with pre-determined LU factors?

[wide spear]

When doing radial basis function interpolation, when is the matrix spd?

Well it's always symmetric

When will it be positive definite?

[wide spear]

Oh I see

[magic spruce]

isnt u already within [0,1]? im not sure why im being asked to convert

[wide spear]

Yes

[magic spruce]

then how am i supposed to convert u

this is what the solution says, but again. i dont understand why its being done in the first place

[wide spear]

https://arxiv.org/abs/2305.10575

[fossil badger]

suppose I've got a dataset (xi, f(xi)) and I wanted to interpolate it using a polynomial. so I want to use Hermite polynomials but I don't have the derivative given. can I use numerical differentiation to find the derivatives at xi of f and then use Hermite polynomials? cause it sounds like it could cause some problems to use hermite polynomials with imprecise derivatives on the data.

[wide spear]

Yes you can

However you do want to make sure that the numerical differentiation scheme is sufficiently high order

[wide spear]

<@&268886789983436800>

[vapid heart]

Also I think Catmull-Rom does that already

Which is a special case of cardinal spline

[maiden wagon]

Hi friends i have a small problem someone can help me plz i need to show that the string method is convergent whenever $f\in C^1([a,b])$ and $\dfrac{f(a)-f(b)}{b-a}f'(\alpha)<2$, my idea is try to proof this secuence is a cauchy secuence so $||x_k-x_{k+1}||=x_k-x_k+\dfrac{f(a)-f(b)}{b-a}f'(\alpha)||=||\dfrac{f(a)-f(b)}{b-a}f'(\alpha)||<2$. Is this idea good ?

[pine jetty]

Altaid

[wide spear]

What is the string method

[wide spear]

How do people feel about the various c++ sparse lin alg libraries?

Eigen, armadillo, trilinos?

[minor osprey]

I’ve used eigen before and I’ve liked it

[wide spear]

I have decided to not go in a direction that would involve sparse linear algebra in c++

[minor osprey]

Oh is this what you mentioned earlier

[wide spear]

Well the thing I mentioned earlier was a sparse spd matrix

But I have decided to not use it

[shy inlet]

how to find it since the function isn't y=Ax+b?

[grave spoke]

Gauss-Newton

[hoary light]

For explicit sparse linear algebra, I don't have too much experience. But if you really just need to supply an abstract routine for computing x->Ax and/or x->A^t x, I enjoyed Eigen due to its emphasis on abstract interfaces

[wide spear]

The desired routine would have been conjugate gradient but the point is moot

[serene veldt]

Can any numerical junky explain why the schrodinger equation behaves so chaotically when doing numerical integration 😭

[wide spear]

Your code has a problem

[serene veldt]

It almost certainly does. I honestly don't know too many schemes other than discrete Laplace (for 2D spacial derivatives) and leap frog for time evolution

I don't know if a symplectic scheme exists

[wide spear]

Ok yeah don't use leap frog

Why don't you use something like RK4

[serene veldt]

I thought rk4 only works with odes

Also I have to take into account boundary conditions

[wide spear]

You spatially discretize to handle boundary conditions and turn your pde into a system of odes

[serene veldt]

So rk4 should work just fine for time derivatives?

[wide spear]

Yes

[serene veldt]

Okay okay. I just didn't know if it being second order in space and first order in time made a difference

[wide spear]

Anyways don't use leapfrog schemes unless you are also using an asselin filter

[serene veldt]

i feel like rk4 is significantly more computationally expensive since I need to gather time and spacial derivatives for all the surrounding points too

[wide spear]

You can worry about computational expense after you have a solution

[serene veldt]

fine...I'll be back

[serene veldt]

ahhh, another thing I just realized that RK4 will most likely allow me to have a bigger time step without worrying about the solutions diverging rapidly

so it may not imply that much more computational expense in the long run

[serene veldt]

so Im just a bit confused over here

with k1 its a function of x and psi
k2 is now just a function of psi + k1*dt/2?

so do I plug in psi + k1*dt/2 over here

but not over here

[grave spoke]

What's the issue with leap frog?

[wide spear]

Well you know how the time derivative is done right

What happens is that you'll end up with solutions at t=0, 2*delta_t, 4*delta_t, ... and delta_t, 3*delta_t, ..... and they can diverge

[grave spoke]

You mean for the position and velocity, respectively?

[wide spear]

There's a version of it for systems with just one time derivative where you discretize u_t using a centered difference

[grave spoke]

Not sure if this is the issue you referred to, but I read this

[wide spear]

Yes

I explained it very poorly lol

[brave crypt]

Don't suppose anyone here is sitting on a pdf of "Chebyshev Series for Mathematical Functions" by Clenshaw (1962)? Can't seem to find this paper anywhere, and figured I would roll the dice and ask here--just in case.

[snow anvil]

I've been been looking into the universal approximation theorems for neural networks this week and this 2013 PhD thesis seems nice
http://www.matfis.uniroma3.it/Allegati/Dottorato/TESI/costarel/Sigmoidal_Approximation_Costarelli.pdf
I say that because the proof is more or less had to get by.
There's this unrelated 2017 blog post
https://mcneela.github.io/machine_learning/2017/03/21/Universal-Approximation-Theorem.html
and tries to recite the original 1989 paper. But that's a bit of a rough read.
Also fun to know that these results have a funny history of investigations. The 1989 proof uses Hahn-Banach to prove the networks not being dense in the space of continuous functions would be contradictory. 👀

[dawn viper]

if we're thinking of the same theorem the usual technique to prove such result is to let \int f dmu=0 for all f in your space (i.e. neural nets), and show mu=0 by various means
in this case it's a lot easier as you have many control over how can change f and stuff which gives extremely strong restrictions on mu

another approach would be to try to approximate all polynomials but that would be pretty hard to construct with certain networks

[deep comet]

I have a question, is it possible to make a numerical stable 3x3 SVD Algorithm without like QR Algorithm or approximation?

I found something to compute the eigenvalues of Symmetrical Matrix which gives me the eigenvalues for A^T * A
I compute the smallest and the largest eigenvalue and solve for the middle one via Crossproduct. The reason i do that is because if the matrix is of multiplicity 2 i would get the same eigenvectors since i just simple slove the 3 Equations which you have for the eigenvectors in a 3x3 Matrix

I then compute the left singular matrix from the right singular matrix (i'll send the formular later and further explanation later)

And lastly the Σ just consists of the square roots of the eigenvalues

Following things of my 3x3 Matrix are known:

• In the optimal case it is of rank 2

• In the optimal case it has Multiplicity 2

• Multiplicity 3 and rank 1 is never the case

What are your thoughts about this?
I implemented it in java if anyone wants to look over the code

[cursive nimbus]

What are you using for the visualization?

[wide spear]

Yes it should be possible, but why do you want to avoid QR? Other thoughts from reading this:

• What is of multiplicity 2? An eigenvalue?
• What does optimal case mean
• Does it give the correct results?

[serene veldt]

It's my own software I'm making with OpenGL. Using compute shaders for parallel computations and then using the rendering pipeline to display results

[cursive nimbus]

I see, thanks

[deep comet]

Ye the multiplicity is refering to the eigenvalues, i dont wanna use approximation with QR Algorithm since i think a direct solution is faster and more accurate

• What does optimal case mean:
  The matrix is built from real world data so the values arent perfect, in the optimal case it should be of rank 2 and have multiplicity 2, most of the time rank is pretty close to 2 and 2 eigenvalues are also really close to eachother, the problem is that this leads to numerical instability
  I couldn't pin it down why it behaves like that, but thats what i found out from my test

More accuratly i wanna decouple my Essential Matrix into rotation and translation, thats what i need the SVD for since, left and right singular matrix are simply the rotation

[wide spear]

Ok if you actually want to be using a SVD in practice, you should not be writing your own routines, especially not in java

You should instead be using something like lapack or another numerical linear algebra library written in a low level language like fortran or C

For example, https://netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga84fdf22a62b12ff364621e4713ce02f2.html

I don't know how much thought you've put into it or how much experience you have with numerical linear algebra but writing your own routines is generally a bad idea because of things like floating point errors

Also lapack routines will be highly optimized

[hoary light]

I was actually working on adding a 3x3 (and perhaps 4x4) extension to Lapack

[wide spear]

Oh fun

[hoary light]

the 4x4 case is not very well-studied, but the 3x3 case has some literature because it is used for when you want really accurate 3D ray tracing

if you search for 3x3 eigensolvers you can find some things

[wide spear]

And I assume the point is that you can do 3x3 faster than general n x n when n=3?

[hoary light]

yes, because you can do a non-iterative algorithm

[wide spear]

Yes right

[hoary light]

you can just directly compute the solution to the cubic (or quartic) in a numerically stable way

and of course once you have the eigenvalues, you can also get the eigenvectors by a variety of methods

assuming everything is SPD then the problem is OK

[wide spear]

Is it getting added

[hoary light]

no, I was just dabbling in writing my own code in the style of Lapack

I haven't submitted anything for distribution

incidentally MKL already has some of these small-n optimizations

but theirs is even better because their small n optimizations are hardware aware

[wide spear]

Oh yeah I mean MKL should have it

Given how many phds intel has lying around

[hoary light]

One of my former colleagues works on such things for MKL as well

my group also found a bug in Lapack around a year ago that had existed for many years, but I guess that routine had basically never been used so nobody noticed

[hoary light]

For the problem at hand, it is clear that there is a direct algorithm since eigenvalues can be computed by cubic solver. But if you just want something that is numerical stable and converges faster than QR iteration, use the Jacobi method.

And it will be a lot easier for you to implement.

[deep comet]

Ik ik, ill use c later for the entire project, im just doing it like that rn bc im more confident in my java skills
When I switch to c i'll build it on performance

Bc I just solve each equation, the 3 that you get when u write it as functions and not in matrices from, in the case of a 3x3 at least

I can show you later

[deep comet]

Nah i wanna solve it with the Cadarnos formulars, i got all that stuff implemented and working, I just have some numerical instability when computing the eigenvectors from my eigenvalues

[wide spear]

Well, how are you computing the eigenvectors?

[deep comet]

Just solving the 3 equations you have at a 3x3 matrix

A * x = lambda * x

Nothing to fancy, its the fastest way, but doesnt seem like the most stable

You have 3 variables in x and 3 equations, just solved for that

Another way would be to triangulate my matrix i know

[deep comet]

For everything else bigger than 3x3, i just use QR algorithm for now, without wilkinson shift tho, just basic QR Algorithm

The bottleneck is the realworld data anyways cuz of the noise

When the base program works tho i'll compute the propagation of uncertainty so i can precisely pin down which uncertainty plays the biggest role

[deep comet]

The thing just is, i hate using libraries + i learn more when I do everything by myself even tho it takes more time
Also in this case when we switch to c later on, i can write the operations specifically for my case which comes in handy when we have to compute a lot

[hoary light]

Is that not how you solve the linear system?

[deep comet]

Nah not through triangulation
I did it the old classy way where you solve each equation for a variable

I tested it also with triangulatipn but both has numeric instability in some cases
Especially when 2 eigenvalues really close to eachother or when they are exactly the same (never happens bc machine precision) cuz then you get back the same two eigenvectors

[wide spear]

I don't know what you mean by "solve each equation for a variable"

[deep comet]

Just rewrite the matrix vector equation as a system of linear equations?

And then solve each for a different dimension of the eigenvector x, y, z

In 3x3 matrix case

[hoary light]

but isn't solving a linear system precisely done via Gaussian elimination?

[deep comet]

There isn't just one way to do it, and gaussian elimination also has some numerical instability

Also its harder to implement and not that efficient than just having straight up equations hardcoded that give you the result

[wide spear]

#odes-and-pdes

[sweet mural]

Wrong channel to post in?

[wide spear]

Yes, off topic for this channel

[sweet mural]

Sorry though it might be fitting since I’m doing numerical analysis for deformation using Kirchhoffs plate theory

[wide spear]

Based on the slide it looked like the question was just about the pdes, but once you start doing numerics feel free to ask questions here

[sweet mural]

Okay, I first need to understand what any of they means

[idle quartz]

Why is the H1 space a bad space to proejct div(u)=0 onto?

This is closely related to why continuous finiiete element methods are a bad approach for the incompressible navier stokes equation

I want to understand the later primerally

[fathom rain]

Maybe check this thesis, defence in September https://urn.kb.se/resolve?urn=urn:nbn:se:uu:diva-501357

[wide spear]

I assume you have some basic familiarity with geometric mechanics formulations of pdes?

[idle quartz]

NO but do your best

[wide spear]

Ok so

There are lots of functions that satisfy div(u)=0 but aren't in H^1 right

Which is why we need to project

However, a lot of these functions that satisfy div(u)=0 (or are in H(div)) but no in H^1 are very fundamental to fluids

Consider the function u(x,y)=[-y/sqrt(x^2+y^2), x/sqrt(x^2+y^2)]

This satisfies div(u)=0 but it isn't in H^1 because of the singularity at 0

At the same time, this is basically the most basic vortex you can have

And you do want to capture it

Just some vague intuition

Even if you restrict H(div) further, the problem still persists

Say intersection of H(div) and H(curl)

This function u is in H(curl) as well on bounded domains

[obtuse furnace]

Where could one learn about this? Is there a good textbook or paper?

[wide spear]

Haha I don't actually know any geometric mechanics for numerical pdes

😵‍💫

[obtuse furnace]

Well thanks for gifting me with the phrase 'geometric mechanics formulations of pde' anyways. I think I've encountered the idea before but didn't have a name for it until now

[grave spoke]

The phrase "structure-preserving discretization of pdes" might also fit

[wide spear]

Oh yeah that's another name

It's certainly related but I don't know if they're the same

Probably

[idle quartz]

I don'tfl feel like that is enough. The flow doesn't need to have a singularity but CG is discouraged in FEM for incompressible flows.

[wide spear]

The flow I mentioned is not singular, but it's derivatives are not L^2 at 0

[wicked hawk]

Any help for algorithms? Merge sort and the likes

[wide spear]

CS server linked in #old-network

[ember nexus]

I want to solve the system of equations in the first picture with the system of equations in the second pictures after converting them accordingly through 4th Runge-Kutta method. I reduced the second-order by introducing another parameter "v". Is my final form in the second picture correct to apply the method?

[wide spear]

Almost but not quite

[ember nexus]

What is the missing part?

[wide spear]

Probably just a typo but check your x3 and v3

[ember nexus]

Oh, yes, a typo! Other than this, my system is correct, right? That is, I do not need to do further manipulation.

[wide spear]

Yes

[strange yoke]

So, say you have N equations, and N unknown variables. If the system has a certain structure to it where the solution at a point depends on the solution at the neighboring points, like it does in all PDEs, then the system of equations expressing the discretization will have a diagonal structure to it where only a small subset of the variables are related to the solution for a given point.

Consequently, if you choose the boundary conditions suitably and perform some clever algebra on the system, then you can compute an iterative solution of the system rather than some kind of matrix inversion process.

[wide spear]

I know how to solve linear systems and PDEs

Manju was not explaining how they were solving Ax=b

[idle quartz]

Let's say I have the weak form of the Laplace equation.
The problem is usually stated as find u in H1 that satisfies a(u,v) for all in H1.

When the we look for an approximate solution with then project on V subspace of H1.
V can be the space of first order polynomials. Why does the basis of V make a difference in the approximate solution? Fir example choosing lagrange basis or the {1,x} basis makes a difference. What theorem states that the basis are important?

[grave spoke]

Céa's lemma states that [\norm{u - u_n}{H^1(\Omega)} \le \frac{C}{\alpha} \inf{v_n \in V_n} \norm{u - v_n}{H^1(\Omega)},]
where (u_n) is the solution to the weak form in the finite-dimensional space (V_n). Using interpolation, you can estimate the fundamental error (\norm{u - v_n}{H^1(\Omega)}) by choosing an appropriate interpolant (v_n) in (V_n). How this estimation exactly looks like depends on your choice of approximation space (I don't think it depends on the basis for the same space though).

C is the continuity constant and alpha the ellipticity constant of the bilinear form btw

[pine jetty]

z...

[wide spear]

The way to think about this is that the more basis functions you have, the more types of behavior you can approximate

For example, say you want to approximate f(x)=x^3 on [-1,1]

You can't approximate it very well with only a constant, and if you expand your basis to {1,x} at the very least you can approximate the increasing behavior

But if you expand your basis to {1,x,x^2,x^3}, then you can capture the behavior exactly

Of course this is kind of a dumb situation but it's like how as your increase the order of a Taylor series then you can capture say sin better

This corresponds to increasing the things in the basis

[idle quartz]

Yes but
2 lagrange basis and the [1,x] basis span the same monomial space. However it seems like there is a difference in the solution when we do something like int grad u grad v = int fv

[wide spear]

No they should give the same solution

However, in practical terms, Legendre polynomials are "better conditioned" numerically and are easier to work with

[dawn mason]

Guys do you have an engineering question solved using bisection

I can't seem to find one

[idle quartz]

Bisection as in finding roots?

[dawn mason]

ye is ok

[idle quartz]

Still on this topic:
if Vh is the space of polynomials of degree n, what does Vh + Vh' represent? Is it a subset or a super set or is it the same as Vh?

[wide spear]

Not sure what your Vh' notation represents

[idle quartz]

derevtives of Vh

[wide spear]

Ok what space if Vh' then

[idle quartz]

Nevermind. I give up on this topic. Thanks alot

[wide spear]

Huh

[simple vigil]

you can find plenty in topology optimization

usually not a bisection in particular, but eg Newton Raphson and related are often in there.

the problem setup is essentially "I have a domain Omega and a set of physics (from which you usually want an easy notion of energy/ underlying energy functional) and maybe a few more constraints; I'll define a (material) density field over this domain, and I want to find the distribution of matter that satisfies the constraints and minimizes the effect of the physics"

[austere moon]

Hi ! Is there a good book for basic numerical analysis ?

[livid gate]

maybe https://www.amazon.com/Numerical-Analysis-Richard-L-Burden/dp/1305253663

my prof recced this to me

[dawn mason]

Thanks tho i already did the assignment just need to make a report

[vapid plume]

I recently read "An introduction to numerical methods and analysis" by James F. Epperson. A very good book with a solutions manual. There are some typos in the book and the solutions manual but they can be easily noticed.

[wide spear]

Oh was Epperson good

My intro class used Burden and it was not very good

[vapid plume]

Yes very good.

[wide spear]

Ok yes I just took a look at the table of contents and it looks nice

The order seems a bit odd to me though

[vapid plume]

I completed the 2nd edition with some sections (including PDEs) skipped. There are a couple of course outlines given for the student to follow.

Well it was my first numerical analysis book so I wouldn't know about the order.

[vapid plume]

May I ask why the order seems a bit odd to you?

[wide spear]

ODEs -> linear algebra -> PDEs seems odd to me

[vapid plume]

Ah okay. Yes I've looked at a couple of books, all do LA at the start.

[wide spear]

But also I think that the content for ODEs/PDEs/lin alg has other books that will be much more in depth and focused

[livid gate]

how do you decide if a book is good or not

[wide spear]

It's mostly the stuff in chapters in 1-5 that I have not found a good book for

[livid gate]

like how do I know I’m dying because the book is bad and not just because I have severe skill issue

[wide spear]

You sort of develop a taste for textbooks over time

Like you use multiple and vibe with some and do not vibe with others

[livid gate]

so in the end it’s subjective, huh

[wide spear]

Of course

[vapid plume]

Chapters 1-5 of Epperson?

[wide spear]

Yeah

[full night]

When transforming a 2nd order ODE to a system of 2 x 1st order ODE, the resulting equations seems coupled. Is it possible to solve the system numerically by considering each equation separately? I.e. do Forward euler on position and then do forward euler on velocities. For example when simulating newton's second law.

[wide spear]

No

[vapid plume]

http://faculty.washington.edu/rjl/fdmbook/chapter1.html ----> matlab/fdcoeffF.m and matlab/fdstencil.m.

In matlab/fdstencil.m, where does the following even come from in the last part?

err0 = c * (j(:).^n) / factorial(n);
err1 = c * (j(:).^(n+1)) / factorial(n+1);
if (abs(err0)) < 1e-14, err0 = 0; end % for centered approximations, expect
if (abs(err1)) < 1e-14, err1 = 0; end % one of these to be exactly 0.

I'm having trouble understanding where the formulas for err0 and err1 come from and why one of them will be 0 if we have a centered approximation.

[hard escarp]

Hi. Not really about a book, but video lectures. Are there sets of video lectures (on an advanced undergrad, or graduate level) on numerical analysis for PDEs? I mean talking about FDM or FEM. I mean from a more theoretical numerical analysis point of view. I've looked for those, but wasn't really successful in finding anything.

[brave crypt]

This is the book used for the Numerical Analysis course at Berkeley. I'm enrolled in it right now and using it and it's very good.

Also, I was wondering about the following question:

Suppose that $f \in C[a, b]$, that $x_1, x_2 \in [a, b]$, and that $c_1$ and $c_2$ are positive constants. Show that a number $\zeta$ exists between $x_1$ and $x_2$ with
$$f(\zeta) = \frac{c_1 f(x_1) + c_2 f(x_2)}{c_1 + c_2}$$

[pine jetty]

kenfps

[odd tusk]

suppose wlog f(x_1) < f(x_2). then the RHS is greater than f(x_1) (replace f(x_2) by f(x_1)). similarly replace f(x_1) by f(x_2) to get an upper bound f(x_2). so use the intermediate value theorem, as f continuous

@brave crypt

[brave crypt]

BOBLESA

FROM IB DISCORD

NO WAY

@odd tusk HELLO

[odd tusk]

WAIT U KNOW ME?

lulz

[brave crypt]

LMAO

[odd tusk]

that was agesss ago

[brave crypt]

DO YOU KNOW GODXEIN?

[odd tusk]

dm me

YES

[gritty rampart]

I was trying an idea I had but it doesn't seem to work. In class we made a program that approximates the brachistochrone numerically using newton's method to find a minimum

the minimum of the time function, that calculates the time it would take a falling object to complete the current curve

So I thought of doing something similar. The idea is to find an approximation for a given function f using other functions

Normally I would use minimal squares, but what if I want to optimize over a set that is not a subspace of functions? Or what if the inner products are too complicated to calculate?

So let $F: C[a,b] -> \mathbb{R}$ such that $F(g) = || f-g ||$

[pine jetty]

casiel368

[gritty rampart]

That norm is the induced norm from the usual integral inner product

Approximate the gradient and hessian matrix of F and apply newton

[wide spear]

Right so usually what you do is discretize your function space somehow by considering some finite basis

[gritty rampart]

So, say I have some function g:[a,b] -> R such that g(x) = ax + b

And I'll evaluate F(g), but now my inputs are just a and b

Let $L:\mathbb{R}^2 \rightarrow \mathbb{R}$ such that $L(a,b) = F(ax + b)$

[pine jetty]

casiel368

[gritty rampart]

Let's consider the next iteration for finding a minimum of the distance:\
$(a_{n+1}, b_{n+1}) = (a_n,b_n) - (H(a_n,b_n))^{-1} * \nabla F(a_n,b_n)$

[pine jetty]

casiel368

[gritty rampart]

The thing is, this doesn't work, and I can't figure out why

Not even in linear subspaces, so I don't think it's an issue on local minima

Gradient and hessian are both okay, I tested them by hand in a sufficiently non-trivial example

But the main problem is that hessian tends to become a tiny matrix, so its inverse just goes boom and my curve diverges

Or converges to a nonsensical one

This is the code:

def funcion(lista):
    l = lista
    f = lambda x: l[0] * (x**2) + l[1] * x + l[2]
    return f


def prod_int(f, g, lista):
    h = lambda x: f(x) * g(x)
    return np.trapz(h(lista), lista)


def norma(f, lista):
    return prod_int(f, f, lista) ** 0.5


def dist(f, g, lista):
    return norma(lambda x: f(x) - g(x), lista)


def grad(f, z, d):
    out = np.zeros(len(z))
    for i in range(len(z)):
        dz = np.zeros(len(z))
        dz[i] = d
        out[i] = (f(z + dz) - f(z)) / d
    return out


def hess(f, z, d):
    out = np.zeros((len(z), len(z)))
    for i in range(len(z)):
        for j in range(len(z)):
            di = np.zeros(len(z))
            di[i] = d
            dj = np.zeros(len(z))
            dj[j] = d
            out[i, j] = (f(z + di + dj) - f(z + di) - f(z + dj) + f(z)) / (d**2)
    return out

a = -2
b = 2
n = 3000
delta = 0.0001
pasos = 300
X = np.linspace(a, b, n + 1)

f = lambda x: x**2 + 3

F = lambda params: dist(funcion(params), f, X)

y = np.array([0.5, -0.3, 4])

for k in range(pasos):
    H = hess(F, y, delta)
    D = grad(F, y, delta)
    try:
        Hinv = np.linalg.inv(H)
    except:
        print(k)
        break
    y = y - Hinv @ D

And this completely breaks. It stops on step 2 in this curve

This is the hessian:

[[0.00109139 0.0003638  0.0003638 ]
 [0.0003638  0.         0.        ]
 [0.0003638  0.         0.        ]]

(in the last step where it breaks)

The cose proposes f(x) = x^2 + 3
Then it tries to approximate that with three parameters that represent some coefficients of the kind ax^2 + bx + c, and applies everything I said to the array [a,b,c]

Oh my, multiplying the hessian by 100 before inverting actually solves it

I'd still love to know why and find a better solution

Nvm, it's not stabe if the initial guess is not close to the real solution

[sage needle]

I'm having some trouble understanding how to treat the convective term in the approximation of non-stationary Navier-Stokes Eq. I'm using BDF2 to approximate the time derivative, and then a semi-implicit formulation for the advective term.
In particular $u^{n+1}\simeq b_0u^n+b_1u^{n-1}$. What troubles me is how this is then used in the simulations i'm reviewing. $\left[(b_0u^n+b_1u^{n-1})\cdot \nabla \right]u^=\left[(b_0u^n)\cdot \nabla \right]u^+\left[(b_1u^{n-1})\cdot \nabla \right]u^*$
Is it because at the time step n+1, the values at n and n-1 are given/known, therefore i can "use" the linearity?

[pine jetty]

thomashawl

[wild surge]

I want to test my eigenvalue solver against matrices with known eigenvalues. Is there anywhere that I can find a large database of matrices and their eigenvalues?

[fathom rain]

not that I am aware of. do you mean known analytically or simply precomputed (with double precision?)?

[wide spear]

There are some online

[simple vigil]

why not just build test cases against an existing implementation

eg, numpy if you've written in python

[wild surge]

Precomputed

[wild surge]

Convenience is all. But you're right, I could do that.

[brave crypt]

s

[fathom rain]

If you want to try your solver use hard edge cases like the Grcar

or simply a Toeplitz with 1s on first subdiagonal and on the second superdiagonal

both will fail for pretty small n using double precision

[simple vigil]

honestly, I think that would be much more convenient than hard-coding N test cases. literally whatever equivalent of

A = np.random.rand(N,N)
# This returns the eigenvalue/vectors pairs in lined-up arrays
np_eigs = np.linalg.eig(A)[1]
your_eigs = eig(A)
assert len(np_eigs) == len(your_eigs)
for exact, test in zip(np_eigs, your_eigs):
    assert exact - test == whatever_pytest_approx(0)

either way

I think sverek is the "real deal" resident numerical analyst

probably just listen to him 🤣

[brave crypt]

Which method is better newton interpolation or lagrange interpolation. , And do they both produce the same answer

[wide spear]

Depends on the problem

[brave crypt]

And do they produce the same answer

[wide spear]

Sometimes

[brave crypt]

Are there any maths techniques that can reduce error in approximations

[wide spear]

Yes

[wild surge]

Interesting, thanks! And thanks to you too @simple vigil

[fathom rain]

How does you eigenvalue solver work? Does it work for Hermitian and non-Hermitian?

[wild surge]

It's just naive power method for a course I'm teaching. I'm putting together a set of test cases that students can use for their homework.

[fathom rain]

maybe take Laplace FDM, then you have the analytical expressions too 🙂 \lambda_j=2-2cos(t_j) where t_j=j\pi/(n+1) j=1...n

[fathom rain]

maybe some good examples but no precomputed spectra https://math.nist.gov/MatrixMarket/

[wild surge]

Whoa, neat

[wide spear]

Oh yes this is what I was thinking of

[wide spear]

.pin

[wide spear]

.unpin

[copper abyss]

✅

[livid gate]

why was that pinned in the first place

[wide spear]

Furry nonsense I guess

Really morally lax

[brave crypt]

What is y(x)

And why does the k+1 the derivative of F(x) have one zero in the interval

[grave spoke]

any (n+1)-times continuously differentiable function interpolating the (n+2) data points

[grave spoke]

repeated use of Rolle's theorem

[serene veldt]

im struggling to come up with a simple FFT algorithm for interpolation in terms of Chebyshev Polynomials. I am using the Chebyshev points which are:
Xj = cos(j*pi/N) : 0 <= j <= N. Writing it in terms of the exponential yields: exp(j*pi*i/N)
Because it's not exp(2*j*pi*i/N) I am having such a brain fart with not being able to recreate a simple FFT algorithm

[odd tusk]

can you not just write it as exp(2ij pi/2N)

and use the algorithm for 2N?

[serene veldt]

thats what I was thinking earlier.

to be more specific Im trying to evaluate this series

with

if I choose an even N, would I just evaluate at k = 1 to K = N, then do k = 0 seperately?

if I do N = odd, for example like N = 3, then Id be evaluating at 4 points which should be able to factor into two seperate transforms

however that N = 3 is ugly in the cosines and makes it annoying for me to try and get a recursive form of the FFT for this case, hence my brain fart

[odd tusk]

if i write $u_i$ as $$\text{Re}\sum_{k=0}^{2N-1}a_ke^{\frac{2\pi i k}{2N}}$$ where $a_k$ are $0$ for $i>N$, then I think you can directly use the FFT algorithm

[pine jetty]

bobles

[odd tusk]

and if you think about it a bit more, you can probably simplify it a bit to avoid arithmetic with 0s

idk if thats the most helpful way tho

[serene veldt]

hmmm, that prob could work and is worth testing, I just have to be careful of memory constraints with doubling my array size (working on GPU)

oh wait

I could just do some conditional

[pine jetty]

bobles

[odd tusk]

then when you sum from 0 to N the second term, you can shift the index

that should work much better than the above

@serene veldt

[serene veldt]

hmmmm

should I still choose an even N?

[odd tusk]

i dont think N even matters right?

[serene veldt]

well if I want to keep cutting the process in 2

[odd tusk]

oh right

yes in the calculation you assume N is a power of 2

so it simplifies nicely

[serene veldt]

so if N is a power of 2, would I have to evaluate from K = 1 to N?

[odd tusk]

you can probably fiddle to get it to work for other numbers but it will get messy

[serene veldt]

the only issue I have with N being a power of 2 is that I'd be evaluating at 2^m + 1 points which is odd

[odd tusk]

yeah i think so and then fiddle with the remaining values separately?

[serene veldt]

okay just wanted to be sure. The algorithms are always messy when you want to do something slightly derivative from the "perfect" scenario lol

[odd tusk]

i think you only have to add/subtract k=0 and k=N terms tho

and those should be nice ones

[serene veldt]

yea not the worst

not going crazy :)

[limber walrus]

Can anyone help with 2.2

Tried everything known to man. Couldn't even find it online. All I know is it needs composite Simpson's rule with nodes x0=1... xm=m

Lecturer said it wasn't hard and I've spent days on it and it's really acting on my nerves

[wide spear]

What have you tried

[limber walrus]

Integrating everything either just left, just right, both. Integrating normally, integrating with Simpson

Also it should be -0.5 instead of -3.5

[night thunder]

Not sure if right channel:
Is there anything that can be said about the singular values of a real dxd matrix X after normalizing (to L2 norm 1) the columns (or rows)?

[fathom rain]

If $A$ is the original matrix and $\tilde{A}$ is the normalized one, look at the gershgorin circles of $AA'$ and $\tilde{A}\tilde{A}'$. probably some relation there

[pine jetty]

sverek

[night thunder]

Thank you kindly

[fathom rain]

since you scale the columns of A with the norm of that column

[night thunder]

I was actually reading on Gershgorin circles as we speak. Intuitively I'd assume that maybe the spread of singular values must decrease (meaning that the radii of the circles also decrease)

But this is my very rough logic rn

[fathom rain]

depends on the entries. if they are small it should increase right?

[night thunder]

Hmm yeah I guess, I'm not sure

[fathom rain]

what is the application of this?

[night thunder]

None really, it just crossed my mind this morning

I've tried the normalization empirically a bit and it looks to improve the conditioning number so I thought it should decrease the radii (in my logic above) but maybe there are corner cases

[fathom rain]

it can worsen the condition number too ```ju-17:53:59> cond(A)
48.374150078708276

ju-17:54:33> for jj=1:10
A[:,jj]=A[:,jj]/norm(A[:,jj])
end

ju-17:54:35> cond(A)
48.42003463895798```

it is an interesting idea 🙂

I work with Toeplitz and Toeplitz-like matrices and this just scales the matrix with a constant except in the corners => we get another matrix algebra and it will shift the whole spectrum

[wide spear]

This idea of normalizing to 1 shows up in a number of other places

For example people in ML normalize to 1 after each layer in deep neural networks if working with low precision so that the numbers stay in range

[night thunder]

Yeah I don't think it's as straightforward as I thought at first

[fathom rain]

you can lookup cauchy interlacing theorem, might be something?

[night thunder]

Looks like this theorem only cares about submatrices

[fathom rain]

probably only relevant for highly structured matrices, like Toeplitz where most submatrices are the same and this could tell you bounds on where the eigvalues can move

[night thunder]

But if there was some (hypothetical) operation that does decrease the Gershgorin circle's radii, does this imply that the spread of sing vals decreases?

Not necessarily I think, but it looks so close to being what I need

[fathom rain]

you need to also decrease the distance between the elements on the main diagonal

[odd tusk]

doesnt that just make it harder to accurately pinpoint where they are?

[night thunder]

Why? It makes the circles tighter

[odd tusk]

sorry i thought u said increase the radii lol

[wispy kelp]

i need help with some hw questions, would anyone be down to help me in dms

[wide spear]

#❓how-to-get-help

[quasi saffron]

can anyone explain on linear shooting method? what is it trying to find?

and how to do it, cause Im stuck at making the 2 IVP equations, u and v

[minor osprey]

Given a bvp,
[\ddot{x}(t)=f(t,x,\dot{x}),,,x(t_0)=x_0,,,x(t_1)=x_1]
You reduce into an initial value problem
[\ddot{x}_a(t)=f(t,x_a,\dot{x}_a),,,x_a(t_0)=x_0,,,\dot{x}_a(t_0)=a]
And you want to vary the parameter $a\in \bR$ until you find a solution to the ivp that matches the bvp.

Solutions correspond to roots of
[F(a)=x_a(t_1)-x_1]
The linear shooting method uses a linear $f$ in order to reduce into two initial value problems (these I think are described on Wikipedia, I’m typing this on my phone and don’t want to tex that much)
Essentially what’s happening is that you’re exploiting linearity of $f$ in order to create 2 ivps that average out to the boundary problem

[pine jetty]

Kirbemily

[minor osprey]

The first few lines are general shooting methods, but understanding the motivation (with root finding of F(a)) makes understanding the linear problem easier

Let me know if you need anything else (I know the second question is on construction, but I just checked and they’re described on the shooting method Wikipedia page)

https://en.m.wikipedia.org/wiki/Shooting_method

[quasi saffron]

thanks 💕

I'll take some time to see its construction i guess

[minor osprey]

Yeah there's a proof referenced too, but I think it might just be some rearrangement of terms

[wide spear]

https://arxiv.org/abs/2307.06787 may be interesting

[void smelt]

I'm having difficulty using the finite difference method to approximate the solution to the PDE $x\cdot u_x=t\cdot u_t$ with initial condition $u(x,t_0)=f(x)$ and boundary conditions $u(x_0,t)=g_1(t)$ and $u(x_n,t)=g_2(t)$. Work is attached- Would anyone be able to give me some advice to get me going in the right direction? Thank you!

[pine jetty]

daniel_the_seer

[wide spear]

Are you trying to do this numerically? Or something

Evidently something is not right because your solution doesn't depend on the boundary conditions

[void smelt]

Technically the presence of V_0(t) in my "solution" makes it reliant on one of the boundary conditions. But yea, not all of them are directly used. I wanted to find a way to make g_1(t) part of V_0(t), but I couldn't find a clean way to extract it from u(x_1,t) without breaking the assembly of V(t).

[fathom rain]

where do you have h_x and h_t in the end I see you call them h_x and h but then h_x disappears

[void smelt]

My bad, that's a typo. There's only supposed to be one h since I converted u_t to its continuous counterpart instead of its discrete counterpart.

[serene veldt]

I hope I am framing this question well enough...
Given that I have a set of values of a function on a line in the complex plane described by z = x + is, where -inf < x < inf and s = constant, alongside a representation in terms of a power series in 1/z, is this enough data to reconstruct my function on the real number line?

[wide spear]

Pretty sure this is #real-complex-analysis ?

[serene veldt]

its in the context of global spectral representations of functions in the complex plane for solving PDEs, thought it would fit here

the functions I am working with are not globally analytic on the real number line, however they are globally analytic (on a line) if I take an imaginary step s into the complex plane

the paper I am reading then does a transformation w = 1/z where w is described by a circle in the complex plane and does a fourier transform in terms of angles on that circle of the function we are working with

[odd cedar]

coming here from #linear-algebra message

re-deployment of the original question:

so, i am having abit of an issue. i want to solve this $x=\frac{1}{(A^{-1}b)\cdot k}A^{-1}b$ using some iterative method (currently trying successive over-relaxation), but am having issues where the matrix $A$ gets close to being singular, as then the iterative method fails. although i do conjecture that the overall equation should still be smooth in that case. and due to that wonder if there could be some other iterative method that solves the entire thing at once, instead of it exploding when trying to compute $A^{-1}b$ directly

[pine jetty]

eriksonn

[odd cedar]

all of k,b,x are n-dimensional vectors, and A is an n-dimensional square matrix

[wide spear]

Do you know anything about the A matrix

Also was your linear system originally in this form or did you manipulate it into this form

[odd cedar]

what do you mean by that?

one could say it was like this from the start

[wide spear]

Well normally you say that you solve Ax=b and not x=A^-1b right

[odd cedar]

yeah, but here i have that in 2 separate places

so swapping back to Ax=b form does not quite work

[wide spear]

A^-1b is still just a vector

[odd cedar]

yeah

one could say that i started with As=b and then get s from that, and then do x=s/(s dot k)

but s fails when A has near-zero determinant, even if x is fine

[wide spear]

Ok do you know anything about the A matrix

Is it symmetric? Positive definite? Banded?

Hermitian? etc...

[odd cedar]

hm, can try finding an example matrix

A can change around quite alot, but it does have some structure still

[wide spear]

If A were orthogonal that would be very nice

[odd cedar]

its definitly not

this might be close to an expected A matrix

[wide spear]

Ok so it's sparse

[odd cedar]

if A were larger (and it definitly can be) then it would be rather sparse

[wide spear]

Right

[odd cedar]

more sparse for larger ones

[wide spear]

Have you tried solving using something like gmres

[odd cedar]

i tried sor on it, and that worked mostly fine, except when it got close to having zero determinant

the thing is that the actual x i want does not seem to care about the determinant being zero

in the limit that is

[wide spear]

GMRES should have better convergence than SOR

[odd cedar]

would it still handle crossing the zero determinant smoothly?

if A changes continiously across that boundary

[wide spear]

Only one way to find out

[odd cedar]

this whole thing is more that i am trying to find another iteration method that gets x directly, instead of needing to compute s first

as x seems more stable than s

so i want to leverage that somehow

without needing anything super complex like gmres (as last time i tried that one it got very weird and it also did not work at all)

doing a 2d example in desmos gives this

[fathom rain]

have you tried doing it in high precision. feasable if A is not too large

[odd cedar]

i only have doubles

[scarlet spruce]

what is the green thing above?

between A0 and A1

is this the parallelogram made from the two column vectors?

[odd cedar]

yes

[scarlet spruce]

this seems to be unsymmetric

[fathom rain]

well your example is integers

[scarlet spruce]

try CGNR

it solves A^TAx = A^Tb

[odd cedar]

the example is simplified abit

[scarlet spruce]

so if even b is not in the range of A due to numerical error

it will still find you the solution that violates L2 least

if you start with an initial guess of zero you will get s^+ = A^+ b

where A^+ is the Moore-Penrose inverse

the code for it is in Saad's book

it's pretty simple

[fathom rain]

where does the entries come from and are you working on a fixed size or a sequence?

[scarlet spruce]

there's LSQR too if you need even better numerical stability

both require you to implement Ax and A^Tx routines

[odd cedar]

all these terms are too fancy for my poor brain

[scarlet spruce]

CGNR = conjugate gradients for the normal equations

Saad's book is freely available

[odd cedar]

in the other channel it was suggested to try solving $\begin{pmatrix} A & b \ k^T & 0 \end{pmatrix} \cdot \begin{pmatrix} x \ c \end{pmatrix} = \begin{pmatrix} 0 \ 1 \end{pmatrix}$ instead

[scarlet spruce]

https://www-users.cse.umn.edu/~saad/books.html

[fathom rain]

you can probably flip your matrix by a antidiagonal and use minres

[scarlet spruce]

go to the part with CGNR

[dense nexus]

$\begin{pmatrix} A & b \ k^T & 0 \end{pmatrix} \cdot \begin{pmatrix} x \ c \end{pmatrix} = \begin{pmatrix} 0 \ 1 \end{pmatrix}$

[pine jetty]

denascite

[odd cedar]

yeah, that one

[odd cedar]

it does not appear to need to be that way for the final outcome to work out

[scarlet spruce]

it doesn't matter, if b is in the range of A you'll get the solution also

just intiialise CGNR with a zero initial guess and you get yourself x = A^+ b even in singular problems

[odd cedar]

does "in range of" mean that b is contained within the span of A?

[scarlet spruce]

takes about 5min to implement

[scarlet spruce]

yes

[odd cedar]

feels like it would still be an issue of it trying to explode to infinity when being close to singular tho

[scarlet spruce]

no

A^+ b is well defined

the only thing you can get from a large condition number is slow convergence

but you will get convergence

not with SOR or SSOR though

[odd cedar]

it still feels abit odd to need to compute that intermediate variable when the final one is so well-behaved anyways

interms of stability

[scarlet spruce]

just implement it please and run it on As = b

[odd cedar]

hm

lets see if i can find that page then

[scarlet spruce]

x0 is initial guess A is your matrix, b is the right hand-side

you don't need to keep vectors around from the previous iterations

you need vectors b, x, p, r, w, z

convergence can be check for after the computation of beta

beta is |z^{k+1}|^2 / |z^k|^2 = |A^T(b-Ax^{k+1})|^2 / |A^T(b-Ax^{k})|^2

so you can store in the beginning |z_0|^2 = |A^T(b-Ax^0)|^2

then you can check how many times the residual has decreased since the beginning,e.g.

|z^{k+1}|^2/|z_0|^2 < tol^2, e.g. you can pick tol = 10^{-6} as an initial guess

there are better stopping criteria

but they require some more work

if 10^{-6} is not satisfactory try 10^{-8} or even smaller

[dense nexus]

I dont think this answers the question. they dont really want to solve Ax=b. they want to find an s with As in span(b) and k^Ts=1. and while you could first find the x and then scale it, this is not good if x itself has huge entries

[scarlet spruce]

they want 1/(k^TA^{-1}b) * A^{-1}b, no?

[odd cedar]

yeah

[scarlet spruce]

the thing in front is just a rescaling factor

if you have A^{-1}b you can easily compute it

[odd cedar]

ye, and that rescaling cancels out any scaling that A^{-1} or b picks up

[dense nexus]

but they want to avoid that

what if A^-1 b has huge entries

[odd cedar]

that is sortof what i mean, but it sounds like this alternate method may avoid that?

[scarlet spruce]

what are you proposing is to be done instead?

[dense nexus]

I think solving this works

[odd cedar]

the issue is that it tries to wrap around infinity and come back with negated values instead

as the determintant sweeps past zero

[scarlet spruce]

just try CGNR and see if you get what you want imo, again the implementation takes 5 minutes if you're on python or matlab, some even have it available as a routine

[odd cedar]

i am nowhere near any such neat math libraries, lol

[scarlet spruce]

how did you derive this

[odd cedar]

but i can try regardless

[dense nexus]

As in span b means As=-cb for some scalar c

so As+cb=0 and k^T s = 1

[odd cedar]

question is still if that is smooth across the det(A)=0 boundary

[scarlet spruce]

what are you programming on?

[odd cedar]

java

[wide spear]

Yikes yikes yikes

[odd cedar]

for minecraft actually

[dense nexus]

at this point you had quite a bit of time to just try it

[odd cedar]

do note that my current setup of things means that just doing the computation of Ax is rather fiddly by itself

[scarlet spruce]

just code it by hand

[odd cedar]

i am trying

gotta go for abit aswell tho

can come back later with results and such

[scarlet spruce]

this seems like an ill-posed thing

c is essentially unconstrained, it doesn't really help the problem imo

[dense nexus]

well c takes care of the scaling that you would otherwise have to do

after calculating A^-1 b

and then dividing

[odd cedar]

hm, this might get rather weird to set up

[scarlet spruce]

it doesn't as far as I can tell

you just introduced a free parameter

You have As+cb = 0 -> s = cA^{-1}b, c= 1 -> c = 1/k^TA^{-1}b , actually this might be fine

if A is non-singular

if A is singular then fat chance

[dense nexus]

if A is singular then its a whole different question what kind of solution you want

[scarlet spruce]

I really doubt that a rescaling will change things much tbh + it makes the system nastier, but they can try

they'll still need to implement CGNR in either case

[odd cedar]

i tried that and it did not work at all

[scarlet spruce]

what is "that"

[odd cedar]

this

[scarlet spruce]

yes, there's a caveat that you need to implement it correctly

[scarlet spruce]

well, does this match the screenshot?

I would say no

[odd cedar]

hm

it gets abit tricky as i dont have good ways of manipulating matricies and vectors directly, so i have to use loops and stuff all over the place

[scarlet spruce]

that's not the issue

you should just be much more careful when implementing things

I also don't know whether your matrix-vector multiplication routines are correct, but you first ought to fix the obvious mistakes

it will also be good practice for programming I guess

you should go line by line, write a comment with the line from Saad and below it implement said thing, and then double check that you implemented it correctly

currently your code is quite wrong, and there are very sloppy mistakes

[odd cedar]

lets see if i can comment some of this then

also that i dont quite get what the $\parallel x \parallel^2_2$ means in comparison to $\parallel x \parallel^2$

[pine jetty]

eriksonn

[scarlet spruce]

it's the 2-norm, i.e. \sum_{j=1}^n x[j]*x[j]

if you put a different index it can be a different norm

but here it is 2-norm, so the usual

[odd cedar]

so they are the same?

[dense nexus]

yes. the general case is $\norm{x}p = \left(\sum{j=1}^n \abs{x_j}^p\right)^{1/p}$ which is why there is the notation with the index, here just $p=2$

[pine jetty]

denascite

[odd cedar]

oh i see what i did wrong, lol

hm, this is atleast better, and does not explode anymore, but am not entirely confident that it is entirely valid

still flips around to negative infinity, but does not freak out along the way

[scarlet spruce]

It wasn't just one thing 😂

And it could also be that your mat-vec multiplies are wrong

Also 5 iterations would likely not be enough

Eitherway, this isn't really numerical analysis

[odd cedar]

comparing to the sor method the solutions seem fairly close when far from the singular point atleast

will continue messing with this later

and perhaps try other methods also

[odd cedar]

i checked it and the solution is correct

[scarlet spruce]

my suggestion is to implement this one correctly before doing that

[odd cedar]

seems correct as far as i can tell regarding its output

input b vector

b2 = A*x

for 10 iterations

larger matrix needed some more iterations, but again converged just fine

[scarlet spruce]

You can set up a stopping criterion based on the residual if you wish, so you don't have to manually set iterations

[odd cedar]

you said something about beta in that regard

think i fixed that

[scarlet spruce]

at the beginning you can compute z0sq = |z|^2, and then after computing beta you can have if (|z|^2/z0sq < tol * tol) return;

where tol is the decrezse in the relative residual that you wish

then the for loop turns into a while(true)

[odd cedar]

that is what i ended up doing, although i kept it as a long for loop just in case

[abstract wing]

i have a kind of dumb question - i've been given a problem to implement 6 digit rounding for floating point numbers where the mantissa has 6 digits, using base 10 (the input is a floating point number x, the output the mantissa m and the exponent e; we're allowed to assume x not equal to 0) - however, i think something is wrong with my code, because using it, i can't show that (for instance) 1 - cos(x) is numerically unstable near x=0, while sin^2(x)/(1+cos(x)) isn't. (i guess this is more a coding question, but i'm guessing my difficulty is arising from not fully understanding the math/setup)

after doing some reading, my tentative function is

`def round_me(x):
'''base 10, 6 digit rounding for floating point # x.
mantissa of x has 6 digits. returns mantissa m and exponent e.'''

digits = len(str(x))
digits_g1 = len(str(x).partition('.')[0]) # number of digits before decimal
e = -1*digits_g1
m = round(x*10**e,6)

return m, -1*e`

but as i said, i don't think this is right - any help/comments would be appreciated.

[wide spear]

This channel is not for debugging, it looks like python so try asking in the python server

[abstract wing]

ah, sorry

[wide spear]

No worries

[abstract wing]

i don't see the python server in the network channel - do you have a link?

[wide spear]

https://discord.gg/python

[abstract wing]

thank you!

[odd cedar]

question: when using CGNR, does performing some number of iterations through it, and restarting it with a single iteration that many times produce the same overall result?

so looping either green or red

so is doing red 10 times the same as doing green 10 times, interms of the overall output for x

the last line of p=z+beta*p seems to conflict with the p=z that happens at the start, so i am not all that confident

also what is the difference between CGNR and CGNE ?

they seem similar in structure apart from CGNE being abit simpler with less variables in it

i might try some of that stuff again

[odd cedar]

further reading has falsified this claim, as the p vectors are changed to be orthogonal after each step, and that is not maintained when restarting (if i understand things correctly)

[plucky kayak]

yes you wouldn't be doing the conjugate gradient part if you were to simply skip the entire conjugation part

[brave crypt]

how can i get the interval in terms of p?

[distant kite]

Anyone know how to solve an r and \phi differential equation ans polarplot in mable?

[wide spear]

#computing-software

[distant kite]

thank you

[serene veldt]

Is there such a strategy to numerically solving time dependent PDE's whose initial condition function has an infinity (singularity) on the real number line but is bounded and finite on a section of the complex plane?

So to solve it on that section of the complex plane then reconstruct real values in post

[wide spear]

Is the solution known to be analytic?

[serene veldt]

On the real number line it's never going to be analytic at 0. On the section of the complex plane the initial function describing the initial condition is analytic everywhere on the section

[wide spear]

Is the solution analytic for time after 0

[serene veldt]

Is there a way to determine if its analytic even if the solution is not known?

[wide spear]

Well

Doing analysis

[serene veldt]

Tfw numerical analysis includes doing analysis numerically

I'm a smooth brain sometimes

[wide spear]

#advanced-pdes and such

[serene veldt]

As far as I know the solution is said to be analytic outside the singularity for all time t.

[wide spear]

Ok well if it's analytic for all time, then at each time step you could compute the extension to the real line

[serene veldt]

What's the procedure for extending it to the real number line?

[wide spear]

Computing the taylor series then evaluating

[serene veldt]

Okay that's what I intuited just wanted to make sure

[wide spear]

I wonder if instead of evolving the solution in complex space you could instead evolve the coefficients of the taylor series

Sort of like a spectral method but taylor series instead of fourier series

[serene veldt]

Spectral methods is what I am doing currently 🙂

The issue is that the PDEs are highly non linear

So I sorta need their value representation

[wide spear]

Unfortunate, sorry for your loss

[serene veldt]

Just to make sure I understand correctly from a "theory" point of view. Is my solution being (complex) analytic the only way to ensure that my evaluation on a line section of C helps create a taylor series that would converge on the real solution?

[wide spear]

One thing you could also do is replace the singularity with something that approximates it

So like if you have a dirac delta you replace it with a very sharp gaussian

[wide spear]

Yes the idea is to use analytic continuation

Identity theorem

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

[serene veldt]

Ohhhh!!! This is exactly what I needed to hear lol. I was always worried that my evaluation in C would yield results that were no longer valid in my purely real solution

[plucky kayak]

how would you solve a set of linear equations with variables being matrices X,Z that has this kind of form? (it's a square matrix) $\left{\begin{matrix}
PXP + A^TZ = -C\
AX = -D
\end{matrix}\right.$

[pine jetty]

transparent_elemental

[plucky kayak]

it's just that as far as I know in order to solve a system of linear equations you must pass matrix of coefficients and right hand side to the library, but here it's not clear how you would pass the matrix because X is being multiplied from left and right by P

[serene veldt]

I think something fascinating worth point out is that my given line segment in the complex plane is given by the equation z = x + i*s where -inf < x < inf, and s is parameterizing how far I want to go into C. Under the transformation w = 1/z, my line transforms into a circle where I can evaluate a Fourier series on!

This only works because my solutions are guaranteed to approach 1 at Re(z) = +/- infinity, so its periodic on that circle

[wide spear]

Solve for X=-A^{-1}D first then solve for Z?

[plucky kayak]

well that's an underdetermined part, the entire set of equations forms a square system, but the bottom part has more unknowns than equations

[wide spear]

Wait what do you mean

What are the shapes?

[plucky kayak]

X is n by n, Z is m by n, A is m by n

P is also n by n

[random quartz]

So n>m? or m>n?

[plucky kayak]

m < n

[plucky kayak]

I think I may be able to solve this system by vectorization of X and Z columnwise and swap that coefficient matrix for tensor product of two matrices, even though I have no mathematical justification for this it seem to work on few test examples

[wide spear]

Is this even linear in the entries of x

[plucky kayak]

yes if you work it out you get something like $\ \sum_j \sum_i p_{kj}x_{ij}p_{li} = - c_{lk}$

(I ommited the Z term for the moment)

[pine jetty]

transparent_elemental

[wide spear]

Oh if it's linear in the components of x and z then vectorizing X and Z should work

[icy fossil]

Hello
Ive been struggling for now 4 days with the golub welsch algorithm especially how the weights of the gauss quadrature are related to the eigenvectors.
I've also now made a reddit post in askmath and hope someone could explain it to me : https://www.reddit.com/r/askmath/comments/155l1tu/need_help_understanding_the_proof_of_the_golub/?sort=new

this is all the work i did trying to prove it

if the equation on the lower left corner is true and what i did to get to the equation bottom middle is correct it would mean that the integrals of w(x) and w(x) *pn(x)/(x-xj) are related by pn-1(xj) which i dont see why

[candid timber]

Does anyone have references for determining a regularization parameter when applying Tikhonov regularization to solve Ax=b, but with exact data b? Most of the methods I see (L-Curve, generalized cross-validation, quasi-optimality stuff, etc) come from inverse problems and they assume noise in the data vector b which causes solutions to blow up when no regularization is applied. In my problem, I don’t have noise but need to solve Ax=b repeatedly where A and b are updated over time, and A is usually ill-conditioned. So I just want a better conditioned linear system, but don’t have the blowup issue due to noise.

[wide spear]

Could you latex this

[wide spear]

How is A updated over time? Are they rank one updates? Have you explored preconditioning techniques? Why do you want to use tikhonov regularization?

If they are rank one updates, then the sherman-morrison formula may be useful

[candid timber]

Well I’m actually solving ODEs for parameters \theta and the matrix/RHS both depend on the parameters. So the system is really:

A( \theta (t) ) \dot{ \theta } = b( \theta (t) )

Preconditioning would probably work, but I don’t have a good idea of how to pick a good preconditioner for this problem so I stopped going that route. Tikhonov regularization at least guarantees we get a better- conditioned problem

[wide spear]

Ah ok so we don't know anything about how the updates will happen

[icy fossil]

im not familiar enough with latex to present it like that i can try to split it up a bit

[candid timber]

Yeah rip

[icy fossil]

this shows that for any J vj = xj vj (xj being eigenvalues, vj are eigenvectors) that J(xj) is related with vj and that xj are the zeros of pn(x)

this is the christoffer darboux identity which is derived from the matrix notation of the recurrence relationship

and this are my notes for the calculation of the weights
from the left is just the definition of the weights while on the right is the christopher darboux identity used at xj

[icy fossil]

the red equation i managed to prove but still do not have the equation that is ultimativly used

[inland orchid]

Im not sure where this belongs but please let me know if its better off elsewhere.
I have been trying to proof to myself number systems and how they fundamentally work.
But I noticed an odd pattern, online it seems like theres no sort of decimal like numbers for any base other than our base 10
theres no base 16 8.A2?
Hoping to find info about this as it feels odd, in searching it only ever brings up base 10 which i assume is since deci tends to refer to 10
would "floating point" be a more valid term for this concept?

Ig our base's correct name is "decimal" but for somereason a "decimal number" to me has always been along the lines "2.3"

[dense echo]

Floating point in computers is almost invariably done in base 2, which is exactly like decimal fractions except to another base.
E.g. 11.101 is the binary fraction for the rational number 29/8.

[inland orchid]

weird, i assumed floating point would be...broken or very weird in base 2
like 1.1?

[dense echo]

That is 3/2.

There are very few situations outside textbooks where anyone bothers writing down numbers in this format for human consumption.

It's much easier to make the computer convert them to decimal, which tends to make more intuitive sense for human users. But there's no inherent mathematical reason why only decimal would work.

[inland orchid]

i want to imagine we can convert any base 10 floating point to another system now but i want to believe there will be an accruacy issue at some point
like converting 1.3 from base 10 to base 2

[dense echo]

(One exception: In advanced number theory, such fractions written in a prime base are part of how the "p-adic numbers" work).

[inland orchid]

ig in base 2's sense that could be almost an irrational number right?

[dense echo]

1.3 = 13/10 would be 1.0100110011001100... in base 2 or 1.4CCC... in base 16.

That's no more irrational than 15/7 being 2.142857142857... in base ten.

[minor crystal]

soooo i don't know where to put this, but i figured this is the best place

i've been mucking about with the remez exchange algorithm recently, and weighing it in different ways, and i've gotten the basic jist of it, except really for one part: and that's choosing the initial points to use when applying it to digital filter design

more specifically, it's choosing how many values to put in the passband vs. the stopband, based on the "weight" of the error in passband vs. stopband - i've tried finding some reading material on it, but it seems a niche enough topic that i kind of just can't find much

the reason that it's important is that the algorithm kind of falls apart when you choose incorrectly, you get more maxima than points...

[vernal narwhal]

For computers though you have finite bits so at some point you have to truncate the number (but it can still be very accurate as long as you supply enough bits).

But if your decimal is not some finite sum of (1/2^n) then you won't get a nice terminating binary number I believe. However note how each bit guarantees us your number to is with some interval that has a length of 1/2^n, n=#bits. So basically this tells you high your error can be from that exact value

,w (1/2^(15)

[pine jetty]

Results provided by WolframAlpha

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

Requested by brandon7716#0

[vernal narwhal]

,w (1/2^(14))

[pine jetty]

Results provided by WolframAlpha

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

Requested by brandon7716#0

[vernal narwhal]

Also realize when working with a decimal in binary you have the number left to the decimal which you can easily convert to a binary number so you're really only concerned with how to write a the decimal portion of it In base 2

Example: 2.5 = 10.10000000 (repeating 0's ) = 10.1 or 10.011111111111 (repeating 1's)

Notice the 10 Base 2 = 2 base 10

[vernal narwhal]

This actually shows up in base 10 as well. 3/10 = 0.3000000 (repeating 0's) or 0.299999999 (repeating 9's).

[manic saddle]

Can someone explain why PLU decomposition is numerically stable whereas LU decomposition isn’t?

[wide spear]

Pivoting

[manic saddle]

But what exactly about pivoting helps with stability?

[wide spear]

Do you know what pivoting is

[manic saddle]

From Wikipedia’s article on LU decomposition, it says a matrix A can be decomposed into one of these 3 categories

[wide spear]

Sure and?

[manic saddle]

So if LU can be “anything” as in there’s infinitely many combos that work

Then you could have ||L|| and ||U|| be really big

[wide spear]

Case 2 corresponds to a matrix that does not have full rank

[manic saddle]

Then your error can blow up

[wide spear]

In which case one of norm(L) or norm(U) will be 0

At any rate, PLU only helps in case 1

[manic saddle]

Oh would that mean you have infinite error?

Wait ok

Only certain matrices are “unstable” for the LU algorithm

[wide spear]

Do you know what the point of LU is

[manic saddle]

The matrix representation of Gaussian elimination?

[wide spear]

And do you have any texts that you are looking at aside from wikipedia

[manic saddle]

Not really just googling for pdfs from random unis and whatnot

I audited a class that covered this last semester but I can’t really remember the exact details anymore

[wide spear]

Ok there are two very well known books

Trefethen and Bau

Demmel

[manic saddle]

I’ll have a look

[wide spear]

Either will cover LU and pivoting in great detail

[manic saddle]

Fantastic, ty!

[little lantern]

I was reading this. Can someone help me with what is H_k in the article. It will be very helpful https://normalsplines.blogspot.com/2019/02/algorithms-for-updating-cholesky.html?m=1

[little lantern]

It is about finding cholesky factor in efficient way after adding/deleting a corresponding row-column.

[scarlet spruce]

Latex fails to render

[ornate plank]

Is this where I can get help for things done in MATLAB?

[odd tusk]

probably #computing-software ?

i think this is more for theoretical questions

[wide spear]

As bobles said

[little lantern]

Go to webpage, scroll down and you will find an option 'View web version', click on that then Latex will render.

[serene geyser]

does anyone know what runge's phenomenon is

I mean I know what it is

but why do most sources use (1/1+25x^2) as the common example for it?

[plucky kayak]

it's just a canonical example, but you can observe it in other cases as well

[serene geyser]

So, no specific reason?

I know some functions like the sine function have well defined derivatives of all orders and is a smooth function, hence it isn't affected by it, but why this function exactly?

[plucky kayak]

suppose that you want to interpolate a set of point (y_i, x_i) where each y_i was obtained via function y_i = a*x_i for some constant "a" and you've got like 10 points, if all your y_i match this formula exactly and you try to interpolate this with a 9th degree polynomial you'll get a function that looks like straight line (as expected)

[serene geyser]

Yeah

[plucky kayak]

but if you add a tiny error to one of the y_i then suddenly you get an interpolation that has same kind of oscillations as in the case of interpolating 1/(1+x^2)

[serene geyser]

oh hmm

https://www.desmos.com/calculator/exd6n33uxt

like this one right

or did you mean it in a different way?

[plucky kayak]

something like this

[serene geyser]

oh

also I was confused about another thing

runge's phenomenon occurs during polynomial interpolation due to equidistant nodes and high degree right?

these are also the factors that affect lebesgue's constant, so does runge's phenomenon affect the lebesgue's constant or is it the other way around?

[plucky kayak]

I don't know what's lebesgue's constant as far as I know oscillations occur even if your nodes aren't equidistant if you choose sufficiently large degree polynomial

[wide spear]

Runge’s phenomenon is what causes the Lebesgue constant to grow for poorly chosen points

[brave crypt]

asked in #book-recommendations but didn’t get an answer so I want to ask here: any refs for either a course or book on numerical approach to odes?

[wide spear]

Iserles

[brave crypt]

thanks

[sweet ore]

Hi I have a question

Suppose I have a function and I want to linearly interpolate it

How would I go about knowing what step size I should be using such that the error is under some bound?

Somewhat confused how to go about this problem

[grave spoke]

equidistantly?

[pine jetty]

Zanarcane

[grave spoke]

where s is the linear interpolant of f with step size h

[sweet ore]

i see

so uh

yeah how do i go about solving this derivative

or i guess is it possibly to manually derive?

[wide spear]

You just take the derivative

[pulsar cipher]

,tex Let $ f : \mathbb{R} \rightarrow \mathbb{R}, \hspace{10pt} \vspace{0.5cm}x \mapsto y = f(x). $ For data with errors $ x + \Delta x $ with small error $ \Delta x, $ it holds that $ \frac{f(x + \Delta x) - f(x)}{\Delta x} \approx f'(x).$

[pine jetty]

Sciencenjoyer

[pulsar cipher]

Why isn't this just equal to the derivative?

How does it matter that this delta x is an error if this is just the definition of the derivative of a function?

[livid gate]

because that’s not the definition of a derivative

[pulsar cipher]

ohh, because delta x is not going to lim->0

[pulsar cipher]

Example 2.9 We analyze Algorithm 2.4 for (naive) calculation of the smaller absolute root $x_2$ backwards. This means we consider the function \
$f(p, q) := p - \sqrt{p^2 - q}$ with $p > 0$. For $|\epsilon_i| \leq \epsilon_{mach}, i \in {1, 2, 3, 4}$, we consider the expression \\
$\leftp-{\sqrt{(p^{2}(1+\varepsilon_{1})-q)(1+\varepsilon_{2})}}(1+\varepsilon_{3})\right$

[pine jetty]

Sciencenjoyer

[pulsar cipher]

Could someone explain how to derive that expression

[pine jetty]

Sciencenjoyer

[pulsar cipher]

Having this information given

The first thing I don't understand is why there are multiple epsilon

[lavish stratus]

Wouldn't the multiple epsilon come from $(x+y)(1+\varepsilon)=x+x\varepsilon + y + y\varepsilon=x(1+\epsilon)+y(1+epsilon)$? I'm not sure about the rest though

[pine jetty]

Kenshin

[pulsar cipher]

I ment, why are there different epsilon_i and not just this one epsilon

[pulsar cipher]

Oh I guess for every operation in that expression, another epsilon should be used, ofc.

[pulsar cipher]

\begin{align*}
\textbf{Task 2.3} \quad & \text{(Stability analysis of addition)} \
& \text{Let the machine numbers } a_1, a_2, ..., a_n \text{ be given and their sum } \
& f (a_1, a_2, ..., a_n) = \sum_{i=1}^n a_i \text{ to be calculated. The computer has the machine accuracy } \epsilon_{mach}. \
& \text{Since rounding must be done after each intermediate result, the actual implementation is given as} \
& f (a_1, a_2, ..., a_n) = \left( \cdots \left( (a_1 \oplus a_2) \oplus a_3 \right) \oplus \cdots \right) \oplus a_n, \
& \text{where } \oplus \text{ denotes addition with rounding errors, defined as} \
& x \oplus y := (x + y) \cdot (1 + \mu_{x,y}) \
& \text{with } |\mu_{x,y}| \leq \epsilon_{mach}. \
& \text{Note that } \mu_{x_1,y_1} = \mu_{x_2,y_2} \text{ only if } x_1 = x_2 \text{ and }
y_1 = y_2, \
& \text{that means, the rounding error at } \oplus \text{ depends on the numbers } x \
& y\text{ that are to be added.} \
a) &\text{ Show that } f\text{ is backward stable, i.e.} \
& f (a_1, a_2, ..., a_n) = f\left(a_1 (1 + \epsilon_1),..., a_n (1 + \epsilon_n)\right) \
&\text{for suitable } |\epsilon_i| ≤ c\epsilon_{mach}\text{ and a constant } c\text{ independent of }a_i.
\end{align*}

[pine jetty]

Sciencenjoyer

[pulsar cipher]

How do I turn the my_x,y into epsilon_i ?

What I mean is: How do I turn the my_x,y form the given definition to what is to be shown, which only uses epsilon_i

Please ask me, if my question is unclear

[void smelt]

I've been learning about using the Finite Difference method to approximate first-order PDE's with independent variables $x$ and $t$, but I'm having trouble understanding the reasoning behind its primary weakness where a larger maximum $t$-value results in a less accurate approximation (no matter how small my $\Delta t$ is). Could someone explain this to me?

[pine jetty]

Daniel

[rich oriole]

Anyone have any interesting areas in numerical linear algebra for a 3000 word research project?

[idle phoenix]

how to prove this

In this article, How is 4.27 proven?

[brave crypt]

Confused on how the solution is 0.625:

[livid gate]

write out what you get for p_0, p_1, p_2, p_3

[serene veldt]

If I have a time dependent PDE where my initial condition function is analytic, how can I be sure the time evolution preserves analyticity?

[wide spear]

This is what the study of well posedness is concerned with

[serene veldt]

Ah okay! Will look into that

[serene veldt]

is it also dependent if we are looking for just real analyticity vs complex analyticity?