7

The issues you're running into now are not a failing of Newton-Raphson, but a question of coupling. You're doing iterated sequential coupling -- solving each equation sequentially and then iterating until (hopeful) convergence. No solver choice in place of NR is going to fix this lack of convergence, as long as you are doing iterated sequential coupling. ...


4

I'm assuming the limitation in 2D and 3D is storing the Jacobian. One option is to retain the time derivatives and use an explicit "pseudo" time-stepping to iterate to steady state. Normally the CFL number you need for diffusive and reactive systems might get prohibitively small. You could try nonlinear multigrid (it's also called Full Approximation Storage ...


3

You will need to write your problem such that the unknowns are a single vector, not a matrix. In your example with $N=2$, you will have an unknown vector $x(t)$ of size $20\times 1$ (not a matrix of $10\times 2$). You will solve a problem of the following shape $$\dot{x}(t) = A(t)\, x(t), \mathrm{ with }\ x \in \mathbb{R}^{n},\ A(t) \in \mathbb{R}^{n\times n}...


3

Yes You should see that there is at least one zero (to machine precision) eigenvalue of that system. The eigenvalue problem can be solved because it is prepared to find the zero eigenvalues. The problem with $A$ and $x$ is different from the problem with $K$ and $u$. $x$ will transport whatever Neumann-ness it has to $b$, but $u$ will satisfy Dirichlet ...


2

As mentioned by Matt Knepley, this is naturally formulated as a system of partial differential algebraic equations. Because you're in Matlab, you could consider doing the spatial discretization yourself (e.g. finite difference, finite volume, finite element) to obtain a system of DAE's, then use the method of lines to step forward in time. You can discretize ...


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

Odespy needs from you a 1st order system of ODEs, namely something of the form \begin{align} \frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y},t), \tag{1} \end{align} where, for your particular problem, $y_i$ are single-variate complex-valued functions that depend solely on time (that means you need to arrive at the semi-discrete form first by discretizing the ...


1

After some digging I could finally find a paper that more or less answered my question: Bermúdez, A., Rodrıguez, R., & Santamarina, D. (2003). Finite element approximation of a displacement formulation for time-domain elastoacoustic vibrations. Journal of computational and applied mathematics, 152(1-2), 17-34. It turns out that the equation for the ...


1

Given the dependency between $z$ and $t$, we can rewrite the first equation with a time derivative: $$\partial_tU=\frac{c}{n}\left(\nabla_r^2U+\rho U\right)$$. Now we define $g=\left(\begin{array}{c}U\\\rho\end{array}\right)$. This allows us to rewrite the system of equations as: $$\partial_t g=\left(\begin{array}{cc}\frac{c}{n}\nabla_r^2&0\\0&0\...


1

The first step is to add one more equation defined as: $$\frac{\partial y_3}{\partial t} = v$$ Then if you substitute equations 1 and 2 and the new equation into your modified equation 3 you get: $$ \frac{\partial v}{\partial t}=(1-y_1)v + \frac{\partial^2 y_3}{\partial z^2} + \frac{\partial y_3}{\partial z} + \frac{\partial^2 y_2}{\partial z^2} + \frac{\...


1

As written, it's not clear to me that it's actually fully coupled. That is, the solution to $B$ depends (via its BC2) on the solution to $a$, but assuming $\omega$ is simply some specified function of position and time, it looks like $a$ doesn't depend on the solution to $B$. If that's the case, then simply solve for $a$ on a 1D mesh first, then you have the ...


1

From my experience with Navier-Stokes equations, one can do very well without fully implicit schemes. If you just want a fast numerical scheme for the solution of the time evolution, take a look at IMEX (implicit-explicit) schemes; see e.g. this paper by Ascher, Ruuth, Spiteri Implicit-Explicit Runge-Kutta Methods for Time-Dependent Partial Differential ...


1

Here is a guess as to what might be going on. ode45 (like all the MATLAB ode solvers) adjusts the step size based on it's estimate of the solution error. All equations are considered in this calculation whether they are coupled or not. The solution errors are compared to two parameters, AbsTol and RelTol. Both of these have default values but sometimes ...


1

From what I understand from your question, the first plot has become unstable, so perhaps if you did a plot with the two methods for the first time period (so up to around 4 seconds) your two solutions will match (looking at your plots they seem to match). If that was the case, then try solving the problem with a different solver, or making your time steps ...


1

Instead of running Newton-Raphson to convergence of each subsystem, try 1 iteration on the first subsystem followed by 1 iteration on the second subsystem. This may keep the subsystems more coupled, not going to "distant" unrelated sub solutions. Repeat this two step iteration and see if it converges.


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