numerical-analysis

[wide spear]

I mean different functions are real analytic vs complex analytic

[serene veldt]

I meant if the condition of well posedness changes if you want to preserve real analyticity vs preserving complex analyticity

[wide spear]

Probably

[frosty stirrup]

Hey, can someone help me with this question? Is it possible to approximate SDE with ODE, by adding to the equation one variable, which represents noise? So that at each time step during numerical simulation, one value from a predefined vector of randomly created gaussian noise would be taken?
Chatgpt says yes, but I would really like to confirm that with an expert 😁 . Thank you very much!

[sage cosmos]

Hey guys, can somebody explain to me how to find the adjoint of a method? I'm trying to solve this exercises

I did not understand how to find it in a general case, not only this, so if somebody could help me, I'd appreciate it

[wide spear]

Sure why not

[frosty stirrup]

thank you!

[uneven thicket]

A little late but this is an actual numerical method called the Euler-Maruyama method

[frosty stirrup]

thanks again! this is also what I was reading about, but was not 100% certain.

[uneven thicket]

Great yeah it's the equivalent of the forward euler method for ito processes

You can also use something like the Milstein method which is related but honestly E-M is pretty good

[brave crypt]

Any lecture videos on numerical analysis?

[sterile vector]

@brave crypt no

No one has made any

It is a niche waiting to be filled

Make your own

[gilded tapir]

Hello. I'd like some help with my intuition that has failed so far regarding finite difference methods for solving a (P)DE:

Let's say we are using a finite difference method to approximate the laplacian and the method is exact for monomials up to order $p$. The truncation error $\tau$ of this then scales as $\tau \propto h^{p-1}$.

[pine jetty]

Faputa

[gilded tapir]

Now let's say we are solving a Poisson equation with such a method.
So we discretize $\nabla^2 u = f$ to get a linear system $Au = f$.
For the solution error $e$ it holds that $Ae = \tau$.

What I can't wrap my head around is the following. How do we expect the solution error $e$ to scale?
My naive reasoning suggests that since $\tau \propto h^{p-1}$ and a typical entry of the matrix $A$ scales as $\propto h^{-2}$ (it contains the weights discretizing the Laplacian operator), the solution error should scale as $e \propto h^{p+1}$. Does that sound correct?

[pine jetty]

Faputa

[wide spear]

No I don't think so

Error should still scale as O(h^p-1)

[dire summit]

nope

[gilded tapir]

How does that work? I guess because not every entry of the matrix A is of the scale $h^{-2}$? (For example, the rows corresponding to the boundary condition)

[pine jetty]

Faputa

[gilded tapir]

Ah I just checked by hand on a simple 1D case and saw that while the inverse matrix has many O(1/h^2) components there are still O(1) components in every row. I guess something like this happens this case as well and the scaling doesn't change. EDIT: Hmm but the O(1) components were only ones associated with boundary nodes..

[simple coral]

Anyone know of any reference which shows how to block-diagonalize one-sparse matrix?

[rocky jungle]

Okay I thing I found a crazy method to calculate pascal triangle diagonals mod N

So imagine this a list of some amount of 0s, let's say 4 of them: 0 0 0 0

Now add one to the first number: 1 0 0 0

Then for every other place in the list add the number is left of him, from left to right

1 1 (1+0) 1 (1+0) 1 (1+0) => 1 1 1 1

