5

Please find below an implementation of the Runge-Kutta 2 method in Python for your problem. This takes care of integrating the differential equation from $-1$ to $5$ for a given value of $y'(-1)=-1$ (fixed) and $y(-1)$ (changing). As you can see, setting $y(-1) = 1$ gives $y'(5) \approx 1.34$ with this method and for $y(-1) = 2$ it yields $y'(5) \approx -5....


5

The short answer is that you need $$\phi_{-1} = \phi_0$$ $$\phi_N = \phi_{N-1}$$ to impose $\nabla\phi=0$. A quick check by making the following change if idx == -1: idx = 0 elif idx == N: idx = N-1 in the code, you have posted shows that the average $\phi$ remains constant up to 14 decimal places. To see why this is the correct boundary condition ...


4

You're approaching this the wrong way by using the product rule of differentiation. Rather, making use of the fact that $$ \frac{\partial u(x_i)}{\partial x} \approx \frac{u_{i+1/2}-u_{i-1/2}}{\Delta x}, $$ the first step in finding what you need to do is to recognize that $$ z(x_i)\frac{\partial u(x_i)}{\partial x} \approx z_i \frac{u_{i+1/2}-u_{i-1/2}...


3

Expanding (sort of) on @MPIchael's answer, you can pick any smooth function you like and plug it into the heat equation to give a problem to then work the other way. In numerical methods, we call this the Method of Manufactured Solutions, and it is used extensively for verifying computer programs designed to simulate PDEs. You'll have to add a forcing ...


3

Dont expand the divergence term, apply integration by parts without doing this $$ \int_\Omega w \nabla\cdot(\alpha\nabla u) dx = - \int_\Omega \alpha \nabla u \cdot \nabla w dx + \int_{\partial\Omega} w \alpha \frac{\partial u}{\partial n}ds $$ In your notation, define a vector field $$ \vec{A} = \alpha \nabla u $$ and then do the integration by parts on $$ \...


2

Here's a minimal solution using Forward Euler method for integration and dichotomy to find $y(-1)$. I found that $y(-1)\approx 1.1926$ which is the same value that I got when I used Mathematica. from math import *; import matplotlib.pyplot as plt def solve(a): #Solve the IVP y''(x)=f(y(x),x) with the initial conditions y(-1)=a and y'(-1)=1 using ...


2

As David said, absorbing boundary conditions won't be completely reflectionless. That said, we can reduce relfections quite a bit, which helps to avoid influence from the boundaries while the particle still travelling inside. Since this is a time dependent problem, one simple choice of boundary conditions will look like this. At the left boundary: $$\frac{\...


2

Your problem was the lower limit of integration. It should have been $-x_e$ instead of 0, since $x_e$ is the equilibrium point for the potential and not the minimum distance. After correcting that, you get the following #%% Solution xe, lam = 1.0, 6.0 # parameters for potential xmax = 10 # Bval, Bval2 = wavefunction values at x = bound1, bound2 bound1, ...


2

It depends on the equation you are solving. For example, consider the following heat equation with a source term: $$ \frac{dT}{dt} - \nabla^2 T = S$$ If $S \neq 0$, then you have a source term that changes the value of the temperature as time evolves. Consequently, even if your initial condition is $T(\textbf{x},t)=0$, your temperature will increase because ...


2

TL DR: $$u_1(x_1) = \cos(2\pi~(\frac{x_1}{L_1}) - \pi) + 1$$ $$u_2(x_2) = \cos(2\pi~(\frac{x_2}{L_2}) - \pi) + 1$$ $$u(x,t) = \exp(-a t) u_1(x_1) u_2(x_1)$$ How to construct it: Sines and cosines are easily differentiable so they make a good starting point to construct such a solution. We chose a section and offset of the cosine which has a derivative of ...


2

Bill Greene presents the point of view of how things have been done traditionally. The "modern" way is to add "constraints" to your linear system which, in the current case, would be $$ u_{3x} = -u_{3y}. $$ These constraints can be entered straight into the linear system without the detour of the augmented linear system you show -- ...


2

The traditional way to handle this type of boundary condition is to create a transformed coordinate system at node 3 so that one of its axes is along the direction you want to constrain. You use the coordinate transformation matrix at this node to transform your stiffness matrix and then constrain the system in the usual way, e.g. by removing the equation ...


2

