4
$\begingroup$

I am trying to solve the following system of differential equations numerically over the domain $x=0$ to $x=D$. The main difficulty is that the boundary conditions are interconnected and depend on the solution variable. I'm not too sure how to approach the finding a numerical solution for that reason.

$$ \frac{dI_1}{dx} = -\alpha(x) I_1(x) \\ -\frac{dI_2}{dx} = -\alpha(x) I_2(x) $$

subject to the interconnected boundary conditions,

$$ I_1(0) = AI_2(0) + C\\ I_2(D) = BI_1(D) $$

where $A$, $B$, $C$ are known constants.

Background

This is an optics problem and describes the multiple reflection of light entering and exiting a dielectric layer with surface reflection coefficients $A$ and $B$, absorption coefficient $\alpha$, $C=(1 - A)I_0$ where $I_0$ is the incident intensity (a known constant).

I have divided the propagation into two streams which propagate downwards and upwards. This is a common approach for solving the radiative transfer equation. The $I_1$ equation is the intensity of light propagating to the right, and $I_2$ is the intensity propagating backwards to the left. So the total intensity is $I(x) = I_1(x) + I_2(x)$.

Two-steams approximation

Analytical solution when $\alpha(x)$ is constant

If anyone is interested there is an analytical solution in the limit that $\alpha(x)\rightarrow\alpha$,

$$ I(x) = I_{1}\left(x\right)+I_{2}\left(x\right)=-\frac{{\left(Be^{\left(2\,\alpha x\right)}+e^{\left(2\,D\alpha\right)}\right)}Ce^{\left(-\alpha x\right)}}{AB-e^{\left(2\,D\alpha\right)}} $$

when $A=B=0$ this reduces to the Beer-Lambert law as expected,

$$ I(x) = I_0e^{-\alpha x} $$

This is a useful expression because this way you don't have to calculate and keep track of the magnitude of the multiple reflections.

Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ These look like Robin boundary conditions for a two-point boundary value problem governing a second-order ODE. I would be surprised if you could not do a multiple shooting method on this problem. $\endgroup$
    – Geoff Oxberry
    Feb 18, 2014 at 4:09
  • $\begingroup$ That's interesting. So can I transform this into a different form in which the boundary conditions are easier to apply? $\endgroup$
    – boyfarrell
    Feb 18, 2014 at 5:18
  • $\begingroup$ No, I'd leave it as it is; I'm just noting a vague similarity. You have a boundary value problem, and as I said, I would be surprised if you couldn't apply standard methods to your problem. That said, I would use the integrating factor approach instead. $\endgroup$
    – Geoff Oxberry
    Feb 18, 2014 at 18:25

2 Answers 2

Reset to default
3
$\begingroup$

Assuming that $\alpha(x)$ has a closed form, I think you can use an integrating factor to integrate this exactly. Even if not, you should be able to write it in a form that's amenable to numeric integration (quadrature) rather than discretization with finite differences.

Improve this answer
$\endgroup$
6
  • $\begingroup$ $\alpha(x)$ is a material parameter so the complexity of the function depends on the situation so in general there is not a closed form. However in practice 99% percent of the time $\alpha(x)$ will be a piecewise constant function over domain. Can we use the integration factor approach in this case? $\endgroup$
    – boyfarrell
    Feb 18, 2014 at 5:15
  • $\begingroup$ Almost certainly. The thing to be integrated will be $e^{\int \alpha(x) dx}$ which can be easily broken and integrated over the pieces. Given that, I suspect that this problem has an analytic solution even for the case you describe. $\endgroup$
    – Bill Barth
    Feb 18, 2014 at 13:10
  • $\begingroup$ Even if the integral isn't closed form, all you need is the numerical value $\int_0^D a(x) dx$. $\endgroup$
    – Geoffrey Irving
    Feb 18, 2014 at 17:22
  • $\begingroup$ @GeoffreyIrving, that's not quite right. With an integrating factor approach, you need, essentially, $\int_0^x \alpha(t) dt$ for all the $x$ values you care about, so you need to be able to integrate over arbitrary domains. In order to nail down the constants of integration, your single definite integral is all that's needed, but in order to plot or evaluate $I_1$ and $I_2$ at arbitrary points, you need to be able to do arbitrary definite integrals. $\endgroup$
    – Bill Barth
    Feb 18, 2014 at 17:55
  • $\begingroup$ Sorry, I was a bit unclear. The issue is solving multiple interconnected ODEs. If you do the single $0$ to $D$ integral, you can reduce one ODE to function evaluation and solved the coupled ODE problem using a nonlinear solver. Plotting the final result does require more than one value, but that's only after the coupled problem is solved. $\endgroup$
    – Geoffrey Irving
    Feb 18, 2014 at 17:58
1
$\begingroup$

You said that for the following problem

$$ \frac{dI_1}{dx} = -\alpha(x) I_1(x) \\ -\frac{dI_2}{dx} = -\alpha(x) I_2(x) $$

with boundary conditions

$$ I_1(0)=AI_2(0)+C\\ I_2(D)=BI_1(D) $$

the coefficient $\alpha(x)$ is piecewise constant. As I can se you can do it in an iterative way. You start with an initial guess for $I_1(0) = I_{1,0}$, and what you can do is iterate in time (with finite differences or with discontinuous Galerkin to have stable solution), until you reach $I_{1,n}\approx I_1(D)$. Then you use that value to calculate the inital guess for $I_{2,n}$ with the boundary condition

$$ I_{2,n} \approx I_2(D) = B I_1(D) $$

and you integrate in the oposite direction, reaching the position $0$ with $I_{2,0}$. With this value as an approximation for $I_{2}(0)$, you update the value for $I_{1}(0)$ usgin the boundary condition at $0$. You then iterate again to reach position $D$, use the approximation for the initial guess for $I_{2,n}$, and go againt until $0$. You perform this recursive scheme until you converge.

I think that the scheme is not garanted to converge, because it depends on the initial guess. My experience is that these kind of schemes usually have good behaviour.

EDIT:

You can also rewrite your problem as follows:

$$ \frac{dI}{dx} = -\overline{\alpha(x)} I(x) $$

with boundary conditions

$$ I(0)=\frac{A}{B} I(2D)+C\\ $$

where the new function $I$ and the new coefficient $\overline{\alpha}$ are defined over the domain $[0,2D]$ as follows

$$ I_{[0,D]}(x) = I_1(x) \\ I_{[D,2D]}(x) = \frac{1}{B} I_2(-x) $$

with boundary conditions

$$ \overline{\alpha}(x)_{[0,D]} = \alpha(x) \\ \overline{\alpha}(x)_{[D,2D]} = \alpha(-x) $$

Then your function $I$ will be continuous over the domain $[0,2D]$ (it is continuous at $x=D$ because the different parts satisfy the boundary condition $\lim_{x\to D-}I(x) = I_1(D) = I_2(D)/B = \lim_{x\to D+ }$, and you still have the other boundary condition conecting the points $0$ and $2D$. So you can integrate in time and refeed the system, or you can build the whole system with for example discontinuous galerkin method.

Remark: If the coefficient function $\alpha$ is piecewise constant as you said, I dont think that you can solve it analyticaly. Only numerical solutions are possible.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ Many interesting suggestions, I need think about this and find the best way to solve. I basically want to be able to specify a structure and reliably calculate $I(x)$. $\endgroup$
    – boyfarrell
    Feb 18, 2014 at 23:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.