17

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


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.


6

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

Using the half angle formulas you can convert the exponent into the form $$a + b \cos 2\theta + c \sin 2\theta$$ which then integrates nicely into Bessel functions. Mathematica gives $$\int_0^{2\pi}e^{a + b \cos 2\theta + c \sin 2\theta} d\theta = e^a \pi I_0\left(\sqrt{b^2 + c^2}\right)$$


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


4

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

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

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}}}^{...


3

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


3

You should definitely use CFIE over EFIE or MFIE, due to the problem of internal resonances. Basically, if you are at/near a frequency where there is a resonance in one of your scattering objects, then the EFIE or MFIE matrix may become singular/ill-conditioned. Using a linear combination of them greatly reduces the chances of singularity. I don't know much ...


2

In the framework of mathematical physics, the fundamental solution is the response of an infinite domain to a point source. E.g. in electrostatics the electric potential field $\varphi$ satisfies the Poisson equation \begin{equation} \nabla^2 \varphi = -\frac{\rho_f}{\epsilon_0} \end{equation} where $\epsilon_0$ is permittivity and $\rho_f$ is the free ...


2

As I understand it, collocation method for partial differential equations is something akin to interpolation. First we characterize the solution space as a linear combination of some set of linearly independent functions $\phi_i(x)$. The appropriate choice of $\phi_i(x)$ depends on the problem. But assuming that the anticipated solution is sufficiently ...


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

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


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


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


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