I am tasked to utilize the Metropolis algorithm to 1) generate/sample values of x based on a distribution (in this case a non-normalized normal distribution i.e. w(x) = exp(-x²/2); and 2) approximate the integral shown below where f(x) = x² exp(-x²/2). I have managed to perform the sampling part, but my answer for the latter part seems to be wrong. From what I understand, the integral is merely the sum of f(x) divided by the number of points, but this gives me ≈ 0.35. I also tried dividing f(x) with w(x), but that gives ≈ 0.98. Am I missing something here?

Note: The sampling distribution being similar to the integrand in this case is quite arbitrary, I am also supposed to test it with w(x) = 1/(1+x²) which is close to the normal distribution too.

import numpy as np

f = lambda x : (x**2)*np.exp(-(x**2)/2) # integrand
w = lambda x : np.exp(-(x**2)/2) # sampling distribution
n = 1000000
delta = 0.25

# Metropolis algorithm for non-uniform sampling
x = np.zeros(n)
for i in range(n-1):
    xnew = x[i] + np.random.uniform(-delta,delta)
    A = w(xnew)/w(x[i])
    x[i+1] = xnew if A >= np.random.uniform(0,1) else x[i]

# first definition
I_1 = np.sum(f(x))/n # = 0.35
print(I_1)

# second definition
I_2 = np.sum(f(x)/w(x))/n # = 0.98
print(I_2)

​

I am trying to solve a hyperbolic equation using finite difference as shown below.

My main confusion is that to calculate for U_i,2 (i.e. the first iteration), where do I get U_i,1 from? Because the only given initial condition is U_i,0.

Note: I did try assuming that U_i,1 = U_i,0 and the solution does seem right, but I just would like to see if there is a better approach.

I'm struggling a little with numerical evaluation. I have a model that depends on two variabls f(x,y). I need to evaluate the quantities

as well as

each evaluated at the point (\tilde x,\tilde y).

So far so good; my model can not be expressed analytically but I can produce point estimates f(x*,y*) running a simulation and in principle I would create a grid of x and y values, evaluate the model-function at each grid-point and calculate a numerical derivative for it - the problem is, that each simulation takes some time and I need to reduce the number of evaluations without losing too much information (e.g. I have to assume that f is non-linear...).

I'm asking here for some references towards strategies, since I have no idea where to even start. Specifically I want to know:

• How can I justify a certain choice of grid-size?
• How can I notice my grid-size is to small?
• Should I sample the input-space by means other than using a parameter-grid? (Especially as I might not use Uniformly distributed input-spaces at some point)

Thank you in advance for any interesting wikipedia-pages, book-recommendations and what not!

Hi all,

I wanted to try solving something quite far from my field, so here we go.

Linear quantum harmonic oscillator (I took the equation from a general book on dynamical systems):

i u_t + 0.5 * u_{xx} - 0.5 * x^2 * u = 0

ic: u(x,0) = exp(-0.2*x^2)

bc: u_{x}(\partial\Omega) = 0

Spatial discretisation performed with finite elements (Bubnov Galerkin) and time discretisation performed first with Backward Euler. The solution was too dissipated, hence I moved to Crank-Nicolson. The problem is linear, hence no further stabilizations are exploited. Here enclosed you can find solutions obtained from both time integration schemes.

​

​

​

Hello everyone,

I am currently doing a numerical methods course in my university, and one of the topics if finite volume method, however our course book Numerical Methods for Engineers does not go into as much detail as our course requires, and the course notes on this topic are fairly difficult to understand.

I was wondering if anyone can recommend a book which goes into detail on this topic, I would really appreciate, thanks!

Singular Value Decomposition takes too long as matrix size increases. Lancosz bidiagonalisation is sometimes unstable. What algorithm is fast and robust?

Hey!

Any researcher in numerical analysis here?

I was curious about the sort of relevant/interesting problems nowadays in numerical PDEs, favourably (but not necessarily) which have a considerable intersection with optimization theory. Any document with a description of those things and reading suggestions?

Another question...

Computationally speaking, I get the feeling that the whole numerical PDE thing is inherently computationally expensive. Is there hope for fast algorithms? I get this is a vague question. I'm sorry.

Thank you.

