9

This is a linear integral equation (because the right hand side is linear in $h$) and presumably the problem is ill-posed (for one, because of the multiplication by $\sin(\theta)$ but also because the "integral kernel" is smooth). Methods to solve this (approximately) are general Galerkin methods or collocation methods. For Galerkin methods you take a ...


6

Seems that you have a duplicate eigenvalue. Thus, you have two eigenpairs $(\lambda_1, x_1)$ and $(\lambda_2, x_2)$ where $\lambda_1 = \lambda_2$. Denote $\lambda = \lambda_1 = \lambda_2$. Let $\alpha$ and $\beta$ be arbitrary complex numbers. Then $$A (\alpha x_1 + \beta x_2) = \alpha A x_1 + \beta A x_2 = \alpha \lambda x_1 + \beta \lambda x_2 = \lambda (\...


6

Preamble: "integers" -> Integers in your problem statement, Mathematica uses CamelType. The method depends on your problem, naturally. In your case it appears to be the Contejean–Devie method. A [buried] excerpt from the Mathematica documentation: For polynomial systems, Reduce uses cylindrical algebraic decomposition for real domains and Gröbner basis ...


6

Actuallly you can scale the variable of the problem into a different interval. Transform $ x \in [0,a]$ where ($a=10^6$) to $\xi \in [0,1]$ by the formula $ \xi = \frac{x}{a} $. Thus, using the chain rule we have $\frac{d y }{ d x } = \frac{d y}{ d \xi } \frac{d \xi }{ d x } = \frac{ 1}{ a } \frac{d y}{ d \xi } $. Hence, replacing $\frac{d y }{ d x }$ by $\...


5

MATLAB always uses the LAPACK libraries to calculate eigenvectors which works on double precision, floating point numbers - MATLAB's default data type. Mathematica's method depends on its input type. For example, when you do TestMatrix = {{1, 2, 3}, {3, 1, 2}, {2, 3, 1}} Eigenvectors[TestMatrix] You'll get an exact answer involving Sqrt[3] and so on. ...


5

By default, Mathematica assumes that all variables are complex numbers, and when working in the set of complex numbers, Sqrt[-21] is well defined. You can tell Mathematica (version 8) that you are working on the set of reals using Solve[Sqrt[2 x - 9] == Sqrt[4 x + 3], x, Reals] which gives no solution. For versions earlier than 8, you need to use Reduce ...


5

Why implement it by hand? Matlab, Maple and Mathematica all have tools builtin to solve differential equations numerically, and they use far better methods than you could implement yourself in finite time. In Matlab, you want to look at ode45. In Maple it's called dsolve (with the 'numeric' option set), in Mathematica it is NDSolve.


5

You could use Runge-Kutta method to solve this system numerically, first rewrite your second order equation as a first order system by doing following substitution trick: $$ \left\{ \begin{aligned} x_1' &= x_2 \\ y_1' & = y_2 \\ x_2' &= x_1 y_2 + g \cos y_1 -(M+m)g/m \\ y_2' &= -g(\sin y_1) /x_1 - 2x_2 y_2/x_1 \end{aligned} \right. $$ Now ...


5

I want to mention another idea in case it helps. The truncation error of replacing $\int_{-\infty}^{\infty}$ with $\int_{-L}^L$ is on the order of $$\begin{aligned} \int_L^\infty \cos(tx)g(x)\,\mathrm{d}x &= \frac{1}{t}\sin(tx)g(x)\big|_{L}^{\infty} - \int_L^\infty\frac{\sin tx}{t}g'(x)\,\mathrm{d}x \\&= -\frac1L g(L)\sin(tL) + \text{asymptotically ...


5

Note: I'm somewhat worried at this point that the integral values Mathematica gives me are bogus. I thought it was working because it gave a sensible-looking result in a short time, but it might be the case that the method it tries to use is buggy or that I did something wrong. So it might be that the code below isn't working at all, I don't know, sorry. ...


4

@tqviet already says this in a some way, but let me phrase it differently: $10^6$ is not a large end time by itself. It seems like a large number because it has 6 zeros at the end, but it all depends on (i) which units you compute it, (ii) how fast the effects are that you are trying to model. To give examples: If I told you that you had to compute to $10^...


4

Julia's macros and string macro are made for doing this kind of thing. You can define a domain-specific language (DSL) within a macro scope and use that as a way to define problems. One example of this is the DifferentialEquations.jl @ode_def macro. You may have to use a string macro (due to parsability into ASTs) but this sort of tool is what allows you to ...


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

The short answer is "It depends". Certainly, you can try to find the Helmholtz decomposition on your sampled data, and find your irrotational and solenoidal components. However, there are certain requirements on your original vector field you started with. In general, it is that the vector field you are trying to decompose has to be sufficiently smooth and ...


3

Are you asking for a plot of $z(d)$ and $w(d)$? Wolfram Alpha might not like the long query, but Mathematica simplifies your number/letter soup just fine: $$w(d)=-\sqrt{3}\frac{ d^2 (x+2 y-1)-2 d \left(x (2 y-1)+y^2\right)+x (4 y-3)-2 y^2+6 y-3}{d^2-2 d y+4 y^2-6 y+3} = \frac{1}{2} \left(\sqrt{3} d-3 \sqrt{(d-1)^2}+\sqrt{3}\right)$$ $$z(d)=\frac{x \left(d^...


2

In Mathematica I recommend do this with Sparse Matrix-Vector products as follows. First I will define a circulant matrix $C$ (with the periodic BCs) which appears on the RHS of your equation. In matrix notation your equation reads \begin{equation} x(n+1) = x(n) + \epsilon C x(n) + \sqrt{\epsilon} \eta(n) \end{equation} or \begin{equation} x(n+1) = (I + \...


2

Create some mesh in $(x,t)$-space, take out the elements containing the origin, and impose nonzero boundary conditions on the resulting boundary that you believe is appropriate for the solution near the singularity. Then solve the problem without the forcing term in a domain with the resulting hole. This will give you an idea how the solution will look like. ...


2

To begin, I would probably advice to consider the steady problem by looking for ground states $\phi = \Phi(x)e^{i\omega t}$. In this way you reduce your equation to the nonlinear Helmholtz-type problem: the elliptic operator $\omega\phi + \phi_{xx}$ + nonlinear part $\lambda\phi^3$. Now, I don't see what you can do theoretically with the nonlinear term. ...


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

Try using the Reduce equation while requiring x to belong to the Reals. See the Mathematica documentation for an example.


2

Seeing what are your boundary conditions, you should use an spectral(Galerkin) method. I am trying to solve sometrhing similar, but in 9 dimentions, And the way to treat the field fading in infinity as a boundary condition is with galerkin methods. Probably using a basis like $F_k(x)=e^{-\alpha x^2}\sin(k x)$. Though you should look more into the probable ...


2

Here is something I threw together, it is definitely incorrect, as the program is throwing errors. In particular, a SolutionVariableNumberError is raised. Likely this is because the equations are written incorrectly. The likely suspect, I'm guessing is in the SourceTerm's. This is a start though for how one might go about solving these equations in FiPy. ...


2

There are two possible issues going on here. Firstly as noted in the answer by knl any linear combination of eigenvectors corresponding to degenerate eigenvalues is also an eigenvector, as shown in that answer. However even for the non-degenerate case eigenvectors may vary between different invocations of a diagonaliser even if the input matrix is the same. ...


1

You got good advices in comments, but as it is only 1D problem, you can manage it using fully explicit method with very small time steps. I checked your code in Mathematica and there are several misunderstandings there. After correcting them, the code works for me. Similarly, your description in the question contains, I suppose, some mistakes. As I guess ...


1

If you want to diagonalize significantly larger matrices on large HPC clusters, I would look at Elemental, which is a state-of-the-art distributed dense linear algebra library. It's a better alternative than other libraries out there (e.g., ScaLAPACK, PLAPACK), and builds upon BLAS and LAPACK implementations.


1

The classic book in B-Splines is de Boor, C. (1978). A Practical Guide to Splines. New York: Springer-Verlag, but I spent many weeks on it and was not able to make his algorithms work. The algorithm in Phillips, G. M., & Taylor, P. J. (1996). Splines and Other Approximations Chapter 6 of Theory and Applications of Numerical Analysis (2nd ed., pp. 131-...


1

How about the C++ (with python bindings) Gamma library for NMR simulations: The Gamma Library. In the GammaDetailedDescription, the authors mention the various components of the library, which in turn would solve many aspects of your question.


1

You could use pdepe and turn Laplace's equation into a BVP, then solve it with a multiple shooting method. This solution is not a turnkey solution, but it does use MATLAB built-in functions. Alternately, you could set up the equations you've written in matrix form; these are algebraic and you should be able to solve them using an LU decomposition (if they'...


1

Given the apparent dearth of quasi Monte Carlo integrators for C++ (or C), I wrote my own implementation to be used with GSL. It's not especially well-tested, nor does it implement Mathematica's algorithm, but it should be better than nothing.


1

Runge-Kutta methods such as (4,5) are also available in the GNU Scientific Library (which is written in C). They also include adaptive time-stepping.


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