# How to solve the following system with NDSolve and FindRoot? Please address

Consider the following coupled systems of equations

$f′′− k^2 f=k\, \mathbf{Ra}$ and $g′′− k^2 g=k\,f$ , where $f$ and $g$ are functions of $z$. The boundary conditions for the problem are $f=0$ and $g=0$ at $z=\{0,1\}$. Here we have two parameters $\mathbf{Ra}$ (Rayleigh number) and $k$ (wave number).

Now, I want to find the minimum of $\mathbf{Ra}$ as $k$ varies. How do we find this? Also, I want to find and plot the eigenfunctions $f$ and $g$ in the interval. I would like to do this using Mathematica.

• Welcome to Mathematica StackExchange. Please take the time to read the tour. You will probably get better answers if you can show that you've made an effort to solve your problem. What Mathematica code have you already tried? Where did it go wrong? Commented Aug 14, 2017 at 8:22
• Don't we need an additional relationship between $\mathbf{Ra}$ and $k$ to learn anything about $\mathbf{Ra}$'s dependence on $k$? I believe your differential equations hold for any functional form of $\mathbf{Ra}(k)$. Commented Aug 14, 2017 at 9:20
• Yes sir, Ra is a function of k, Ra=Ra(k). Then I have two alternates. 1. Put outer loop k, for each k find Ra and then plot Ra vs k. This gives min of Ra as k varies. 2. Can we directly find critical k at which Ra is minimum. I think yes because at min or max of Ra, dRa/dk=0. we introduce this equation to the above system. Then I dont know how to make modifications to my code and plotting f and g. Would you mind helping me in this coding.
– PAL
Commented Sep 21, 2017 at 6:32

Update: Slowing down a little bit on relationship between $\mathbf{Ra}$ and $k$.

I believe we need more information from the original poster (perhaps standard fluid physics, but not mentioned in the post) to establish a dependent relationship between $\mathbf{Ra}$ and $k$. Any nonsingular function $\mathbf{Ra}(k)$ for $k>0$ could and will still satisfy the differential equations as written.

So I don't think I can help with finding the minimum of $\mathbf{Ra}$ with the information provided.

## Plotting $f(z), g(z)$

I think it's easiest to approach the problem (as currently stated) analytically.

eqns = {f''[z] - k^2 f[z] == k Ra,
g''[z] - k^2 g[z] == k f[z],
f[0] == g[0] == 0,
f[1] == g[1] == 0};

solns = Flatten[DSolve[eqns, {f[z], g[z]}, {z, 0, 1}]] // FullSimplify


yields

 {f[z] -> (Ra (-1 + Cosh[k z] - Sinh[k z] Tanh[k/2]))/k,
g[z] -> (1/(2 (1 + E^k)^2 k^2))  E^(-k z) Ra (2 E^(k z) +
2 E^(k (2 + z)) + 4 E^(k + k z) +
E^(k + 2 k z) (-2 + k (-1 + z)) + E^k (-2 + k - k z) +
E^(2 k z) (-2 + k z) - E^(2 k) (2 + k z))}


To get numbers out we need to choose $k$ and $\mathbf{Ra}$. In some units perhaps $\mathbf{Ra}=1$ and $k=1$. I'll plot with this choice.

Plot[Evaluate[{f[z], g[z]} /. solns /. Ra ->1 /. k -> 1], {z, 0,
1}, Frame -> True,
FrameLabel -> {z, TraditionalForm[Row[{f[z], ", ", g[z]}]]},
DisplayForm[SubscriptBox[f[z], k -> 1]]],


yielding

## Solving numerically

To solve numerically (using NDSolve) we need to choose values for the parameters.
As in the analytic plot above will set $\mathbf{Ra}=1$ and $k=1$ in some units.

solnNumeric =
NDSolve[eqns /. k -> 1 /. Ra -> 1, {f[z], g[z]}, {z, 0, 1}] //
Flatten


and

Plot[Evaluate[{f[z], g[z]} /. solnNumeric], {z, 0, 1}, Frame -> True,
FrameLabel -> {z, TraditionalForm[Row[{f[z], ", ", g[z]}]]},