# 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 ...

9

I don't have experience with Octrees, but whenever there is some nice C++ library that I want to use in Fortran, I simply write a simple C driver --- typically a few C functions that do exactly what I need. Then I call them from Fortran using the iso_c_binding module. This has the great advantage that you reuse a well tested library with a community around ...

7

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}}\... 6 You seem to be very set on using Fortran. Octrees, when implemented efficiently, are rather complex data structures and, as such, better suited to programming languages that have more support for this, such as C/C++. There are a number of very high quality implementations in C/C++ that you could use. 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-... 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 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. 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 First, it's always good to further debug the issue and seeing exactly where (what term, for what parameters) the overflow is coming from. This was not clear to me from the question. Once you know precisely what the problem is, you can better diagnose the issue. For example, lets take the overflow issue. This is somewhere above 1e300 if I remember ... 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 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 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 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 ... 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 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. 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

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 ...

1

Look into Anderson mixing. It is used for self consistent quantum chemistry applications which are nonlinear. It is faster and more robust than simple damping. The headline is that you keep a history of n previous iterations and calculate the residual for each, then perform a linear solve on the set of residuals using the approximation that the solution ...

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