Let's take the example of the unsteady one-dimensional heat equation inside a solid on a domain $x\in[0,1]$: $$\partial_t u - D\partial_{xx} u = 0$$ with the initial profil $u(0,x) = u_0(x)$ at $t=0$, and Dirchlet boundary conditions enforcing that the wall at $x=0$ (respectively $x=1$) is at temperature $u_{L}$ (respectively $u_{R}$). If we use discretise ...


1

Have you plotted the solution using a visualization package or just a cut through the middle of your domain? Sometimes just looking at the solution you're generating will give you an insight into where it's arising from. Based on your description, I'd say there's a mistake in your implementation of your boundary condition. ETA: The last sentence is probably ...


1

You simply make sure that the initial condition satisfies the boundary conditions and that you don't add anything to the elements of $u_{n+1}$ corresponding to the boundary. In other words, you drop the rows and columns of $M$ that correspond to the boundary nodes and solve only in the interior nodes. Let me add that inhomogeneous Dirichlet boundary ...


1

As the mathematica.se thread shows, the solution of $$ \begin{aligned}\frac{\partial}{\partial x}\left( \operatorname{sign}(x) u(x) \right) + \frac{\partial}{\partial x} \left( \sqrt{u(x)} \frac{\partial u}{\partial x}(x) \right) &= 0 & &\text{in } \Omega = (-6,6), \\ u &= 0 & &\text{on } \partial \Omega = \{-6,6\} \end{aligned}$$ is ...


1

It depends on the implementation. More than the connectivity, you need to focus on the orientation of the normals. The common convention is to define the connectivity such that the normals on the boundary elements point away from the solid. This is relatively easy in 2D since there are only two nodes for one (linear) element. For 3D, choose the connectivity ...


1

It's just a linear 1D problem, you can easily do implicit time stepping here so that numerical stability would not be a problem. The accuracy should not be an issue either since for such a simple problem you should be able to use any spatial resolution you need. More specifically, let the equation be discretized in space by any scheme, e.g., low-order ...


1

The boundary conditions don't depend on the choice of your basis but on the formulation you have for your problem. If you have a "standard" finite element formulation, you don't need to do anything to apply (homogeneous) Neumann boundary conditions, they are already satisfied by your system. In the most common formulation, Neumann boundary conditions are ...


1

There's a really nice book by Li and Ito on the immersed interface method, which was designed to solve problems like yours. Chapter 2 describes 1D problems. Basically, the finite difference method works well when the coefficients are smooth, but when they're rapidly varying or have discontinuities things can go to hell very quickly. You can make it work by ...


1

It looks to me like you want to look at the step-62 tutorial program of deal.II, written by Daniel Garcia. It covers pretty much exactly the application you have in mind.


1

Maybe this isn't a helpful response but the reason this happens for the matrix form of Laplacians is because this actually happens for the true infinite-dimensional Laplacians in some settings. In recangular coordinates, $N$-D Laplacians act exactly like their discretized counterparts in this way. For example, let $\Delta_3$ by the Laplacian on $[0,1]^3$ ...


1

This is my attempt at providing some intuition. Everything I state might be obvious, moreover it doesn't have much to do with physics, so this could be a non-answer. I will ignore boundary conditions. Imagine we have a 2D grid with function values $u_{ij}$ for $i=1\ldots m, j=1\ldots n$. Let the $x$ axis point downward, the $y$ axis to the right. N.b. $z$ ...


1

I have also had this problem and spent a lot of time on various forums and I have finally come up with a good solution. This solution will allow you to specify a boundary condition on any whole subdomain within your geometry, and each subdomain can have its own separate boundary condition. I am posting this solution here as the old official Fenics QA forum ...


1

This seems to be more of a Physics question than a Computational Science one. Due to the symmetry of your problem, you can conclude that the solution is of the form $$\mathbf{u} = u_r \hat{\mathbf{e}}_r(r)\, ,$$ since the selection of the zenithal and azimuthal angles is arbitrary. This turns the PDE system \begin{equation} (\lambda + \mu) \...


1

Essentially going one step back in the derivation of the wave equation you take the divergence of the (simplified) Euler equation $\mathbf{\nabla}.(\rho_{0} \frac{d\mathbf{v}}{dt}=-\nabla p$) and next, you let $\rho_{0}$ out of the operator, assuming it is a constant. Now at the inferface between two media with $\rho_{01} \neq \rho_{02}$, $\nabla \rho_{0}$ ...


1

What you find is indeed correct. It is known that positivity is lost if very small time steps are chosen, see https://doi.org/10.1515/cmam-2015-0018 This loss of positivity happens even for semi-discrete scheme. The analysis for 1-d case is given in section 6 of this paper.


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