# Tag Info

18

There seems to be quite a bit of confusion about how to apply multi-step (e.g. Runge-Kutta) methods to 2nd or higher order ODEs or systems of ODEs. The process is very simple once you understand it, but perhaps not obvious without a good explanation. The following method is the one I find simplest. In your case, the differential equation you would like to ...

8

In my opinion, it is not a good neither a bad mesh. It clearly depends on the PDE you are considering. The finite space to which the PDE is projected is your mesh, where your operators, e.g. $\vec{\textrm{grad}}$ (gradients), $\textrm{div}$ (divergences), $\triangle$ (Laplacians)... strongly depend on that mesh and become matrices: $$\vec{\textrm{grad}}\... 5 Quadrature points x_1,x_2 are also unknown. You need two more equations corresponding to x^2 and x^3. Then you solve for a_1, a_2, x_1, x_2. 5 A common strategy is to employ a damping strategy, i.e., to compute \vec{w}^{\ast} = F \left( \vec{w}^{(n)} \right) and then set \vec{w}^{(n+1)} = \alpha \vec{w}^{(\ast)} + (1-\alpha)\vec{w}^{(n)} where \alpha\in[0,1]. You typically choose \alpha in such a way that it minimizes \|[\alpha \vec{w}^{(\ast)} + (1-\alpha)\vec{w}^{(n)}] - F(\alpha \vec{w}^... 5 You haven't specified boundary conditions, but it sounds like you want to solve the Cauchy problem (i.e., the initial value problem on the whole real line \infty < z < \infty). In that case, there is no need to use finite differences. You can simply use Fourier analysis (separation of variables) to write down the exact answer. Once you have that, ... 5 A classic paper for evaluation of the integrals commonly present in computational electromagnetics (EM) is: D. R. Wilton, S. M. Rao, A. W. Glisson, D. H. Schaubert, O. M. Al-Bundak, and C. M. Butler, "Potential integrals for uniform and linear source distributions of polygonal and polyhedral domains," IEEE Trans. Antennas Propag., vol. AP-32, no. 3, pp. 276-... 4 First of all, it is impossible to completely rewrite NIntegrate because usually, it evaluates (parts of) the integrand symbolically to check for certain properties. This means, you probably have to re-implement Mathematica completely. Nevertheless, if you are tackling one specific integral, I assume chances are very good that you can get similar numerical ... 3 I think that both meshes are good. But depending of the problem at hand, one might be better than the other. In the domain that you are mentioning, one thing to consider is that you probably don't want element sizes to be too different to satisfy the sampling criterion. Regarding mesh quality I would suggest that you check this post: Commonly-used metrics ... 3 Let's first rewrite this. From your formulas, you have that$$ A_{t+1} = A_t^\rho \exp(e_t) $$where A_t is just the previous value and, consequently, just a fixed parameter. So what you need to compute is then$$ E_t[\phi(A_{t+1},\eta_{t+1})|A_t] = \frac 1C \int_{-\infty}^\infty \int_{-\infty}^\infty \phi(A_t^\rho e^a,b) e^{-a^2/(2\...

3

Your description "the function $h(y,x)$ gives the proportion of the mass $u(y,t)$ at position $y$ that moves to position $x$ in space in a small unit of time" seems to indicate a slight misunderstanding of derivatives. Since the term on the left of your equation ($\partial_t u$) is a time derivative, the "small unit of time" must be infinitesimally small. ...

3

I am turning my comments into this answer. First of all, the analytic expression that you proposed does not satisfy the differential equation. According to Maple, the solution is $$f \left( x \right) ={\frac {k_{{2}}\sqrt {\mu}}{\sqrt {a{k_{{2}}}^{2}- a+\mu}}{\rm sn} \left( {\frac { \left( \sqrt {2\mu - 2a}\ x+k_{{1}} \right) \sqrt {\mu}}{\sqrt {a{k_{{2}}}^{... 2 The question asked is too broad, in my opinion. Assessment of BIE methods vs. finite difference methods (or other domain methods) require careful analysis of many different points. Among them is D \subset \mathbb{R}^2 or D \subset \mathbb{R}^3 or even D \subset \mathbb{R}^N with N>3? is D finite or infinite? is D simply connected? which is ... 2 The decay properties of the Green's function depend, among other things, on the coefficients in your equations. For example, lower-dimensional wave guides typically transport information for long distances, and you may have to add up a lot of terms for good accuracy. Alternatives are to write the Green's function as a sum of sines and cosines that already ... 2 Well, it seems that there is a size requirement to get a good resolution I'm not sure if that is the properly way to solve this, but it worked for me. What I did was to interpolate the x and y components of the vector field. Talking about Python, this can be done using the scipy.ndimage.zoom function. As an example, I had this vector field (100 x 100): ... 2 LIC is very good for visualizing such fields, but what is the texture warped here? First of all, use LIC with white noise in the very beginning. You could then try Perlin noise and similar. There are also other variants of information, such as the coherence directions, which would give you a better idea of the flow field. You could check out this work of ... 2 You can use scipy.integrate with a given method (Simpson's rule or composite trapezoidal rule for instance): from scipy import integrate x = array1 y = array2 int = integrate.cumtrapz(y, x, initial = 0) See here for more. 2 Like what @Kirill says, write a script that defines f(u) to be a function that approximates the integral. Thus all that you need to is solve f(u)=0 which can use standard root finding tools like Newton's method. You're best choice for this is probably Julia. Using the built in quadgk function for Gauss-Kronrad quadrature, you can define the f function ... 2 I think that you can use convolve() from scipy.signal. As mentioned in a previous question, you can take advantage that the Fourier Transform of a convolution represents a product. 2 You can usually solve these kinds of equations via a transformation. Shampine discusses how Volterra integral equations can be transformed into an ODE which is then solved with a stiff ODE solver. If you discretize u(x) into a system of ODEs first then you can maybe do something similar in that case. If you want to handle a general functional ODEs directly, ... 2 The problem is that you are integrating an oscillatory function over an infinite interval. The MATLAB website doesn't give specifics on the algorithm behind their integral function (it just says 'globally adaptive'). The older quad function used to be adaptive Simpson, so I can assume that integral is the same. What probably is happening is that the code ... 1 I didn't manage to find a suitable formulation for the diagonal. But doing the convolution of the column then the convolution of the rows was enough computation time gain for me h_vec = gaussienne(X) convol_row = scipy.signal.convolve(p, h_vec, 'same').T adj = np.diag(scipy.signal.convolve(convol_row, h_vec, 'same')) 1 There is no need to use numerical methods. The integral has a simple analytical solution once you know the trick. Hint: First solve$$ \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} e^{-(x^2 + y^2)} dx dy $$by change of variables. 1 That depends on what kind of integral transform you are looking at. Your comments suggest that you are looking at a Fourier transform specifically, so I would recommend the FFT implementation of NumPy. Using this approach you can also tackle Laplace transforms. For the most general case you will have to evaluate your convolution using a brute force ... 1 In general, there are three major method families that are used (excluding locally corrected Nyström, as I am not an expert on this one): Treat the singularity. In this Q/A, I describe a canonical example in electromagnetics community with an EFIE integral, where for a singular integral, the singular part is treated analytically, while the rest is computed ... 1 Let's try some alternative avoiding the forward and backward differentiation of the np.gradient function. You need to apply the boundary conditions in each computation of the time derivative, as that determines the derivative at the boundary points. Thus you have to use there either$$ \rho_{zz}(x_0,t)\approx 2\rho[x_{-1},x_0,x_1]=\frac{ρ(x_{-1},t)-2ρ(x_0,t)+...

1

I'm not gonna edit your code to make it work but based on your code, I think you have some misunderstanding about finite difference method. Your PDE equation is: $$\frac{\partial \rho}{\partial t} = D \frac{\partial^{2} \rho}{\partial z^{2}} + G - B(N_{D} \rho + \rho^{2})$$ Your boundary conditions are: $$-D \frac{\partial \rho}{\partial z} \Bigg |_{z=0} =... 1 If a homogeneous Neumann problem is considered, i.e. f=0 in \Omega and \mathbf{n}\cdot\nabla\phi on \Gamma, one solution is given by the constant function$$ \phi\left(\mathbf{x}\right)=\phi_0\,. $$Then the boundary integral equation reads$$\int\limits_{\Omega}\phi\left(\mathbf{x}' \right)\delta\left(\mathbf{x}'-\mathbf{x}\right)\mathrm{d}V' =\...

1

You can use odeint, which is part of the Boost library, here's a long list of examples odeint

1

Though I have never implemented this myself, the conventional wisdom I've heard is that the spatial summation converges very slowly and it's preferable to use the Ewald transform, e.g. http://w3.uniroma1.it/lovat/giampiero/Documentipdf/IJ24.pdf There's probably boatloads of papers on this subject in IEEE transactions.

1

When you call q = integral(@(r)LO_scheme2a(r,kk), 0, inf), the integral function choosesr`, which is a vector that discretizes the integration bounds. The documentation states: For scalar-valued problems, the function y = fun(x) must accept a vector argument, x, and return a vector result, y. This generally means that fun must use array operators instead ...

Only top voted, non community-wiki answers of a minimum length are eligible