I have the following elliptic PDE (describing temperature in a plate, w in thermal contact with two fluids h and c): $$\lambda_h \frac{\partial^2 \theta_w}{\partial x^2} + \lambda_c V \frac{\partial^2 \theta_w}{\partial y^2} - (\beta_h+V\beta_c)\theta_w+(\beta_h\theta_h+V\beta_c\theta_c) = 0 \tag 1$$

pde = Subscript[λ, h] D[θ[x, y], x, x] + Subscript[λ, c] V D[θ[x, y], y, y] - (Subscript[β, h] + V Subscript[β, c] [Theta][x, y] + (Subscript[β, h] θh[x, y] + V Subscript[β, c] θc[x, y]) == 0

with the boundary conditions ($\beta_h, \beta_c, V, \lambda_h, \lambda_c$ are constants)

$$\theta_w(0,y)=1, \theta_w(x,0)=0$$ $$\frac{\partial \theta_w(1,y)}{\partial x}=\frac{\partial \theta_w(x,1)}{\partial y}=0$$

where \begin{eqnarray} \theta_h(x,y) &=& \beta_h e^{-\beta_h x} \int e^{\beta_h x} \theta_w(x,y) \, \mathrm{d}x,\\ \theta_c(x,y) &=& \beta_c e^{-\beta_c y} \int e^{\beta_c y} \theta_w(x,y) \, \mathrm{d}y. \end{eqnarray}

It is also known that $\theta_h(0,y)=1$ and $\theta_c(x,0)=0$.

Trying with separation of variables:

Substituting $\theta_h$ and $\theta_c$ in $(1)$ and using the ansatz $\theta_w(x,y) = e^{-\beta_h x} f(x) e^{-\beta_c y} g(y)$ we obtain two linear third-order ODEs with constant coefficients for $F(x) := \int f(x) \, \mathrm{d}x$ and $G(y) := \int g(y) \, \mathrm{d}y$: \begin{eqnarray} \lambda_h F''' - 2 \lambda_h \beta_h F'' + \left( (\lambda_h \beta_h - 1) \beta_h - \mu \right) F' + \beta_h^2 F &=& 0, \tag 2\\ V \lambda_c G''' - 2 V \lambda_c \beta_c G'' + \left( (\lambda_c \beta_c - 1) V \beta_c + \mu \right) G' + V \beta_c^2 G &=& 0,\tag 3 \end{eqnarray} with some separation constant $\mu \in \mathbb{R}$.

The boundary conditions take the form:

For G: $G'(0)=0, G(0)=0$ and $\frac{G''(1)}{G'(1)}=\beta_c$

In order to get another homogeneous condition of $F$, we use $\theta_{h1}(x,y):=\theta_h (x,y)-1$. This does not change Eq. $(2)$ or $(3)$ but gives us another condition on $F$ i.e. $F'(0)=0$