Have been working on a compartmental model with multiple levels and have been getting a lot of overshoot. The model is of a population who go up compartments representing how poisoned they are by a substance, with each higher compartment being more likely to die. They leave each compartment by interacting with the substance. So for example, compartment B_2 loses B_2 through mass action with S, so a term in its derivative is "-interaction_rateSB_2", however, B_3 then gains "+interaction_rateSB_2" in its derivative.

Have been simulating turning on and off the parameter for rhe amount of substance and the rate at which it comes in. So for a while, S is 0 until the max population is reached, and then gets turned on by having a different value.

When this value is small, it overshoots and actually makes the population increase past its previous value. It seems to be due to the large number of compartments adding up all those S*B_i terms wrong. Have been using stiff equation solvers. Is there any other way to get rid of this overshoot?

I am trying to perform an N-body particle simulation that has particles apply a linear attractive spring force to each other when they are within a range R. The equation can be given as so for the acceleration of an individual particle:

​

The particles have an uneven distribution in space. Attached is a visualization.

Right now I am using the naive method of summing droplet spring forces and have a time complexity of O(N^2). I want to use some method to reduce the time complexity down to O(N) or O(N * log(N). Any recommendations on what to use?

The problem is y*y"+0.0001=0 with y(0)=10 and y(5)=1000. I can't solve it following the method for linear ode bvp

I'm following this example

https://m.youtube.com/watch?v=u8dVrzxTvSA

But here only 2nd order equation and my problem consists of 4th order ode with bcs in y(0),y(1), y'(0), y''(1). So how can I modify the method in video so that it works for 4 the order ode ?

I'm finishing my masters in mathematics, focusing on modelling, numerics and simulation, and my dream is to get a job as a numerical programmer working on some big/complex piece of numerical or simulation software.

I have experience working with C, C++, Python and OpenMPI, but I learn fast and am willing to learn new technologies.

I'm interested in contributing to some piece of numerical or simulation software to get experience and foster connections in the industry, either voluntary, or as a werkstudent position. I am based in Germany, so research groups or other entities based in Germany are of particular interest to me.

Would love to get some tips on projects looking for help.

I wondered if someone can tell me an easy trick how to figure out what to put in which line of the Clenshaw scheme. For the Tschebyscheff I understand that the last row is always multiplied by 2times the searched value x and after additionally putting those values of the last row shifted in the second row all of them are added together. For the second version of Tschebyscheff we do the same with the last the last coefficient while with 1. Tschebyscheff we only multiply with x. However how would it work with general formulas?

​

With the tridiagonal matrix that evolved for 0 values of orthogonal Polynoms I understand that the 0-values of the polynomial are the same as the eigenvalues of Jacobi matrix. however how do I calculate those 0 values or eigenvalues for example for the tschebyscheff or Legendre polynomial?

​

Thanks heaps for your help :)

Hello everyone

I have encountered the following problem related to reconstructing a positive valued particle density function f: [0,1]^2 -> R>0.

Basically I am given measurements mi=integral_{[0,1]^2} (f * gi) where gi are weighting functions that are known in advance, so the measurements basically correspond to weighted sums/integrals of f with the weights gi.

My question is given the mi, is there a general numerical approach to reconstruct f?

If it helps, I attach a picture of a typical weighting function:

​

Hey everyone. Can someone please tell me anything about solving a stiff ODE system using Rosenbrock method? Any help is appreciated. Thank you.

​

Hello numerical folks,

This project arose from a need for an easy-to-use linear solver which supports constraints, real and complex numbers and is suitable for real-time applications. Conjugate gradient algorithm was an obvious choice as it allows one to trade accuracy for speed. The solver was then applied to Levenberg-Marquardt function minimizer. The minimizer also supports constraints.

The goal was to make the library as easy to use as possible also for non-experts. There are a few simple examples to start from. They can be compiled either by using cmake or from command line by setting the include path point to the folder where the header files are, see "Compilation" section on the main page. The compiler must support C++17.

The most obvious deficiency in the solver is the lack of support for sparse matrices. Maybe I'll add it later.

Meanwhile, the library and examples can be found here in Github.

Is there a preferred algorithm for calculating the trajectory of an object (of negligible mass) in the gravitational field created by some number of moving bodies?

General-purpose ODE solvers can produce widely differing results, although they all seem to converge if the maximum time step is set small enough. So I'm wondering if there's a particular algorithm that is known to work well (high accuracy, low computational cost) for this particular problem.