Then apply the same thing again and again, but remember that the addation is in mod N (Let's say 10)

so then we get 2 3 4 5, 3 6 0 5, 4 0 0 5

So it's basically going down the addition diagonal in bulk... which is probably something people already kenw.

[wide spear]

No

Wrong channel

#discrete-math

[wide spear]

No

[open sable]

Hello, it recently read a paper about some estimation method, which requires to convert an ODE model to a state-space model

However, in my case I have a Delayed Differential Equations model

Is it still possible to convert the DDE model to a SSM?

[chrome heath]

guys

how hard is complex analysis compared to real analysis?

because real analysis was brutal for me

[fathom rain]

wrong channel

[wide spear]

#real-complex-analysis or #math-discussion

[fathom rain]

Date: August 07, 2023
Subject: PhD Positions, Scuola Normale Superiore, Pisa, Italy

Four PhD fellowships are available at SNS for students interested in
pursuing a PhD degree in Computational Methods and Mathematical
Modeling for the Sciences and Finance. Areas of specialization include
Numerical Analysis and Scientific Computing, Quantitative Finance,
Stochastic Analysis, and Computational Physics/Chemistry/Biology.

The application deadline is August 24, 2023; start date is November 1,
2023.

For application information, please see
https://www.sns.it/en/phd-admission-competition/phd-calls

Further information about the PhD program is available at
https://www.sns.it/en/computational-methods-and-mathematical-models-sciences-and-finance-phd-programme

Informal inquiries are welcome.```

[dense glen]

was working on some asymptotics and can't get this technique (see this talk by Zagier at around 17 min ) to work - would appreciate some advice

briefly - given a sequence $a_n \sim c_0 + \frac{c_1}{n} + \frac{c_2}{n} + \cdots $, to extract the $c_0$ coefficent

you instead calculate the sequence $ \frac{1}{k!} \Delta^k (n^k a_n) $ for some integer $k$

[pine jetty]

rubixcyouber

[dense glen]

from math import factorial, exp, sqrt, e, comb as C

X, Y = [], []
A = []
for n in range(1, 20):
        A.append(sum([C(n,k)**2 * C(n+k,k)**2 for k in range(n+1)]))
        if len(A) >= 2:
                X.append(n)
                Y.append(A[-1]/A[-2])

def get_asymptotic_constant(x, y, power=8):
        # Given y(x) ~ c0 + c1/x + ...
        # This function calculates c0
        #
        delta = lambda k: lambda arr: [sum(C(k, j) * pow(-1, j+k) * arr[i+j] for j in range(k+1)) for i in range(len(arr)-k)]
        transform_seq = delta(power)([pow(x_, power) * y_ for x_, y_ in zip(x,y)])
        return [i / factorial(power) for i in transform_seq ] 

for i in range(3, 6):
        print(get_asymptotic_constant(X, Y, power=i))

here im testing the technique on the Apery numbers as per the talk but I am not getting the correct result lol

[wide spear]

Is this latex correct

[pine jetty]

rubixcyouber

[dense glen]

should be good now

all out of order now but

[wide spear]

Are the denominators supposed to be like that

[dense glen]

oops

[wide spear]

Or should the powers increase

[dense glen]

yeha

asymptotic series

briefly - given a sequence $a_n \sim c_0 + \frac{c_1}{n} + \frac{c_2}{n^2} + \cdots $, to extract the $c_0$ coefficent, you instead calculate the sequence $ \frac{1}{k!} \Delta^k (n^k a_n) $ for some integer $k$

[pine jetty]

rubixcyouber

[dense glen]

that should be correct

[wide spear]

Wait is this a coding question

[dense glen]

yeah because $ \frac{1}{k!} \Delta^k (n^k a_n) \sim c_0 + \cdots + \frac{c_{k+1}}{n^{k+1}}$

[wide spear]

Or a theory question

[dense glen]

it is both i guess

[wide spear]

Ok this is not the place for coding questions

But what is your theory question

[dense glen]

im pretty sure the theory is correct

or

i may be misunderstanding the theory

never mind

[wide spear]

Ok

I guess

[dense glen]

its that the algorithm gives the wrong result when im pretty sure the theory works

wait it works lmoa

[wide spear]

Nice

[dense glen]

tangential question: is there a fast but accurate way to evaluate $f(k) = \frac{k^k}{k!}$ for integers k

[pine jetty]

rubixcyouber

[plucky kayak]

for large integers you can approximate factorial by stirling's approximation, which will cancel out k^k, then it depends what exactly do you consider fast/accurate

[wide spear]

You can do fast and you can do accurate and there is a tradeoff between them

[fathom rain]

for how large k? In julia just @. f(k)=(BigInt(k)^k)//factorial(BigInt(k))

it is quick for k=1E6 at least

(and exact)

[dense glen]

how fast for 1E6?

I’m using python

[fathom rain]

most the time was priinting

let me see

[wide spear]

Python has arbitrary precision integer arithmetic by default

How are you computing this at the moment

[fathom rain]

benchmarked it was 1.5s on a macbook air m2

and then no approximation, you get it exact

[dense glen]

doing some weird stuff with product trees

calculating f(k) = prod([k/i for i in 1..k]) with a product tree

[dense glen]

actually i think the delay is coming from me calculating a sum with the terms:

$s_m = \sum_{k=1}^{m-1} \frac{k^{k-1}}{k!} \frac{(k+m-2)!}{(k+m)^{k+m-1}}$

[pine jetty]

rubixcyouber

[dense glen]

wanted to calculate the asymptotics of this sum (from a paper of Knuth if anyone is wondering)

calculating s_1 to s_k is already O(k^2)

[pulsar cipher]

Exercise 2.4 (Stability analysis of the scalar product)
The calculation of the scalar product $\sum_{i=1}^n x_i y_i $ with $ x_i, y_i \in \mathbb{R} $ and $ n \in \mathbb{N} $ shall be considered. The calculation is done successively with the help of the partial sums \
$S_{r},=,(x_{r}\boxdot,y_{r}),\boxplus,S_{r-1},\quad r\in{1,2,\dots,n}.$
\
If $\varepsilon_{\text{mach}} $ is the machine accuracy, then $ S_r $ can obviously be written as \ \
$S_{r}=\left(x_{r}y_{r}(1+\eta_{r})+S_{r-1}\right)(1+\xi_{r})$ \ with $|\eta_{r}|,|\xi_{r}|\leq\epsilon_{\mathrm{mach}}.$

[pine jetty]

Sciencenjoyer

[pulsar cipher]

And I don't understand the "obvious". Why is there no "error constant" for the addition? Why is it not:
\$S_{r}=\left(x_{r}y_{r}(1+\eta_{r})+S_{r-1}\left(1 + \zeta_r\right)\right)(1+\xi_{r})$
\ with $|\eta_{r}|,|\xi_{r}|,|\zeta_{r}|\leq\epsilon_{\mathrm{mach}}.$

[pine jetty]

Sciencenjoyer

[pulsar cipher]

The box-operations return the rounded value to the nearest machine number after performing the exact operation.

[warm otter]

the error constant for the addition is the (1 + xi_r)

[pulsar cipher]

HOW did I instantly understand this after reading your message but didn't get this idea thinking about it myself for so long??

my confidence is dropping so low doing math

Thank you nevertheless🥲

[warm otter]

lmao

[warm otter]

dont worry i think the rounding stuff is a topic that confuses most people at first

[pulsar cipher]

thanks for cheering me up

[vapid plume]

Can someone confirm if this is correct? It is supposed to be the second-order approximation for the fourth-order mixed derivative on the LHS.

[wide spear]

Can you present your derivation

[vapid plume]

Kinda did it rough but kept applying the central second-order approximation on each partial derivative

Sec I'll show

[vapid plume]

Anyone?

[maiden wagon]

Hi friends small question here i need to solve this, my idea is solve this equality y(xk)-(n{k+1}-n_k)/h =0 is this ok or i need other approach? other thing i think is i can use this theorem

[vernal narwhal]

Anyone have any good resources (introductory level) on the topic of numerical-analysis? Let me know (books, notes, lecture videos, ect)

[maiden wagon]

http://physics.ujep.cz/~jskvor/NME/DalsiSkripta/Quarteroni-SkaitMetod.pdf best book

[errant wagon]

The book looks amazing

[heady wedge]

Does numerical methods fall under num analysis?

[wide spear]

Yes

[wooden kiln]

I was given the following Runga-Kutta integration method

Though I'm a little confused about the notation. It seems that $u(t)$ is the actual function we are trying to obtain. But what does $f$ represent?

[pine jetty]

jacobjivanov

[wide spear]

Well

With an ODE

You write it in the form u_t=f(u,t) right

[wooden kiln]

Sure

[wide spear]

So that’s what f is

[rich oriole]

The notes my university produced are quite good I can send if you still want them

I pretty much did the course with them I didn’t really need the lectures

[vernal narwhal]

Sure if you would like to send them

[grave spoke]

I'm also interested in the lecture notes if you don't mind sending them

[rich oriole]

@vernal narwhal @grave spoke sent

[worldly coyote]

so I've been recently getting into the HP-15c calculator, and it can perform factorials, but only so large. I was able to find a program that approximates factorials and stays within the limits of the calculator by keeping everything as logarithms, but I have no idea why this approximation works as well as it does. doing a quick calculation on Desmos shows that, at least for the calculable amounts Desmos can work with, the limit of f(x)/Γ(x+1) approaches 1 incredibly quickly and stays there. f(x) definition to follow.

$f(x)=\frac{(51840x^3+4320x^2+180x-139)x^x \sqrt{2 \pi x}}{51840x^3e^x}$

[pine jetty]

GoldenPhoenix

[livid gate]

stirling's approximation

[worldly coyote]

huh. interesting

I'll look deeper into it.

[wide spear]

Stirling's approximation is fundamental

[worldly coyote]

you'll have to forgive me, this is my first dealing with higher level approximations of things. Does fundamental here mean that it can't be broken down very meaningfully in terms of what it's approximating?

[wide spear]

Fundamental as in very important

[worldly coyote]

ah

I'm noticing that only the last little bit here seems to be matching up with the definition on wikipedia. I'm most curious what that rational polynomial is doing there.

[livid gate]

simplify it

[worldly coyote]

1+1/12x+1/280x^2-139/51840x^3 was the original form in the calculator program.

(at least after deciphering the code as best I could)

[livid gate]

yes this is an extension of stirling's approximation using an asymptotic series

[worldly coyote]

huh, you learn something new every day

[wide spear]

https://en.wikipedia.org/wiki/Stirling's_approximation#Speed_of_convergence_and_error_estimates

The first few terms are listed here

[worldly coyote]

oh cool!

only takes about 3 seconds to run on my calculator, too, which is quite impressive for how much heavy lifting it's trying to do (the program itself calculates the logarithm of the approximation in parts, to keep the numbers small, and then displays it in scientific notation by calculating 10 raised to the fractional portion of the common log, and then simply keeping the integer portion for the characteristic).

[urban lodge]

I was reading paper about ETDRK4 method

and I saw matlab code about solving KS equation with ETDRK4

but I cannot understand the purpose of LR

https://doi.org/10.1006/jcph.2002.6995

[upper hull]

hey would anyone be able to help with a pretty simple analysis question, i just started analysis and would like some clarification

can anyone give a good explaination of digit chopping arithmetic

[wide spear]

What in particular are you confused by

[midnight current]

In deriving the error bound for $|f^\prime(x_1) - q^\prime(x_1)|$, what should the Taylor expansion be centered at?

[pine jetty]

e5l7@1

[wide spear]

What have you tried

[midnight current]

$f^\prime(x_1) = f^\prime(x_1) + f^{\prime \prime} (x_1) (x-x_1) + \frac{1}{2!}f^{\prime \prime \prime} (x_1) (x-x_1)^2$

[pine jetty]

e5l7@1

[midnight current]

however, i don't think that makes sense because the interpolating polynomial vanishes at $x_1$

[pine jetty]

e5l7@1

[wide spear]

This is not correct as written

[midnight current]

i feel like I should be using the variable spacing in determining where to center it at, but i'm unsure how

[wide spear]

Ok but does f'(x_1)-q'(x_1) vanish

[midnight current]

sorry, but how? was it just the missing $+$?

[pine jetty]

e5l7@1

[midnight current]

i'm unsure- looking at $f^\prime - q^\prime$, it's unclear to me how to manipulate them into common terms

[pine jetty]

e5l7@1

[midnight current]

i get $f^\prime(x_1) - q^\prime(x_1) = 2 f[x_0,x_1,x_2] (x - x_1) + \frac{1}{2!} f^{(3)}(x_1)(x-x_1)^2$, but don't know how to move on from here

[pine jetty]

e5l7@1

[dry ether]

we're supposed to show that newton's method converges linearly for this function, but I can't find a way to demonstrate that

nvm I think I got it

[wide spear]

Nice

[dry ether]

cool so, I was wrong

my classmates have no idea wtf is going on

why does newton's method converge linearly on even roots?

I've seen stuff on stackexchange, it's too terse and I can't follow all of hte assumptions. The stuff the grouptext for the class has is wrong

so this SPECIFIC case is $(x-1)^2e^x$

[pine jetty]

mr.mseeks

[dry ether]

they are arguing that, $\frac{(x-1)^2}{(((x^2-1)}$ has terms of the same order so it must be linear

that's... wrong

[pine jetty]

mr.mseeks

[dry ether]

because by that logic $e^x - c$ converges linearly, which is obviously untrue

[pine jetty]

mr.mseeks

[dry ether]

I am both extremely disheartened by getting disemboweled on my last homework, and unironically spending a couple days trying to get an answer even with full "google it" engaged and still putting "I don't know" on this hw

I know newton's converges linearly on even roots! I know it has to do with the concavity of the function, and newton's method failing to be well defined at our fixed point

but beyond that I'm kind of lost

[vapid plume]

@dry ether Download "An introduction to Numerical Methods and Analysis" by James F. Epperson. Go to page 133, read Theorem 3.7 and then section 3.11.4.

Hopefully that will help.

[mortal pier]

Is there a simple way to take a symbolic function as an argument in a function in matlab?

[wide spear]

#computing-software

[fathom rain]

>> syms x;
>> f=@(x)x^2;
>> g=@(x)x^3;
>> h=@(f,x)f(x)^2;
>> h(f,x)

ans =

x^4

>> h(g,x)

ans =

x^6

[vapid plume]

What would be the 2D version of this 1D problem?

-u''(x) - s*u(x) = f(x), 0 < x < 1, s >= 0
u'(0) = u'(1) = 0

[wide spear]

What do you think

[vapid plume]

I think -Uxx(x,y) - Uyy(x,y) - s*U(x,y) = f(x,y), 0 < x < 1, 0 < y < 1, s >= 0

Not sure about the Neumann condition though

[wide spear]

Yes this is fine

What have you considered for the Neumann condition

[vapid plume]

Hmm

Nothing yet

Maybe Ux(0,y) = Ux(1,y) = 0 and Uy(x,0) = Uy(x,1) = 0?

[wide spear]

Yes that seems reasonable

[vapid plume]

Thanks

[still forum]

Is it true that given an integrable weight w(x) and a function f in L2_{w(x)}, the expansion of f in orthogonal polynomials converges in L2 to f?

I know it's true for all Classical OPs (iirc), but couldn't find any results stating for more general weights

(ignore question it was answered satisfactorily in #advanced-analysis )

[rich oriole]

I’ve got a question on my assignment that says show that the optimal implementation of sin(x) is backwards stable

But sin(x) is not backwards stable?

Are these two different?

[tall solar]

Does anybody even have an idea how to estimate the interval of convergence of newton's method?

E.g. say f=(x^2)/(x^2+1) with simple root x=0.

I believe the interval of convergence of newton's method is (-1,1).

I saw an idea based on xn+1=xn-f'(xn)/f"(xn) implying some contraction idea about the errors en+1<= C en^2. From Taylor's theorem en+1=f'(En)/(2f"(xn)) en^2.

En is that weird epsilon from the remainder

I don't know what to do from here. Do I need to inspect this specific f and some up with a bound?like hard analysis

[bitter dirge]

Yeah it's very specific and requires knowing the root apriori.

A contraction is formed if the error ratio $|x_{n + 1} - s|/|x_{n} - s| < 1$ is achieved.

[pine jetty]

Cute Pink Unicorn

[dry ether]

is this a reasonable thing to say about divided differences? I'm worried the professor wants more but it's kind of all I've got

also I guess it's actually at the kth point but

[bitter dirge]

Secant method is root finding method. Newton divided difference is finite difference scheme for non equal distance grid…

It’s approximating the derivative at x0, the very first point. Divided difference is a forward difference scheme.

[wheat mist]

So i'm comparing my textbook and the solution to a homework problem, and in all the solution's i've seen, they solve for lambda while completely ignoring alpha. they say that this is linearly convergent because the limit evaluates to 1, yet the textbook says that lambda must be less than 1 to be linearly convergent. Is there something I am missing? am i looking at the wrong equation? i'm very confused

is alpha plugged in beforehand?

[sand heron]

can the squeeze theorem be used as a numerical method?

[wide spear]

What problem would you use it solve

[sand heron]

Wdym

[wide spear]

Well, you come up with numerical methods to solve problems right

[sand heron]

Well. This theorem doesn't find values of functions but just limits...

[dry ether]

did I do this correctly? it feels wrong that there isn't an (x-4) in there somewhere

[brave crypt]

can someone explain this?

[next garden]

Does anyone here know how to construct the matrix for polyharmonic funcitons with Dirichlet boundary conditions?

Like if you have a given tehcnique to discretize the laplacian and you are searching for:

$$\Delta u = 0$$

Subject to a set of bc?

[pine jetty]

Makogan

[next garden]

I know one can use, e.g. the laplace bbeltrami operator to dsicretice the interior portion as a matrix $L$ where the entries correspond to the cotan weights. But I really do not know how to deal with the BC.

[pine jetty]

Makogan

[wide spear]

Can you do this in 1D

[plucky cipher]

hello, my professor gave us this example for forward difference:

i don't get where the (e^2-1)^3 came from

[languid jackal]

You applied Delta three times

[plucky cipher]

shouldnt it be just (e-1)^3 since h is equal to 1?

[languid jackal]

c=e^2

[plucky cipher]

ohhh ((e^2*h)-1)^3?

[languid jackal]

Ye

[plucky cipher]

thank you

[keen gazelle]

any ideas im just tryna help a friend

prove that for all n of N* the limite=0

[coarse creek]

Try it for n=1

What does it mean for f to be continuous on R?

[wide spear]

#real-complex-analysis

[smoky flare]

cos(x) + sin(x) is bounded over R, lying in [-2, 2], and since f is continuous f(cos(x) + sin(x)) is also bounded. Since the top is bounded and the bottom goes to infinity, the whole thing goes to zero

[wide spear]

Bad

[gilded willow]

Can u guys solve this for me.

[wide spear]

1. I can barely read this
2. No

Please read #❓how-to-get-help

[dry ether]

am I stupid or does the (1,1,2,2) not make any sense

that's not our function values, and it's certainly not our derivative values

[wide spear]

So with Hermite interpolation you are also using derivative information right

[dry ether]

yeah but the derivative of ln(x) at 2 is 1/2 not 2

[wide spear]

So I assume that means that x_i=(1,1,2,2) means use the values of f(x) and f'(x) at x=1 and x=2

[dry ether]

yeah I guess it's my best shot

[wide spear]

Form the hermite interpolating cubic

Compare with the other interpolating cubic

And the exact value at x=1.5

[wide spear]

#❓how-to-get-help

[pulsar stream]

Any books that talk about catastophic cancellation?

[bitter dirge]

Catastrophic cancellation has a wiki I think.
It happens when taking difference between 2 numbers under IEEE float64 with relative magnitude ratio beyond 2^(54).

I have questions. It's kinda long...
A signal recovery problem assumed observing a signal $\hat u$ with noise, the task is to recover the original signal. The following theoretical framework has been proposed by assuming a signal $u:[0, 1]\mapsto \mathbb R$

$$
\min_{u}\int_0^1 (u(t) - \hat u(t))^2 + \alpha |u(t)|dt
$$

Using trapezoidal rule and finite difference on a the sampled, discretized signal $u_i, 0\le t \le N$ on equal time interval simplify integral into optimization objective function

$$
f(u) = (h/4)((u_0 - \hat u) + (u_N - \hat u_N))^2 + \sum_{i = 1}^{N - 1}(h/2)(u_i - \hat u_i)^2 + \alpha \sum_{i = 1}^N |u_i - u_{i - 1}|,
$$

Now, according to my knowledge, L1 regularization, under the context of LASSO, has Baysian interpretation of giving regression parameters Laplace Prior distribution.

My question in this case would be, is there similar Baysian interpretation for this type of reguarlization under the context of singal recovery as well? What are we mininizing... exactly?

[pine jetty]

Cute Pink Calculator

[bitter dirge]

Ping me if you have the right references that might lead to the answer of my questions.

[livid gate]

no

[brave crypt]

The minimization problem is just a continuous analogue of LASSO

[bitter dirge]

Shit, I should have written $|u'(t)|$ in there.

[brave crypt]

Oh okay

[bitter dirge]

Lasso is $|u(t)|$, not $|u'(t)|$

[brave crypt]

Then its not

[bitter dirge]

I will try rewriting it...

[brave crypt]

Can somebody help me verify my answer for the third question

[bitter dirge]

Fixed.

A signal recovery problem assumed observing a discrete signal $\hat u$ with noise, the task is to recover the original signal. The following theoretical framework has been proposed by assuming a signal $u:[0, 1]\mapsto \mathbb R$

$$
\min_{u}\int_0^1 (1/2)(u(t) - \hat u(t))^2 + \alpha |u'(t)|dt
$$

Using trapezoidal rule and foward finite difference on a the, signal, discretized into $\hat u_t, t = 0,\cdots, N$ on equal time interval simplify integral into optimization objective function

$$
\begin{aligned}
f(u) =& (h/4)((u_0 - \hat u_0)^2 + (u_N - \hat u_N)^2) +
\
& \quad \sum_{i = 1}^{N - 1}(h/2)(u_i - \hat u_i)^2 + \alpha \sum_{i = 1}^N |u_i - u_{i - 1}|,
\end{aligned}
$$

Now, according to my knowledge, L1 regularization, under the context of LASSO, has Baysian interpretation of giving regression parameters Laplace Prior distribution.

My question in this case would be, is there similar Bayesian interpretation for this type of reguarlization under the context of signal recovery as well? What are we mininizing... exactly?

Ping and Point me to some references and stuff would be great.

[pine jetty]

Cute Pink Calculator

[vernal narwhal]

For gaussian quadrature, can anyone point me to some resources for (1) why our choice of x_i's are always interior and (2) why they symmetrically placed. I have only worked through some basic cases where I followed the "rule" and solved for my weights and nodes for accuracy of degree 5

like how do I find why it is an "open" rule

[wide spear]

Gaussian quadrature points are given by the roots of the legendre polynomials which are always even/odd functions and thus have symmetric roots

[vernal narwhal]

But what about it being only interior points. I need to look more into this topic and my book never brought up this property

I just heard that this property Is important because we can evaluate Integral numerically where bounds might not be defined on the f(x) if we were using a different method

[wide spear]

All the roots of the legendre polynomial are between -1 and 1

[vernal narwhal]

okay 👍

[vernal narwhal]

can someone help me explain why these are equivalent approaches in computing the weights

Like this system of equations will be exactly the same as the Lagrange polynoimal if we have degree n polyniomals and n+1 points

[vernal narwhal]

by equivalent I mean why is doing this integration over the lagrange polynomials to compute w1, w2,...wn the same as solving the linear system

where our y vector is the what the true integral for 1, x^2, x^3, x^4,.... and our x1, x2, x3,... are the roots for the legrende polynomial

[upper hull]

this might be vague but can any give a good explaination of cubic spline interpolation because im trying to get this and it makes no sense, it just looks like the polynomials and matrices the book is making are generated randomly

[wide spear]

Do you have specific questions about it

[bitter dirge]

I might still remember.

The magic is comming from the vandermond matrix. The Lagrange interpolation, Hermite Interpolant(?), the Newton's divided diff, all has coefficients linked to the inverse of the Vandarmonde matrix. The Vandermonde matrix is what is displayed there in your lecture slides.

The actual story that answer you question it's a long one I believe.

[bright palm]

Can I ask kind of an embarassing question? I'm working on a method to simulate burgers equation in 1d. Ive gotten up to this step with understanding the code. I'm really curious about line 3. Unfortunately I've lost where the hell I came upon this code in the first place. Is this a valid method? Should I be using something better on that nonlinear term? This seems so much like cheating

[pine jetty]

jan Niku

[bright palm]

the results of the code don't look uncharacteristic, but this step has me very concerned

[celest leaf]

is there anyone i could talk to about approximating PDEs with initial and boundary conditions? Ive been writing a university thesis on a topic in that area, and struggling quite a bit with some parts and how to convey the ideas

[errant wagon]

Just ask your question.

[celest leaf]

well, its not one particular question i have, is the thing. well i guess the one question is "how do i write this"

I do have actual questions, but theyre more open-ended areas that i would find it helpful to understand better

like, Ive learned how to prove an energy bound on the solution under a semi-discretised system, where the spatial dimension is discretised but time is left continuous. thats the method of lines. how do i know this bound still holds in the fully discretised problem though, the model that my computer will actually run?

[wide spear]

This is indeed approximately true but also this approximation is low order so you can do better

[bright palm]

this is going to sound silly but what do you think is the next-simplest thing to implement here? Predictor Corrector?

[wide spear]

Why do you want something more complicated

[fathom rain]

if you want more complicated treat time as a space dimension, then do 2D discontinuous galerkin and then use pseudotime to go to steady state 🙂

[bright palm]

I'm not sure, I still have doubts about this method. I guess it works just fine for a low enough time step.

[wide spear]

Have you implemented it

[bright palm]

I have, yea

[wide spear]

Does it work

[bright palm]

it's a little ringy to be honest

but I can just lower time step to make it more reasonable

[wide spear]

Unsurprising

[bright palm]

yea

maybe I'm putting work in the wrong place then

I could address the ringing instead

[wide spear]

Given that it's burgers equation and blah blah blah

[bright palm]

yea

maybe stop here, I could choose to expand on it for a full project instead.

I appreciate the input, this really seemed like an invalid choice to me.

🙇‍♂️

[fathom rain]

example of what I mentioned (but with a non-conventional dg method) https://www.2pi.se/publications/pdf/p1.pdf page 645

[real palm]

does anyone own a version of magma calculator that i can run a few lines on?

the online calculator needs a bit more time

[fathom rain]

stop spamming everywhere

[real palm]

Im sorry, there is no specific channel for math codes and applications so I didnt know where to send it

[odd cedar]

are there any resources here on how to find closest-match solutions to ill-conditioned and possibly overdetermined matrix systems?

a possible version of my matrix looks like this, where all the capital letters are 3x3 block matrices

ones in red may or may not be singular

but the green ones are definitly nonsingular

[plucky kayak]

you can just call standard least squares function or choose an appropriate sparse solving method for better efficiency

[odd cedar]

would a least squares function be good on an ill-conditioned square matrix aswell?

[lavish stratus]

#computing-software

[odd cedar]

also that i dont know how to find what an "appropriate method" even is

apparently i have to QR factorize the matrix

[plucky kayak]

well yeah there's a few options to do that, it's best if the QR factorization offers least squares solving as a built-in option

at least QR if done right doesn't square the condition number, it will be as bad as original matrix was

[odd cedar]

but will least squares die horribly if the matrix has bad condition number then?

[plucky kayak]

depends on the implementation of the least squares function and how bad the matrix really is, at least the QR decomposition formally exists for any matrix including singular and rectangular

at least eigen in c++ provides documentation on available decompositions and least squares options
https://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html
https://eigen.tuxfamily.org/dox/group__LeastSquares.html
https://eigen.tuxfamily.org/dox/group__TopicSparseSystems.html

if you want to add more regularization to get better condition you can solve $||Ax-b||_2^2 + \lambda ||x||_2^2$ instead of $||Ax-b||_2^2$, which is equivalent to $||\begin{pmatrix}A \ \sqrt{\lambda}I \end{pmatrix}x - \begin{pmatrix}b \ 0 \end{pmatrix}||_2^2$, $\lambda > 0$ regularizes the system for larger values of lambda

[pine jetty]

Transparent Elemental

[odd cedar]

so i can tune that lambda to get this more stable at some potential loss of accuracy?

[plucky kayak]

well you will loose out on the accuracy of the obtained solution when applied to the problem, e. g. if you're fitting the polynomial it will be slower changing than before which may be better or worse depending on the problem

but numerically it's probably more accurate because smaller condition number

[fathom rain]

use high precision data type. your matrix is small

then the condition number is not a problem

why complicate things, code it in julia and use \

🙂

[odd cedar]

it is more that i am worried about a scenario where all the matrices in some row ending up as projection matrices onto the same plane, while the corresponding a constant vector approacing zero

as i have a feeling that such a setup would mess up somehow

[fathom rain]

no idea what your sub matrices are, but why not just try?

[odd cedar]

probably a good idea, yes

atleast gonna try the qr decomposition and get a feel for how that is ment to work

as i have never coded up that before

might come back at some later time with more questions, but this seems like a good plan for now atleast

[fathom rain]

where does these matrices come from?

[odd cedar]

they come from constraining rigidbodies to be aligned with each other

solving for what nudges to momentum and angular momentum of each body will result in the proper relative offsets in position and angle

so this, but with more constraints at once

and also having constraints on relative orientation aswell as relative position

the possible reasons for some of the sub-matrices to be singular is how i deal with having the constraints behave like a bearing joint, where it does not care about rotation around a given axis

and that ends up behaving like a projection

[winged snow]

Can we express the analytical solution of $$\int_0^{ k} x \sin\big(x \cos(x) + \sin(x)+x\big) dx$$ in terms of bessel functions or other ways?

[pine jetty]

Chalk

[grim bobcat]

#latex-testing

[dry ether]

sorry to double post a bunch of dense ass math, if you know what the second order approximation of f'''(x) should be you can just skip to the end and tell me if I screwed up my algebra

this disagrees with what my classmates got, but I can't spot my error

also there's a 7 in there and that worries me

[pine jetty]

mr.mseeks

[wide spear]

I’ve just woken up and can’t formally check it but yes that looks wrong

[long lake]

Please don't post off-topic stuff in here

Sigh I'm giving you one last warning

Please stop

[fathom rain]

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

the 6/7 looks weird

the rest looks reasionable

(have not looked at the derivation)

[fathom rain]

you are missing a h in the f' term in second last equation

[dry ether]

I ended up leaving my messed up scratchwork and just putting up the right answer, I'm pretty sure I set it up wrong

[eternal bronze]

i will use this channel in 10 months

and finally be in the advanced mathematics category 😄

[toxic grove]

Anyone has notes about Runge's Phenomenom?

[amber slate]

huh how would i approach this problem numerically

let's say i'd use RK4

[wide spear]

Where is your uncertainty

[pine jetty]

jan Niku

[bright palm]

However my professor did this in maybe the most confusing way possible and I'm totally lost on how CIF is applied here or what we are calling the function we are even applying it to. I have no idea how standard or widespread this calculation/method is but maybe someone can provide some guidance?

As much as I've followed is that we need to create some region in the complex plane that encloses all of the eigenvalues of L

[pine jetty]

jan Niku

[wide spear]

What is CIF lol

Oh cauchy integral formula

[bright palm]

yea

I think (?) he said its alright to assume that L is diagonal, or that we can always diagonalize the linear portion to end up with a diagonal matrix here

so then you have a bunch of decoupled applications of CIF

Even in this case the application is not clear to me

[wide spear]

If L is diagonal then what would you use CIF for

[bright palm]

I'm not certain

[wide spear]

Right let me think

[bright palm]

He'd said that its possible to compute $L^{-1}\qty(e^{\Delta t L} - I)$ using a taylor expansion but that this might create computability problems due to floating point issues

[pine jetty]

jan Niku

[wide spear]

Yes this is possible

[bright palm]

and that CIF has been shown to resolve these issues

[wide spear]

For CIF, I vaguely recall that it is possible to evaluate matrix exponentials with a contour integral of some kind?

[bright palm]

right, I believe thats what were doing here, in a really general sense I understand this part

[wide spear]

Anyways the correct way to evaluate the matrix exponential is to diagonalize L

[bright palm]

He'd said this leads to numerical issues or something

there is some kind of issue that motivates CIF

Sorry i know i sound confused, I think he has some kind of issue because when he writes he switches letters around and frequently misses stuff and I have the same issue so its like playing a game of telephone

what is the pathway to using CIF to evaluate a matrix exponential

where do you start

maybe I can work forward from there

[wide spear]

Diagonalizing L=UDU^* and then e^(delta t*L)=U e^(delta t*D)U^* is a very standard thing to do

[bright palm]

Sure, i completely understand this

How do you stick CIF into this?

[tall solar]

Is it possible to interpolate segments of curves? Rather than just getting the right value at a finite set of notes.

[amber slate]

the max thing

[wide spear]

Ok how would you do this for forward euler

[bold violet]

Anyone know a decent resource for understand butchers tableaus? They've suddenly appeared in my coursework but aren't mentioned at all in the textbook or any of the published PDFs with the course

[wide spear]

Have you looked at the wikipedia article on runge-kutta methods

[bold violet]

hmm i did yesterday and only confused myself more but somehow it makes some sense looking at it now

ye okey i get it

Excuse my penmanship but would you mind sanity checking my understanding here?

[wide spear]

Yes that's correct

[bold violet]

Alright thanks

[wooden kiln]

If I have the following PDE:

$$ \frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} + b \frac{\partial u}{\partial y} = 0 $$

[pine jetty]

jacobjivanov

[wooden kiln]

where a, b are constants

I can solve this system, iff I have periodic boundary conditions for u, by using the Fourier Transform to obtain the derivatives.

I can then timestep u. However, how can I maintain stability?

[wide spear]

Depends on how you’re time stepping

[wooden kiln]

I would like to use an RK method. I currently have an RK3, but it seems to "blow up"

As can be seen here

[wide spear]

What happens if you shrink the time step

[wooden kiln]

dt is currently a tenth of dx = dy. How low should I go?

[wide spear]

What happens if you halve the time step

[plucky kayak]

does anybody recognize the decomposition A = USU^T where U is orthogonal and S is tridiagonal? not sure if A being symmetric positive semidefinite plays a role here

[wide spear]

Real schur decomposition maybe

[plucky kayak]

I thought that had upper triangular S

[wide spear]

For schur decomp then S is upper triangular

[fathom rain]

https://www.phys.uconn.edu/~rozman/Courses/m3511_18s/downloads/householder1.pdf

maybe better https://en.wikipedia.org/wiki/Householder_transformation#Tridiagonalization

[plucky kayak]

I hoped there would be some name for that decomposition that I could maybe call one day

[prisma burrow]

I hope this is the right place to ask this, I was teaching someone finite element for generic differential equations and as an example I used y' = y = e^x to estimate the value of e, and it spat out this series:

$e\approx\left(\frac{2N + 1}{2N - 1}\right)^N$

[pine jetty]

Copper_irl

[prisma burrow]

this is extremely similar to Euler's definition of e, so I'm confident that it converges to e. However, I can't seem to figure out how to find the convergence rate and order. Experimentally it approaches e linearly at a rate of 1, but is there an easy way to prove this?

the actual finite element process is honestly more interesting than this but man the problem is a brainworm

[rich oriole]

Shrinking the time step generally wont help if the method is unstable

For PDE’s

[raven dew]

let $a_n = \Big(\frac{2n+1}{2n-1}\Big)^n$, then $\log a_n = n \log\Big(\frac{2n+1}{2n-1}\Big) = n \log \Big( 1 + \frac{2}{2n-1}\Big) = n \cdot \frac{2}{2n-1} \to 1$

[wide spear]

Oh I forgot about this

[raven dew]

@prisma burrow

[wide spear]

potato that's much more complicated than it needs to be

[raven dew]

How lol

[pine jetty]

potato

[raven dew]

tbf okay sure

[wide spear]

$\frac{2N+1}{2N-1}=1+\frac{2}{2N-1}=1+\frac{1}{N-\frac12}}$ so the convergence when raised to the N power to e should be clear

[pine jetty]

守沢千秋
Compile Error! Click the reaction for more information.
(You may edit your message to recompile.)

[raven dew]

That doesn't really constitue a proof right

I mean that is the intuition

But yeah one thing would be like

[wide spear]

sparkle numerical analysis sparkle

[raven dew]

$\Big( 1 + \frac{1}{N-1/2}\Big)^N = \Big( 1 + \frac{1}{N} + O(1/N^2)\Big)^N = \Big(1 + 1/N\Big)^N + O(1/N)$

[pine jetty]

potato

[raven dew]

I think

lol

[wide spear]

Yeah looks good

[raven dew]

No u

[prisma burrow]

Woah neat

[prisma burrow]

I was able to get to this initial state, but finding the convergence rate stumped me

Proving that it was e from there was real comfy because that’s just e_{N + 1/2}

Lim n-> inf, sequence becomes Cauchy

Proving it has linear convergence (same with e’s original def) is still tricky. I know of this formula but it seems sus:

\newcommand*\abswrap[1]{\left|#1\right|}
$\lim_{N\to\inf}\frac{\left|e_{N+1} - e\right|}{\left|e_{N} - e\right|^p} = r \leq 1$```

[raven dew]

Oh okay by convergence rate do you mean like the error term?

Maybe I'll try to estimate the error better than what I did there lol

[pine jetty]

Copper_irl

[prisma burrow]

I think it’d take me less time to just do the latex manually than fiddle with command syntax debugging rip

[raven dew]

Yeah so I think probabyl easiest is like

If you want to make it more rigorous that is lol

[prisma burrow]

$\lim_{N\to\infty}\frac{\left|e_{N+1} - e\right|}{\left|e_{N} - e\right|^p} = r \leq 1$```

There we go that bug was not chill in the preview lmao

[raven dew]

We can bound $(1 + 1/N)^N \le (1 + 1/(N-1/2))^N \le (1 + 1/N + A/N^2)^N = (1 + 1/N)^N + \sum_{k=1}^{N} {N \choose k} (A/N^2)^k (1 + 1/N)^{N-k}$.

[pine jetty]

Copper_irl

potato

[raven dew]

And then you can probably just take any good enough bound on N choose k to see that that sum is O(1/N), i think

[prisma burrow]

Oooh

[raven dew]

Or tbh easier just to consider the ratio lol

[prisma burrow]

So you don’t even touch that original def? That’s just like a “general form” you wanna land at?

[raven dew]

And show that it's bounded

Oh nah I was just saying like

I'm showing it is very close to the original thing

And then we know error on that thing

So we just add the errors

Nvm that's overcomplicated

[prisma burrow]

Hey if it works haha

[raven dew]

Note that $(1 + 1/N + A/N^2)/(1 + 1/N) = 1 + O(1/N^2) and so (1 + 1/N + A/N^2)^N = (1 + 1/N)^N (1 + O(1/N^2)^N$

[pine jetty]

potato

[raven dew]

But then that last bit is 1 + O(1/N)

So that works nicely

[prisma burrow]

I’ll have to poke around with this more since I’m not that familiar with using O notation outside of my 10 year old intro CS brain cells, I suspect this is where it came from heh

Analysis is too fun college has cursed me lmao

Anyway, thanks! I can explain where I got this weird formula from if you’d like, finite element is trippy as hell but is too “applied” for most math nerds

Although that’s prob something for #modeling

[raven dew]

Yee lol I would be interested

[wide spear]

Finite element stuff is good for this channel

[bold violet]

If anyone with understanding of boundary value problems & finite difference methods want to have a look at #help-36 that would be very much appreciated. Alternatively if you could reccommend any resources for learning the concepts

[wide spear]

You can repost the question here

[bold violet]

Running out of time for this weeks hw so skipped the last one, gonna ask a TA about it next week. Have another question though;

After reducing an RK method to $y_{n+1} = y_n \times R(h)$, $R(h) = 18h^2 - 6h + 1$

[pine jetty]

Joachim

[bold violet]

Am asked to determine for what values of h(step sizes) the RK method is stable

[wide spear]

Ok do you know how to do this

[bold violet]

It's an ERK, specifically Ralstons method, no guarantee that i calculated R(h) correctly but method should be the same. The wiki page for RK only talks about stability related to IRKs

No not a clue

[wide spear]

https://en.wikipedia.org/wiki/Stiff_equation#A-stability

Apply the method to the equation y'=ky for k<0

See for what values of h the solution goes to 0

[bold violet]

ODE given for the task is $y' = -6y$ $y(0) = 1$, so $y=e^{-6t}$, a bit confused which method you mean i should apply here.

[pine jetty]

Joachim

[bold violet]

$y_n = (R(hk))^n \cdot y_0$, solve for $R(hk) < 1$?

[pine jetty]

Joachim

[wide spear]

Ok so k=6

And set y_0=1

And see for which values of h this goes to 0

[bold violet]

Using the equation in my last msg?

[wide spear]

Yes

[bold violet]

Tried plotting it and playing around with values for n, limit seems to approach infinity as n does, guessing that either means i calculated R wrong or the method is never stable which I doubt is the case

[wide spear]

For a fixed h?

[bold violet]

plotted as a function of h

,w plot (18 (6h)^2 + 6(6h) + 1)^100, 0 < h < 0.2

[wide spear]

Why are you starting h at -0.2

Can you have a negative time step

[bold violet]

no that's fair just zoomed it to the relevant area

[pine jetty]

Results provided by WolframAlpha

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

Requested by jofredh

[bold violet]

Though this breaks down as i increase n any significant amount

[wide spear]

Anyways I don't know why you're plotting this

If you want (18*36h^2+36h+1)^n to converge to 0 then you need 648h^2+36h+1<1

Right

So you just solve for h

[bold violet]

That seems way simpler yes

That gives me the range $-\frac{1}{18} < h < 0$ though

[pine jetty]

Joachim

[wide spear]

Ok so what does that tell you about the stability of this scheme for this problem

[bold violet]

Unstable for any positive timestep?

[wide spear]

Yes

[bold violet]

Hmm I think I'll have to think about this some more and probably solve a few more problems for it to make much sense

But now at least I know how to get there

Thank you so much

[eternal drift]

where did the - with the 6 go to?

[eternal drift]

this one here

[wide spear]

Oh whoops

648h^2-36h+1<1

[bold violet]

Nice catch

Yes flipping the B coefficient mirrors the graph around the Y axis so we now get 0 < h < 1/18

[wide spear]

So now we have some small stability region

[bold violet]

Alright thank you again, this was very helpful

[chrome tendon]

Hi, can someone help with this problem? I treid to calculate abserr(f(a), f(a*)) and stuck with this equation. Can't undetstand how to modify further to use the given abserr(a, a^x)

[wide spear]

You know that abs(a-a*)<=0.001

[chrome tendon]

Yeah, but what we can say about abs(a^2 - a*^2)?

Is it just abs(a-a*) squared or what?

[prisma burrow]

"Galerkin's method of weighted residuals" was what i used for this, but my god its astonishingly hard to find any text about it that isn't already compressed and given to mechanical engineers

i wanted to test it using an ODE i know really well: y' = y

[wide spear]

$\abs{a^2-a*^2}=\abs{a-a^}\abs{a+a^}\leq0.001\abs{a+a^*}$

[pine jetty]

守沢千秋

[chrome tendon]

Ahh, yeah

Thanks

[wide spear]

And a=-0.7314256 and a* is within 0.001 of it so the smallest a^* could be is -0.7324256 so abs(a+a*)<=1.4638512

[prisma burrow]

the method is pretty hairy to solve, but setting it up isn't too bad. Computers can do most of the heavy lifting. You start with a model "element" of your system, usually a polynomial approximation. You then reformulate the constants of the model to be in terms of the outputs as an interpolating function. Galerkin provides an additional constraint on this model that causes the system to become solvable in a pretty damn stable way.

ODEs can be rewritten in the form of Operator(y) = 0, in my case y' - y = 0. We define the right side to be the Residual of the ODE, which is ideally 0 for the true solution, so y' - y = R. We then apply the model to the ODE to find the residual. Integrate R(model) over the domain of the element and set the result to 0 to find the constraint.

i believe this also means the error of the entire ODE has the same convergence as euler's definition of e

not a great way to approximate e^x but you can do better with higher order models

[chrome tendon]

If given function f(x) = x^3 - 3x^2 - 3, how I can estimate and show this?

[wide spear]

What do you know about the convergence of the bisection method

[chrome tendon]

Not much

Ahh, guess I undetstood

[wide spear]

Ummm

Ok

[chrome tendon]

You were referring to this equation?

[wide spear]

Yes that's part of this problem

[chrome tendon]

In this case n = 10 or 11?

I think 10 but not sure

[wide spear]

Why do you think 10

[chrome tendon]

Well, from the equation, as |r-Cn| = |r-X10|

[wide spear]

Ok sure

[chrome tendon]

And what should I do with |f(x10)| part?

f(r) = 0, so I got that |f(x10)| = |f(x10) - f(r)|

[wide spear]

Well have a bound on how far x_10 is from r

So you can also bound how far f(x_10) is from 0

[chrome tendon]

So its necessary to calculate f(x_10)?

Can it be estimated without calculating f(x_10)?

[wide spear]

It can be estimated without calculating f(x_10) I think

[chrome tendon]

Any hints how to do it?

[split belfry]

Is anyone able to help me over understanding bisection and newton's method error analysis?

Question:
A) Let f be a continuous function on the interval [a,b] and f(a)f(b) < 0. Suppose that a = 1, b = 1.1, and c(delta) = 1. Estimate the error in approximating the root after running two iterations of the bisection method. Do the same for two iterations of the Newton's method.

B) Let f be a continuous function on the interval [a,b] and f(a)f(b) < 0. Suppose that a = 0, b = 100, and c(delta) = 1. Estimate the error in approximating the root after running two iterations of the bisection method. Do the same for two iterations of the Newton's method.

C) Suppose that you can use ten iterations of either the bisection method or Newton's method. Explain which method would you choose and why for the values of a, b, and c(delta) given in the part a) and b)?

My understanding:
For bisection method, I knew that b - a/2^n < epsilon, so that's what I did for part A and B. And I have this information about Newton's Theorem and it's error is bounded by ∣r−x_n+1∣ ≤ c(δ)∣r−x_n|^2

I also have this example from my professor:
|r-x_0| <= 0.01 c(delta) = 1

|r-x_1| <= c(delta)|r-x_0|^2<= 1(0.01)^2=(10^-2)^2=10^-4=0.0001

this is better than bisection because you would need to do 10 iterations of bisection.
|r-x_2|<= 1(10^-4)^2=10^-8

I don't understand how to apply newton's method to this question if I do not have the function, and I do not know what it is bounded to get |r - x_0| like how it is in the example from my professor. Am I missing something? How do I approach this problem?

[wide spear]

What is c(delta)?

Anyways the point is that you don't need to know what the function is to bound the error in approximating the root

Oh ok I guess that's what c(delta) is

So you have that abs(r-x_0) <= 1

Then abs(r-x_1) <= 0.01*1^2 =0.01

Then abs(r-x_2) <= 0.01*(0.01^2)=1e-6

For (a) that is

However, for (b) we have that abs(r-x_0) <= 100

So abs(r-x_1) <= 0.01*100^2=100

And abs(r-x_2) <= 100 as well

So for (a) we have that newton's method converges much faster than bisection

But for (b) we don't know that newton's method will converge at all so bisection is better in this case

[split belfry]

sorry, c(delta) is just c(δ). I didn't have the symbols at first so I typed it out. Unless you mean in the problem then c(delta) = 1 and I think c(delta) is the boxed red part in the image here of Newton's theorem (? again sorry I am still trying to understand this topic)

[wide spear]

Yeah sure c(delta) is not standard notation but it doesn't matter too much

[split belfry]

where did you get 0.01 in 0.01*1^2? Did you base this off of the example? that's what I am having trouble with on understanding how I know I needed to use 0.01

[wide spear]

c(delta)=0.01

Oh wait c(delta)=1

Lol

Whoops

Ok so in the first one abs(r-x_0) <= 0.1

So abs(r-x_1) <= 0.1^2=0.01

And abs(r-x_2) <= 0.01^2 = 0.0001

[split belfry]

how do you know that abs(r-x_0) is <= to just 0.1? what happened to the part where |r - x_n|^2 in the formula for this step?

[wide spear]

Well there is no x_{-1} right

x_0 is the initial guess

And it is in the interval [a,b]

And the furthest it can be from r is b-a

[split belfry]

Ohhhh okay, that makes sense. I didn't know what to do with [a, b] when it was given for Newton's method. Thank you!!

[lavish palm]

I'm needing some help showing the proof for the flops count for Gauss-Jordan elimination. Everything im finding is for just Gauss elim

[wide spear]

It’s the same

[near geode]

The only time I've seen problems like this is in physics. Wondering what the variables are to plug in.

[chrome tendon]

Can't really remember but the thing was that calculating 11 approximation by hand is really hard and considering that this was an previous year midterm example, were calculators are prohibited, I was sure that calculating it was not the way to solve it 😅

[wide spear]

?

[near geode]

When I was a child I had a Physics textbook. I got some understanding of the basics and saw some problems. Mathematics is considered the queen of the sciences.

[wide spear]

OK？

[bright palm]

embarrassing but I need a hint as to how to start on this thing. I've spent a pretty great deal of time reviewing the pre-req knowledge for this problem but I must be missing something.

We would want our basis functions to be bessel functions, right?

Should it be clear to me what role this matrix is playing or how it's even connected to the problem? Why should a matrix pop out as a reasonable thing to construct here?

beyond just, it's a numerical method so constructing a matrix is probably going to happen at some point. I understand that.

I understand that 9.9 is an eigenvalue problem. What is the relationship between this continuous equation and the eigenvalues of some matrix and what would be the bridge between the two? I'm so lost for a way to even start here.

[fathom rain]

You use some spectral method (I assume you know which this refers to) to discretize the problem so you get a system like $A_N\mathbf{u}_N=-\omega^2\mathbf{u}_N$. Then $\lambda_j(A_N)=-\omega^2$ for $j=1,\ldots,6$ is computed using a numerical solver. Then, I assume you check that it is correct by looking at the roots of the Bessel functions.

[pine jetty]

sverek

[brazen zenith]

Hello, I have a question:
When giving a set of points and asked to approximate the area under the curve, can u use the simpson's method?
Also what does solving for the integral that is promotional to the lower and higher bounds mean?

[opaque crag]

I was wondering if anyone with an HP Prime G2 would be able to help me as I can't find the answer anywhere online and I do not know what to search. I also do not know any other channels/mediums to ask this question. I am trying to transition from a TI-NSpire. I have the calculation ((x^2)6e^2x-3e^2x * (2x))/(x^2 +1)^2 @ x=3 1.275.

On the NSpire, I would type out the equation ((x^2)6e^2x-3e^2x * (2x))/(x^2 +1)^2 then type | x=1.275 after it in the same entry as it, and it would solve at that X. Is there any way to replicate this on the HP Prime? The only way I can figure, is setting my x with an 1.275 sto-> x, before typing my equation, but if I forget to do this, it is a hassle to work backwards.

[wide spear]

#computing-software

[bright palm]

I'm not certain which method it's refering to. I have really minor exposure to converting to a matrix with just like a super straightforward wave equation but it's not clear to me what the process would even look like for something this complicated.

[wide spear]

How did you discretize for the wave equation?

[bright palm]

turn it into a system of first order equations

so a super small matrix, 2x2

then used backwards euler

is a similar thing supposed to happen here?

choice of a spectral method lets you turn it into an ode

then ???? then you have a system of first order ode?

[wide spear]

No not quite

Are you familiar with the ways of solving laplace's equation

[bright palm]

I guess if its not turned into a system im not sure what it means to turn a pde into a matrix

[bright palm]

youd have to be more specific

i havent seen anything but the vanilla wave equation in the context of numerical

and i guess i used some horrible method on burgers

or maybe someone can point me to a resource? I've tried reading the chapter (and this book in general) but I'm finding it incomprehensible

[wide spear]

Which book

[bright palm]

trefethen spectral methods

is that bad

[fathom rain]

see the hint

[thick onyx]

I'm trying to do Richardson Extrapolation from solving a differential-equation using a Runge-Kutta method (RK4). I have a hard time in general with understanding the reasoning behind Richardson Extrapolation, but I think I'm starting to get somewhat of an understanding of how to use it.

What I have a hard time understanding, is when to use which coefficients. My textbook some times uses E = 2E - E = 2A(h/2) - A(h) + O(h²) with E being the exact value, h being a step size and A(h) being an approximation with step size h. Now other times my textbook also likes to use 3E = 4E - E ≈ 4A(h/2) - A(h). From what I can read online it has to do with the truncation error of the approximation. How do you see the order of this error, and in turn which coefficients to use for the Richardson Extrapolation? In particular, which ones would I use for my differential-equation that has been solved using a Runge-Kutta method of 4th order?

[fathom rain]

just so you know, all spectral methods where you are interested in the eigenvalues fail spectacularly in my experience (you need high precision to make it good, so it is super slow)

but for just 6 eigenvalues it is probably fine

[warped heath]

Would you say that it is an implicit scheme because y_n+1 depends on f(y_n + y_n+1 / 2) is enough ?

[wide spear]

Yes

[bright palm]

ah it's not really critical to me whether it works or not I just want to understand what it's saying to do.

The hint seems to address more how to work the boundary condition into the matrix and not how to construct the matrix in the first place.

is it necessary to understand chebyshev differentiation matrices, and that should seem like the obvious choice of A_N here, or is it something else

[wide spear]

I think this should be in #real-complex-analysis

[warped heath]

It's in my numerical analysis course

but ok

thank you

[warped heath]

We have this cauchy problem :

and

how to show that the function X

satisfy that?

[warped heath]

someone ? 😢

[wide spear]

What have you tried?

[warped heath]

There is a method to solve this, I knew how to do a lot 3 years ago, but I forgot

I'm looking for my old notebooks...

[wide spear]

What is X'

[warped heath]

X'(t) = F(t, X(t))

[wide spear]

Ok but what is X'

[warped heath]

$ X'(t) = \begin{pmatrix}x'(t)\ x''(t)\end{pmatrix} = \begin{pmatrix} -x''(t) - tx(t)\ -x'(t) - tx(t)\end{pmatrix} $

[wide spear]

Ok and what does x'' equal

[pine jetty]

Pg

[wide spear]

Ok and what does F(t,X) equal

[warped heath]

they don't say that

[wide spear]

They do

[warped heath]

F(t, X(t)) is not F(t, (y z))

[wide spear]

X is a vector

[warped heath]

we have F ( t , ( x(t) x'(t))

not F (t , (y z))

[wide spear]

y=x(t), z=x'(t)

[warped heath]

$X'(t) = \begin{pmatrix}x'(t)\ x''(t)\end{pmatrix} = \begin{pmatrix} -x''(t) - tx(t)\ -x'(t) - tx(t)\end{pmatrix} = \begin{pmatrix} x'(t) \ -x'(t) - tx(t)\end{pmatrix} $

[wide spear]

Yes you already wrote this

[pine jetty]

Pg

[warped heath]

like that ?

[wide spear]

Yes

[warped heath]

what is important to see here?

[wide spear]

The point of this exercise is to see that you can turn an ODE of order>1 into a system of first order ODEs

[warped heath]

but it does not satisfy that

the first line ok, mais not the second line

[wide spear]

?

X does not satisfy the same ODE

[warped heath]

ok

thank you for you help!

[warped heath]

in fact the point is to show that any solution of

is solution of

and this allows us to deduce that there is a unique solution

but not sure to understand why this allows us to deduce that there is a unique solution

but not matter

[warped heath]

Hello

Someone Know about Newton method ?

[wide spear]

[warped heath]

"Write the recurrence relation of a sequence approaching the solution
of the implicit Euler scheme (with Newton's method)."

Correction :

I don't understand how they did that, on the internet I can't find anything about Newton's method applied to implicit Euler

[visual dune]

With an implicit scheme, you generally end up with a nonlinear equation (or system of equations) in the form of F(X) = 0, which you would solve using Newton's method. So, implicit Euler leads to the use of Newton's method.

[brave crypt]

how do I prove the rate of convergence and asymptote for a general iterative method $x_{n+1} = f(x_n)$

[pine jetty]

! matthewzz

[bold loom]

I am having problem with this example. I understood how to find tan(0.12) & tan(0.26) but for tan(0.40) & tan(0.50), I have no clue where to start.
I only have been taught the newtons interpolation formulas (forward & backward)
As it says in last ss that its extrapolation, which idk how to solve it, idk where to start. I tried looking up online & found nothing that was explaining this.
our uni has asked this specific question many times in past exams so it has been bugging me for a while.
Any help would be appreciated

[wooden kiln]

I have a Low-Storage 3rd-Order Runga-Kutta Method that I found in a old journal. I'm trying to reproduce it. For each step, we will evaluate the following, where $\frac{du}{dt} = f(u)$ (time derivative is autonomous).
Defining $u_n = u_0$, $u_{n+1} = u_3$:
$$ u_1 = u_0 + \Delta t \frac{8}{15} f(u_0) $$
$$ u_2 = u_1 + \Delta t \left( -\frac{17}{60} f(u_0) + \frac{5}{12} f(u_1) \right) $$
$$ u_3 = u_2 + \Delta t \left( \frac{5}{12} f(u_1) + \frac{3}{4} f(u_2) \right) $$
My goal is to create an equivalent Butcher Table for this system. How can I do so? I've already shown that it converges as well.

[pine jetty]

jacobjivanov

[wide spear]

What have you tried

Do you know how to turn a butcher tableau into the corresponding rk scheme

[wooden kiln]

I tried rewriting it completely out in one equation, as in this form. However, I don't believe it is correct because I think it would mean that my b2 value would be 0, which doesn't make sense to me

Also, please ping. Otherwise I might not see it

Secondarily, since I know it is autonomous, the c1, c2, ... values are irrelevent

[shut forum]

How would I get the system of equations from this

[fathom rain]

look for circulant matrices

[shut forum]

I got it now

Thanks though

[brave crypt]

How do i find the order of convergence for x_k+1 = (10/x_k)^(3/2)

Apparently i use something to do with atkinsons theorem

[night swan]

I’m trying to solve an IVP using spectral methods (MHD equations - A, v, t, rho). One of my main issues right now is that my initialization for my density is probably inconsistent with my initial vector potential (which I get from an LBVP using known current densities). I’m trying to write an NLBVP with the same physics as the IVP to get an initial density field that is consistent with that A, but I’m not sure what changes need to be made to the equations/formulation of the problem to make it the correct corresponding NLBVP. So far the only thing I know is that all time derivatives should be set to 0. Anyone know what else I need to change? For instance, should I take all terms with velocity to be 0, and are there equations that should be removed since I assume the vector potential is known, so there are fewer unknowns that need solving?

[wide spear]

Ummmm

For incompressible fluids there is a way to solve for the pressure from the velocity as an elliptic problem

Wait

Wait why don't you have pressure in any of your equations

Have you already used an equation of state to relate pressure and density

[night swan]

yeah, I think there are assumptions that pretty much remove pressure terms to just have temperature or density expressions

[wide spear]

Hmmmm

Ok so we want to find a balanced density state right

What if you set all the time derivatives to 0 and plug in all your initial conditions and an initial density

[bright palm]

Hmm could anyone recommend a good resource for understanding chebyshev matrices and boundary conditions? I have found trefethen spectral methods incomprehensible over the past week and am pretty lost on how to move forward, not finding a good resource through google. I am pretty new to numerical methods for differential equations.

I feel like trefethen treats it like the understanding of how to apply boundary conditions should just immediately follow from understanding how to construct the chebyshev differentiation matrices but I am not sharing this experience

[tight geode]

I'm finding conflicting info on the ability of the qr algorithm - if the input matrix has been converted to upper hessenberg, does the product of all the intermediate Q's always give the eigenvectors? If not, what's the next best thing?

[waxen hazel]

I Need solution some questions

[wide spear]

?

[wide spear]

What functions do people use for high order finite elements on triangles

[fathom rain]

I guess mostly Lagrange, Hermite, B-splines

[wide spear]

How do you do lagrange on triangles

Doing it on quadrilaterals is fine but the tensor product construction doesn't work on triangles

[fathom rain]

on the reference element and a mapping

[fathom rain]

If it is an unstructured triangular mesh you need to scale every integration with the area of the triangle. If it is a structured triangular mesh you can probably do a tensor product construction (by combining 2 triangles to a square). I have never tried

[wooden kiln]

We can use the streamfunction potential to remove pressure from the Navier-Stokes equations, but only in 2D

[wooden kiln]

If I have an array u[t, x, y] and am comparing the values in this array to an analytical solution for u(t, x, y), I can obtain a new array e[t, x, y]. What would be a good way to quantify the overall error in this array in order to prove convergence with a finer mesh? Would I be able to use a standard Ln-Norm, or would I have to use some type of mesh-dependent variation?

[wide spear]

That's not what they did, they used an equation of state to replace pressure by density

[wide spear]

Relative L^2 error is a pretty standard choice, depending on the true solution

[wooden kiln]

Would you be able to point me to a source for that? Additionally, I've previously been using the following, which my PI recommended, but I have no clue where it comes from.

[wide spear]

So this L^p is the absolute error

And is the discretized version of the L^p norm

[vernal narwhal]

Can someone point me to a resource or something that dicussed the local error analysis for (heun's method ) improved euler method. My book doesn't provide any analysis of it and It seems harder than applying taylor series unless I am missing something.

I followed along with this

but trying to apply it with the new method it seems I am stuck

I am not certain how to handle the f(x_n+1,y^*_n+1) term

the book just mentions that locally I should expect O(h^3) so I am somehow supposed to cancel up to the quadratic term in y(x) taylor series

[fossil halo]

This is a typical predictor-corrector method. https://en.wikipedia.org/wiki/Predictor–corrector_method#:~:text=In numerical analysis%2C predictor–corrector,satisfies a given differential equation.

Equation(4) is for predictor. You need to predict yn+1 first using Eq.(4) and correct it using Eq.(3).

I guess O(h^3) refers to the order of errors.

[idle phoenix]

how

[wide spear]

#advanced-pdes

Oh you posted there

Posting your question in multiple channels is discouraged

[vernal narwhal]

here it refers to the local error (true value vs computed after 1 step) order can how it is proportional to h "the step size". So if it linear (h), quadratic (h^2), cubic (h^2), etc. I ended up finding a video on computing the local error.

So here locally it happens to be h^3 so if I were to half my step size I would expect the error ~(h/2)^3=(1/8)h^3. So 1/8 of previous error

Anyways, I found the answer to my question. It seems you have to apply a taylor series of 2 variables and you can cancel the h and h^2 terms in the taylor series

🤷‍♂️

[waxen hazel]

i need solution question

[worn night]

what does it mean for a spline to be natural?

and what must be true about a spline for the endpoints?

[waxen hazel]

cubic spline

[worn night]

yes but natural spline has a specific meaning

[wide spear]

What have you tried

[waxen hazel]

??

[waxen hazel]

hello

[waxen hazel]

i need solution

[wide spear]

Please read #rules , no one will give you a solution

[fathom rain]

do your own homework

[waxen hazel]

i do not understand

[fathom rain]

then you get the grade you deserve.

nobody here will give you a solution. if you want to discuss and understand that is another thing

[waxen hazel]

I want to understand

i do not the solution direct

[wide spear]

Is there a way to precompute a svd for solving underconstrained least squares in lapack

Like how you can do dgetrf and dgetrs

[bright palm]

I gotta ask for help on the chebyshev thing again

still trying to do this problem, im really not sure how to set the boundary conditions

I asked my prof to help me out and he mostly talked about how to fix the solution, but we arent trying to find a solution, so im not sure

we also werent able to find an example in the text that was helpful with this problem

the program is like 5 lines of code so far, it just spits out really slowly converging eigenvalues plus a bunch of garbage

i just need resources or an explanation or something because i have not been able to find any resources for this, not even with my professors help-

m = 0; % periodicity in theta component

[D, r] = cheb(N); 
r = (r+1)/2; % transform to [0,1]
D2 = D^2; % second derivative

invr = 1./r; % these are bad, at this point
invr2 = 1./(r.^2); % both have Inf in them (at end)

M = diag(invr)*D + D2 - diag(invr2*m^2);
M = M(2:N, 2:N);

[~,l] = eig(M);
sqrt(-diag(l))```

this is literally all the code

output, for example, on N=7

   11.5295
   11.1846
    2.4046
    2.4306
    5.3966
    5.4003```

I'm not doing anything to enforce the boundary conditions here because i have no idea what to do, or even to what object I should be doing things too, since i can find no examples like this problem in the textbook or googling

🙏

[lunar widget]

Anyone here interested in entropy stable schemes?

[silk drum]

I think dgeqrf is to dgesvd as dgetrf is to dgetrs. https://netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga84fdf22a62b12ff364621e4713ce02f2.html

[wide spear]

dgesvd isn't for linear least squares though, dgelsd is

Or dgelss

[silk drum]

Looks like dgelqf via efficient path, A = L*Q.

https://netlib.org/lapack/explore-html/d7/d3b/group__double_g_esolve_ga94bd4a63a6dacf523e25ff617719f752.html

https://netlib.org/lapack/explore-html/d7/d3b/group__double_g_esolve_gaa6ed601d0622edcecb90de08d7a218ec.html

[wide spear]

My problem falls under path 2, not 2a

[vernal narwhal]

Oh you too are dealing with numerical analysis stuff

[bright palm]

yea, and the problem is still fucked

i cannot for the life of me figure out how to apply these boundary conditions

[vernal holly]

Can someone help me understand how to solve this for Gaussian quadrature weights etc

[pure geyser]

I also need some help on how to start this problem, do I need to do something with Jacobian matrix and eigenvalues?

[idle phoenix]

How did you learn the finite element program for magnetohydrodynamic equations

[wide spear]

Are you familiar with finite element for simpler problems

[idle phoenix]

I learned one-dimensional Poisson with Teacher He Xiaoming

I have programs written by others using fealpy, but just looking at the code makes me feel like I can't understand it

[wide spear]

Do you understand the theory

[idle phoenix]

understand

yes

May I ask what kind of programming do you use

I have to go to bed now. We'll talk tomorrow. Thank you for your reply

[subtle hemlock]

Is this the channel where I can ask about numerical linear algebra topics? I need help understanding how to code an LU decomposition of a tridiagonal matrix

[wide spear]

This channel is generally intended for questions related to theory

Programming questions should go to #computing-software

[glad fossil]

hello

[brave crypt]

i have no idea what equations to use

[prime flare]

haiiiiii

is anyone doing simpsons1/3 rule?

i need help to solve this.

[tight geode]

you just take the equation https://en.wikipedia.org/wiki/Simpson's_rule#Simpson's_1/3_rule, plug and chug for xi and xi+1 to get the integral of one step, and then add all the steps together to get the total
it's just tedious

[tight geode]

Is there a way to make gaussian elimination more robust? I'm using a linear algebra library and it sets a default error threshold of 10^-11. I'm working with a matrix that has numbers on order between 10^7 and 10^1. Trouble is, it scales the matrix down by its max value when solving for the eigenvectors, so it becomes 10^1 to 10^-7. When solving, it finds the rref, and one of the row operations produces a pivot value on order of 10^-17 which gets treated as a zero. If I ask for the rref of the original, non-scaled matrix, the pivots are within tolerance and it produces the correct rref.

[wide spear]

Which linear algebra library are you using

[tight geode]

this: https://github.com/markrogoyski/math-php

[wide spear]

Can you use lapack

[tight geode]

I wish but there's one person who wrote bindings for it 7 years ago

[wide spear]

https://pecl.php.net/package/lapack oh this?

[tight geode]

my mistake 12 years ago, but yes

[wide spear]

Can you call some C/C++/Fortran scripts from PHP

Like you write your matrix to a file, then you have some code that takes that file as input and does whatever you want, then writes the output to another file, which you then read in again

This is really borked but lapack is the gold standard for NLA

You could also write your own LU code

[tight geode]

Oh this thing does have LU decomp, I haven't tried substituting that in yet

[fathom rain]

What kind of elements do you have? Solve is symbolically or use some rational data type if possible

[tight geode]

just decimal numbers, i can't solve either of those two ways though. It's looking like i can divide the matrix by the smallest value to scale it, solve for the ref, and then multiply by the smallest value to avoid the loss. I think the margin of error becomes relative to the scale of the matrix values rather than absolute that way.

Actually do I need to re-multiply by the smallest value? Do the same principles for solving a linear system still apply? If so you'd get the same rref no matter how it's scaled.

[fathom rain]

why are you doing this in the first place?

and why not use rational represenations, it seems to be part of this math-php

[tight geode]

long story
i tried adding a pseudoinverse so i could solve a matrix system with a row of 0's in the A matrix. It worked until i tried solving a system with a particular 6x6 A matrix. The SVD given by the lib for this A matrix is wrong. After a lot of debugging, I've found that it broke both the eigenvector and rref methods in the library because it can't cope with the disparity in magnitude this particular matrix has.
I fixed the eigenvector thing, and now i just need to fix this. Then I can be on my merry way

It does but it needs some abstraction to get it in a usable state https://github.com/markrogoyski/math-php/pull/475#issuecomment-1804241830

I did try scaling by the lowest element and it just kinda worked

[fathom rain]

you wont get away from limitations of floating point arithmetics

[tight geode]

it's a self-imposed limitation, the default error tolerance is 10^-11 which is fine. I only need like 10^-3 honestly in my final answer
The issue is that a matrix with numbers spanning just two orders of magnitude (10^1 -> 10^3) leads to some matrices with elements on the order of 10^-17 that get zeroed out opaquely.

[fathom rain]

why do you even consider doing this by gaussian elimination? and why php? it is just bizarre

[tight geode]

solving for force/moment equilibrium on a gearbox using timeseries data in a database managed by a backend written with symfony 🤷
the lib internally is using gaussian elimination to find the rref, not me

[fathom rain]

then somebody with no clue has implemented it...

[tight geode]

would you use the LU decomposition instead?

[wide spear]

Gaussian elimination is LU

[fathom rain]

if the matrix is ill-conditioned I would use high precision, symbolic (since you said it was just 6x6) or rationals. lu is one choice yes

you solve this system many times? doe the matrix change? does the matrix have structure?

[wide spear]

Wait if you're trying to use a SVD why is it computing the LU decomp?

[tight geode]

the matrix has structure, i could technically write it out by hand but that would create a mess
it goes pseudoinverse -> SVD -> SVD computes eigenvectors of MTM and MMT -> eigenvector method uses qr alg to get eigenvalues and then solves underdetermined system using an algorithm involving rref

[fathom rain]

what is MTM and MTM ...I hope you are not computing svd using the normal matrix...

[tight geode]

MT = M transpose, I'm not good at latex; here: https://github.com/markrogoyski/math-php/blob/ea4f212732c333c62123c6f733edfb735a4e3abd/src/LinearAlgebra/Decomposition/SVD.php#L95

[fathom rain]

that is a bad idea since you square the condition number...

[tight geode]

If you can point me to a better way of computing it I'd be happy to try it
I'm doing the best I can with what I got, and what I got is this library and whatever's left in my brain from that one numerical methods class I took 6 years ago

[wide spear]

https://en.wikipedia.org/wiki/Singular_value_decomposition#Numerical_approach

[fathom rain]

see also slide 12 https://cims.nyu.edu/~donev/Teaching/NMI-Fall2014/Lecture5.handout.pdf

[tight geode]

thanks! i'll look into it if i hit a wall

[silk drum]

Is this a structural problem? Is the matrix dense or sparse?

[cerulean rover]

Hey can someone help me with a question?
If A is matrix that is symmetrically positive definite,
If you use the cholesky method for determining decomposition A=LL^T , does it fail in some instances where you don't use pivoting? I thought the cholesky method always worked for symmetrical positive definite matrices

[wide spear]

I think it should work without pivoting unless you're running into floating point errors

Also pivoting would unsymmetrize the matrix wouldn't it

[cerulean rover]

Yea

So I'm confused rn

Because I have this question

One of my answers here is incorrect, and I know for a fact that by my textbooks definition that the cholesky method is O(N^3), and I tested out the first two statements and I found that the diagonal was consistently positive, so I was thinking the only possible last thing would be the case of pivoting

[wide spear]

The screenshot is a dox

[cerulean rover]

Better?

[silk drum]

LAPACK allows for complete/full pivoting so it's likely beneficial numerically in some cases. Partial pivoting would introduce a symmetry issue but obviously not complete pivoting. https://netlib.org/lapack/explore-html/da/dba/group__double_o_t_h_e_rcomputational_ga31cdc13a7f4ad687f4aefebff870e1cc.html

Stability of the computation due to round-off addressed/mentioned in https://en.wikipedia.org/wiki/Cholesky_decomposition

[wide spear]

You need pivoting for when the matrix is only semidefinite

[tight geode]

definitely sparse, this is the structure
it's the sum of forces and moments about an x-axis aligned shaft that solves for the reactions of a bearing pair. The x/y/z are the offsets of the bearings on the shafts. The x's are the only non-zero offsets.

[wide spear]

You can do a block SVD then

[wide spear]

Compare with https://www.netlib.org/lapack/explore-html/d1/d7a/group__double_p_ocomputational_ga2f55f604a6003d03b5cd4a0adcfb74d6.html

[silk drum]

It's worth considering (iterative) spare solvers rather than a direct solver. If the matrix is ill-conditioned, preconditioners can be considered. The PETSc library has historically been very good, and easy to use, and can provide inspiration even if you need to use another library or strategy. PETSc is especially convenient for large problems but is still a nice starting point for small problems.

[wide spear]

Writing a petsc interface for php sounds like a nightmare

[waxen hazel]

what is meaning Hermite poly?

[fathom rain]

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

[idle phoenix]

Why does proving this indicate that the format is stable

[long lake]

Sorry we can't allow discussion of piracy on this server

[cerulean rover]

What are the arithmetic complexity (big O notation) for a general invertible matrix using the 1. Overall step method (Jacobi) 2. Single step method (gauss-seidel) 3. Sor procedure

I can't find anything in my textbooks about them, I assume they are probably relatively simple but idk

[wide spear]

Trefethen and Bau or Demmel should include these

[cerulean rover]

Thanks

[cerulean rover]

Pain

[wide spear]

Is this Trefethen and Bau

[cerulean rover]

Yes

Looking through demmel rn

I found the answers too sorta

These are what it's supposed to be

Now the thing is I know my question doesn't want the space column because last time I answered O(n) and it marked it wrong

But there isn't any fractional

Just n^1 2 3 or 4

So it's serial time

Also that is for reccomending demmel

Honestly the best book I've read on numerical linear algebra

My professors book is all mumbo jumbo

Filled with definitions and nothing else

[wide spear]

Space is how much memory it needs, not the runtime scaling

[cerulean rover]

Can someone clarify this for me, are they all a possibility, none or what

[cerulean rover]

Also this one, it says widest range of convergence, so if

1. |C-1|<3 was -3<c-1<3
2. |C+1|<3 was -3<c+1<3
   And 3. |c| was -2<c<2

If either 1 or 2 was "wider" than 3, the answer should be one of those two, correct? And in that case it would be probably 1? Or am I mistaken

[idle quartz]

Lets say I have a large matrix A and i can solve the system Ax=b.

I have another system
transpose(B) *A^-1 * B u = v
How can I solve this when B is not a square matrix?

I don't have an explicit form of A^-1. I m just using conjugate gradient for the first system

[wide spear]

Is A sparse

What do you know about B

[idle quartz]

A is sparse

Nothing other than B is rectangular

The dimensions match of course

[wide spear]

How far from square is B

How fast do you need to solve it

[idle quartz]

Not necessarily fast

Like how can I manipulate my systems to do this?

It s not like I can say let Bu = w and
Now I can write B^t M^-1 w =v and then say M^-1 w = s and now I can write
B^t s = v
I could solve that and back solve and substitute but B^t isn't square to invert

[wide spear]

Computer the SVD of B, then you also know the SVD and pseudoinverse of B^T, left multiply by the pseudo inverse, solve for Ax=B^-Tv for x=Bu, then solve for u using the SVD

[idle quartz]

There has to be another way tho right?

B is sparse too

[silk drum]

How large is B?

[wide spear]

What sort of way do you want

[crystal verge]

QR-decomp?

[dusky junco]

Does the following equation have any real solutions: bxc+axd +ayc-byd=0

[dense echo]

x=0, y=0, at least.

[brave crypt]

hi
Someone here recommends a good book about Sobolev spaces?

[wide spear]

#book-recommendations

[edgy flax]

i am trying to find the hessian matrix of $\lVert Au-Bv\rVert_2^2$ respect to u and v

[pine jetty]

亜城木 夢叶

[edgy flax]

so far, I have $\begin{bmatrix}A^TA & -A^TB\ -B^TA & B^TB\end{bmatrix}$

[pine jetty]

亜城木 夢叶

[edgy flax]

it is not symmetric

[wide spear]

What do you have for the first derivatives?

Also not sure what you mean with respect to u and v

[plucky kayak]

looks symmetric to me

B^T * A = (A^T * B)^T

and the matrix ((A,B), (B^T, C)) is symmetric provided A and C are

[mortal pier]

$\text{How do you choose a suitable function to approximate} \frac{1 - e^{-\frac{x^3}{64}}}{3x^3} \text{so that it can be numerically integrated from 0 to infinity? I was thinking of expanding it into a Taylor polynomial and then truncating the tail, but it seems like I would have to expand it to a very high-degree polynomial.}$

[pine jetty]

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

[wide spear]

Please fix your latex

[fathom rain]

julia> quadgk(x->(1 .-exp.(-x.^3/64))./(3*x.^3),big"0",Inf)
(0.02790560973653903785058013480181955879300718102038854576688064921304477674929857, 1.14255890558531101826176060193660020211371754083048451332118798824097631927917e-40)```

[dense glen]

When using the Euler Maclaurin summation formula, to what accuracy do you need to calculate the Bernoulli numbers? Seems to be the main overhead when using the formula when I'm using it

[fossil glade]

sin(0.5x) can't be expressed as any finite trigonometric polynomial, can it?

[wide spear]

Yes, it has an infinite fourier series, only with cosine terms though

[wide spear]

What accuracy do you want your sum to be?

[cerulean rover]

isnt the second and 5th answer the same?

[wide spear]

A^TA and AA^T are not the same

[cerulean rover]

care to explain?

[wide spear]

What do you want me to explain

Why they're different?

[cerulean rover]

yes

[wide spear]

Think about the SVD A=UDV^T

[cerulean rover]

the eigen vectors and eigen values are ordered differently?

[fathom rain]

example: ```julia
julia> A
2×2 Matrix{Int64}:
1 1
0 1

julia> A*A'
2×2 Matrix{Int64}:
2 1
1 1

julia> A'*A
2×2 Matrix{Int64}:
1 1
1 2```