Imposing normalization in numerical solution of of ODE

[jorgdv]

Hello all,

I would like to know how to impose a normalization condition to numerically solving an ODE. For simplicity let's consider the example
\frac{dy}{dx}=y
You could use different methods using an initial value, but if you consider the interval [x_0,x_1] and \int_{x_0}^{x_1} y(x)dx=1 How would you impose that condition in a numerical scheme (instead of the initial value condition)?
This arises in the context of PDEs over PDFs, such as Fokker-Planck equations. I consider the ODE case for simplicity.

Thank you very much

 

[Delta2]

It must be a diabolical coincidence but in this very example the integral condition is equivalent to y(x_1)-y(x_0)=1

But,well i don't know a generic method to impose such conditions into the numerical solution process.

 

[the_wolfman]

It depends on the problem. In general you can introduce a tunable parameter, and then iterate on the parameter unit the integral constraint is satisfied.

This particular ODE is invariant under the change of variables y = \alpha y^*. However, the constraint is not invariant: \int y dx = \int \alpha y^* dx = 1. This allows us to solve for \alpha using \alpha= 1 / \int y^* dx.

Solving this problem then takes three steps. First, solve the ODE y^* = \frac{dy^*}{dx} for some guess initial condition. Use something simple like y^*\left(x_0\right) = 1.

Second, integrate \int y^* dx to find \alpha.

Finally, use y = \alpha y^* to find y.

 

[Strum]

If you solve in Fourier space you will only have to make sure that your zero component of the transformed data is preserved in your numerical scheme.

 

[jorgdv]

Thank you all for your responses. The analytic approach of invariant equations under change of variables is very interesting, and would work in this case, however I am looking for more generic approaches (if they exist). And about preserving your zero component in Fourier space, wouldn't it only preserve an additive constant in the solution?

The thing is that this solution would be a stationary solution of a PDE, say for example
\frac{\partial y}{\partial t}=\frac{\partial f(x)y(x)}{\partial x} -K\frac{\partial ^2 y(x)}{\partial x^2}
Then the stationary solution would be something like the example I posted for f(x)=1 applying the condition lim_{x \to \infty}y(x)=0, so if it was not stationary you normally could use a Finite Elements or Finite Differences scheme, so a question would also be: how would you implement the normalization condition in the FE and FD approach to preserve the normalization across time?

Thank you all

 

[Strum]

Well the I did assume that the inteterval in the first post was the boundary of your system. Perhaps that was not the case?
$lim_{x\rightarrow \infty} y(x) = 0$ probably isn't the easiest boundary condition to implent in a finite simulation. Have you considered transforming x such that infinity is moved to a finite value?

 

[jorgdv]

Sorry, it was not the case, maybe I did not explain myself well. There's no boundary to the system. In fact the normalization condition would be \int_{-\infty}^{\infty} y(x) dx=1, but for purpose of computation it could be aproximated with \int_{x_0}^{x_1} y(x) dx with an interval that's big enough, that's why the condition \lim_{|x| \to \infty} y(x)=0 (yes, there's also the simetric limit condition) helps, since you could find better approximations by widening the interval.

But that conditions comes to place mainly when reducing the previous equation to a stationary sistem in the way
0=\frac{\partial}{\partial x} (f(x)y(x)) -K \frac{\partial ^2 y(x)}{\partial x^2}
Solving the first derivative means
f(x)y(x)-K\frac{\partial y(x)}{\partial x}=A
But since that should be valid for all x, then A should be 0 to satisfy the limit condition, and we are left with a similar system to the first I posted.
But that is only to give a context, I am intereseted in the question of implementing the approximated normalization condition in a numerical scheme.

Thanks

 

[the_wolfman]

First, as a general comment. The study of differential equations is a huge field. There is no general approach that allows one to solve all differential equations. Likewise there are no generic approaches that allow you to treat the integral constraint for all differential equations.

how would you implement the normalization condition in the FE and FD approach to preserve the normalization across time?

I'm not sure that I'm following you. But this is not a well-posed problem. If you give me an initial condition and two boundary conditions then you've completely specified the PDE. By then requiring that the integral over the domain has to be equal to unity, you are over-constraining the problem. There is no guarantee that a solution exists.

I am looking for more generic approaches

My statement "you can introduce a tunable parameter, and then iterate on the parameter unit the integral constraint is satisfied" is a very general approach. I then showed you how to use it with the ODE that you supplied.

In general consider a 2nd order ODE with 1 boundary condition and 1 integral condition. I can solve this by using the supplied boundary condition, and and guessing the second boundary condition. For instances I could suppose that the second boundary condition is y(x_2)=\alpha. Then the trick is to iteratively adjust \alpha until the integral constrain is satisfied.

Another approach is to reformulate the differential equation into a minimization problem (this works particular well for finite element methods). You can then use a penalty methods. The idea is to introduce a penalty function of the form c \left(\int y dx -1\right) where c is a large positive parameter. Since c is large, the minimization will converge on a solution that approximately enforces \left(\int y dx -1\right) = 0.

 

[jorgdv]

Okay, thank you!