For F: $\frac{F''(1)}{F'(1)}=\beta_h, F'(0)=0$

Can anyone suggest a way forward in this problem ? I know I need to find the Eigenvalues $\mu$ and use the non-homogeneous b.c. and orthogonality to complete the solution, Krazug has already provided a method to find the $\mu$ values in an answer below. I need to understand how to use these with the $F$ equation to complete the $\theta_w$ solution

CONTEXT The integro-differential equation from which the above problem arises is $D$:

I had a system of three PDEs $$\frac{\partial \theta_h}{\partial x}+\beta_h (\theta_h-\theta_w) = 0 \tag A$$

$$\frac{\partial \theta_c}{\partial y} + \beta_c (\theta_c-\theta_w) = 0 \tag B$$

$$\lambda_h \frac{\partial^2 \theta_w}{\partial x^2} + \lambda_c V\frac{\partial^2 \theta_w}{\partial y^2}-\frac{\partial \theta_h}{\partial x} - V\frac{\partial \theta_c}{\partial y} = 0 \tag C$$

On eliminating $\theta_h$ and $\theta_c$ from $C$ I reach $$\lambda_h \frac{\partial^2 \theta_w}{\partial x^2} + \lambda_c V \frac{\partial^2 \theta_w}{\partial y^2} +( -\beta_h - V \beta_c )\theta_w +\beta_h^2 e^{-\beta_h x} \int e^{\beta_h x} \theta_w(x,y) \mathrm{d}x + \beta_c^2 e^{-\beta_c y}\int e^{\beta_c y} \theta_w(x,y)\mathrm{d}y = 0 \tag D$$

The BCs are

$$\theta_w(0,y)=1, \theta_w(x,0)=0$$ $$\frac{\partial \theta_w(1,y)}{\partial x}=\frac{\partial \theta_w(x,1)}{\partial y}=0$$ and $\theta_h(0,y)=1$, $\theta_c(x,0)=0$.

  • 2
    $\begingroup$ Do $\lambda_h$ and $\beta_h$ have values? And what are the other boundary conditions on $F$? $\endgroup$
    – SPPearce
    May 14, 2020 at 16:08
  • $\begingroup$ @KraZug $\lambda_h-\lambda_c$ and $\beta_h=\beta_c$. The other boundary conditions on $F$ could not be determined because the two other conditions left i.e $\theta_w(0,y)=1$ and $\theta_h(0,y)=1$ are non-homogeneous in nature. I was hoping that once the eigenvalues are determined from the homogeneous $G$ equation, maybe they could be used to determine the function $F$. That is how I have seen PDE equations being solved in textbooks. But am I missing something here ? Can you suggest some alternative approach for this problem ? $\endgroup$
    – Avrana
    May 14, 2020 at 17:29
  • $\begingroup$ Edit: $\lambda_h=\lambda_c$ $\endgroup$
    – Avrana
    May 14, 2020 at 18:15
  • $\begingroup$ Oh, I think I understand a bit more of what you are trying to do. I didn't realise you had an integro-differential equation. What values for mu do you currently get, because when I evaluate your code then I FindRoot errors regarding Accuracy/Precision, and it just returns the original start point without searching for a root. $\endgroup$
    – SPPearce
    May 14, 2020 at 19:32
  • 1
    $\begingroup$ I've put a quick semi-answer to one of your questions. It by no means solves all your problems. If you knew another two BCs conditions on $F$ you could include them in this method, to solve for both sets simultaneously. $\endgroup$
    – SPPearce
    May 14, 2020 at 20:09

1 Answer 1


This is not a full answer to the question!

I have a package which helps solve eigenvalue BVPs by calculating the Evans function, an analytic function whose roots correspond to the eigenvalues. Some details are available at these two questions, or this PDF. Or search for CompoundMatrixMethod to see my previous answers here. It comes much better than the determinant method for finding roots accurately.

Install the package (also available on my github page):

    "Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"]

Load the package and setup the system:


eqn = {V λc G'''[x] - 2 V λc βc G''[x] + ((λc βc - 1) V βc + μ) G'[x] + V βc^2 G[x] == 0};
bcs = {G[0] == 0, G'[0] == 0, G''[1] == βc G'[1]};
sys = ToMatrixSystem[eqn, bcs, {G}, {x, 0, 1}, μ] /. βc ->
        0.921 /. λc -> 1.775*10^-4;

Now we can evaluate the Evans function at a given value of $\mu$ and $V$ with e.g. Evans[μ, sys /. V -> 1]. Roots of this function correspond to eigenvalues of that third order equation.

Plot[{Evans[μ, sys /. V -> 1], Evans[μ, sys /. V -> 2]}, {μ, 2, 3}]

enter image description here

You can see there is a periodic set of eigenvalues coming out for positive values of $\mu$, and it is worth noting that the eigenvalues for $V=2$ are exactly twice those for $V=1$.

FindRoot[Evans[μ, sys /. V -> 1], {μ, 2}]
{μ -> 1.97219}
  • 1
    $\begingroup$ Thanks a lot for this. So I get the mu eigen-values for the $G$ equation using your method. Can these be not utilized to determine the other function with non-homogeneous condition ? I am just drawing on analogy here from standard PDE problem examples. $\endgroup$
    – Avrana
    May 14, 2020 at 20:26
  • $\begingroup$ Or probably I should try to convert this problem into a homogeneous system using some form of transformation (something like $\bar{{\theta_h}}(x,y):=\theta_h(x,y)-1$). If I could do that, then both F and G could be solved simultaneously using your method. Subsequently, what should be the next step in that scenario ? $\endgroup$
    – Avrana
    May 15, 2020 at 2:41
  • $\begingroup$ I have found (in your code) that in your definition of sys, you have used eqns instead of eqn which defines the $G$ equation. When I plot the same system with same parameters the eigenvalue plots are different. I have added the plot as an addendum to my original question. $\endgroup$
    – Avrana
    May 15, 2020 at 4:10
  • 1
    $\begingroup$ Ah, sorry, that was my mistake, I had a lingering definition and didn't check anything before posting last night. I'm not a PDE expert, so I'm not sure what steps you'll need to make next, but yes, you can solve the two eigenvalue BVPs simultaneously if they are homogeneous. $\endgroup$
    – SPPearce
    May 15, 2020 at 5:46
  • 1
    $\begingroup$ Thanks for the bounty award, even though it isn't a complete answer, I hope that my method is at least helpful for finding the eigenvalues. $\endgroup$
    – SPPearce
    May 21, 2020 at 8:24

Your Answer

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