Differential equations

I have 4 differential equations. Two 2nd order ODEs

eqns3 = θ1''[y] + Q1 == 0, θ1[0] == θh, -θ1'[0] == Qh
eqns4 = θ2''[y] + Q2 == 0, θ2[1] == 1, -θ2'[1] == Nc (θ2[1] - 1)


and two 4th order ODEs.

eqns1 = θs''''[y] - Bi (k + 1) θs''[y] - Bi k (wf + ws) == 0
eqns2 = θf''''[y] - Bi (k + 1) θf''[y] - Bi k (wf + ws) == 0


The boundary conditions for the 4th order ODEs

θ2[y2] == θf[y2],
θf[y1] == θ1[y1],
θ1'[y1] == ke1 θf'[y1] + k ke1 θs'[y1],
θ2'[y2] == ke2 θf'[y2] + k ke2 θs'[y2]

θs[y1] == θ1[y1],
θs[y2] == θ2[y2],
θ1'[y1] == ke1 θf'[y1] + k ke1 θs'[y1],
θ2'[y2] == ke2 θf'[y2] + k ke2 θs'[y2]


I guess, because the boundary conditions for the fourth order are coupled, Mathematica has been unable to solve for theta f and θs (they are still blue), which does not allow me to plot the graphs.

Also, I have the general solutions of these equations.

θs[y] = Es y^2 + Fs Cosh[y * Sqrt[Bi (k + 1)]] + K1s y + K2s
θf[y] = Ef y^2 + Ff Cosh[y * Sqrt[Bi (k + 1)]] + K1f y + K2f
θ1[y] = A1 y^2 + B1 y + C1
θ1[y] = A2 y^2 + B2 y + C2


Where Es, Fs, K1s, K2s, Ef, Ff, K1f, K2f, A1, B1, C1, A2, B2, C2 are all unknowns. Any tips on how I can solve this problem.

My code

eqns1 = θs''''[y] - Bi (k + 1) θs''[y] - Bi k (wf + ws) == 0
eqns2 = θf''''[y] - Bi (k + 1) θf''[y] - Bi k (wf + ws) == 0
eqns3 = θ1''[y] + Q1 == 0
eqns4 = θ2''[y] + Q2 == 0

DSolve[
{eqns2,
θ2[y2] == θf[y2],
θf[y1] == θ1[y1],
θ1'[y1] == ke1 θf'[y1] + k ke1 θs'[y1],
θ2'[y2] == ke2 θf'[y2] + k ke2 θs'[y2]},
θf[y], y]

DSolve[{eqns3, θ1[0] == θh, -θ1'[0] == Qh}, θ1[y], y]

DSolve[
{eqns1,
θs[y1] == θ1[y1],
θs[y2] == θ2[y2], θ1'[y1] == ke1 θf'[y1] + k ke1 θs'[y1],
θ2'[y2] == ke2 θf'[y2] + k ke2 θs'[y2]},
θs[y], y]

DSolve[{eqns4, θ2[1] == 1, -θ2'[1] == Nc (θ2[1] - 1)}, θ2[y], y]


Solving for θ1 and θ2 was okay. But for θs and θf it could not solve it. Mathematica gave me an answers for θf and θs but those answers had θf and θs in them.

• Can you supply the code in a workable form? I.e. what you actually tried yourself. – Feyre Aug 21 '16 at 16:02
• Yes I can. Is there a way I can attach my code file here? – Dammy Ige Aug 21 '16 at 16:22
• Why not just paste the code in your question? – Feyre Aug 21 '16 at 16:28
• Just pasted the code into the question – Dammy Ige Aug 21 '16 at 16:35
• You can use the {} button to format code. – Feyre Aug 21 '16 at 16:50

The systems of equations you have in your question are not sufficient to solve without two generated constants. When you look at the equations for θf and θs it looks like there are four pieces of information for each, eight total.

However since they are coupled and two of the equations are identical, there are really only six pieces of information. As a result you get two generated constants.

The work flow path I took was to solve the two 2nd order ODEs first.

DSolveValue[
{
θ1''[y] + Q1 == 0,
θ1'[0] == -Qh,
θ1[0] == θh
},
θ1[y],
y]

(* 1/2 (-2 Qh y - Q1 y^2 + 2 θh) *)


Define a function for it

θ1fun[y_] := 1/2 (2 θh - 2 Qh y - Q1 y^2)


and

DSolveValue[
{
θ2''[y] + Q2 == 0,
θ2[1] == 1,
θ2'[1] == Nc (1 - θ2[1])
},
θ2[y],
y]

(* 1/2 (2 - Q2 + 2 Q2 y - Q2 y^2) *)


and a function for it as well

θ2fun[y_] := 1/2 (2 - Q2 + 2 Q2 y - Q2 y^2)


The two coupled equations should be written in one DsolveValue. Use the functions from the solution of the two second order ODE.

{θfy, θsy} = DSolveValue[
{
θf''''[y] - Bi (k + 1) θf''[y] - Bi k (wf + ws) == 0,
θs''''[y] - Bi (k + 1) θs''[y] - Bi k (wf + ws) == 0,
θf[y2] == θ2fun[y2],
θf[y1] == θ1fun[y1],
θ1fun'[y1] == ke1 θf'[y1] + k ke1 θs'[y1],
θ2fun'[y2] == ke2 θf'[y2] + k ke2 θs'[y2],
θs[y1] == θ1fun[y1],
θs[y2] == θ2fun[y2]
},
{θf[y], θs[y]},
y];


This takes a long time. The output is so long that you can't look at it without pain.

You can check however that each has two generated constants.

Cases[θfy, C[i_Integer], Infinity]


The output here has 226 cases of C[5] or C[6].

The second function has fewer occurrences of generated constants.

Cases[θfyθsy[[2]], C[i_Integer], Infinity]

(* {C[5], C[5], C[5], C[5], C[5], C[5], C[6], C[6], C[6], C[6],
C[6], C[6]} *)


See if it is possible to re-formulate your problem so the ambiguity is removed. It might make the output simpler.

• FWIW, Variables[{θfy, θsy}] should yield all the degrees of freedom, including the constants of integration -- is that what you mean by "ambiguity"? Tho, you may have guessed right, that the OP won't expect or want them, once it's clear that they're there. – Michael E2 Aug 22 '16 at 9:47
• Thank you guys for your contribution. you are right that the boundary conditions are six instead of eight cos it is duplicated. This problem is from this paper sciencedirect.com/science/article/pii/S0735193316300537 . The writer used the the general solutions { θs[y] = Es y^2 + Fs Cosh[y * Sqrt[Bi (k + 1)]] + K1s y + K2s θf[y] = Ef y^2 + Ff Cosh[y * Sqrt[Bi (k + 1)]] + K1f y + K2f θ1[y] = A1 y^2 + B1 y + C1 θ1[y] = A2 y^2 + B2 y + C2 } and somehow arrived at a set of simultaneous equations and used Solve for the unknowns. But it is not clear what he actually did. – Dammy Ige Aug 22 '16 at 12:04
• @DammyIge What if the C1 and C2 are the constants of integration? (Probably not, but I don't have time to look into it right now.) – Michael E2 Aug 22 '16 at 16:51

I think the problem is that some of the equations are coupled, but the OP has written DSolve commands as if they are independent. Try this:

{θfθsSol} =
DSolve[DeleteDuplicates@{eqns2, θ2[y2] == θf[y2],
θf[y1] == θ1[y1], θ1'[y1] == ke1 θf'[y1] + k ke1 θs'[y1],
θ2'[y2] == ke2 θf'[y2] + k ke2 θs'[y2],
eqns1, θs[y1] == θ1[y1], θs[y2] == θ2[y2],
θ1'[y1] == ke1 θf'[y1] + k ke1 θs'[y1],
θ2'[y2] == ke2 θf'[y2] + k ke2 θs'[y2]},
{θf, θs}, y];


One can delete the duplicates by hand, too, if one wants to.

We can check that the solutions are now free of θf and θs with FreeQ:

FreeQ[{θf[y], θs[y]} /. θfθsSol, θf | θs]
(*  True  *)


Note that ReplaceAll (expr /. θfθsSol) is the standard way to substitute the solution into an expression, which one may observe throughout the documentation for DSolve, especially in the plotting examples.

• Oh well, maybe all the equations should be put into one big system, although the 3rd & 4th seem to have perfectly self-contained solutions. They could be plugged back into θfθsSol if desired. – Michael E2 Aug 22 '16 at 1:43
• the 4th order equations were gotten from two second order ODEs [(1/k) [Theta]f''[y] + Bi ([Theta]s - [Theta]f) + wf == 0 [Theta]s''[y] - Bi ([Theta]s - [Theta]f) + ws == 0] that were coupled. After uncoupling, it became a 4th order ODE. I ran the answer you pasted above and did not get any output. and the FreeQ said true. But what does that mean? – Dammy Ige Aug 22 '16 at 12:44
• @DammyIge Did you read the (linked) documentation page for FreeQ? True means the expression contains no instances of θf or θs (i.e. it's "free" of those symbols). My code gave no (visible) output because of the semicolon. The output is stored in the variable θfθsSol. It a set of two replacement rules that are very long. Execute θfθsSol in a cell by itself to see the output. – Michael E2 Aug 22 '16 at 12:49
• Now I understand. Executed θfθsSol and still have unknowns C[5] and C[6] in the whole expression just like @JackLaVigne and what I have been getting. Is there by any chance another method mathematically that I can explore? – Dammy Ige Aug 22 '16 at 13:30
• @DammyIge The reason you have two free constants of integration C[5] and C[6], as Jack points out, is your equations have more degrees of freedom than constraints (boundary conditions). If there are two boundary conditions that should be imposed (that you forgot about, say), then including them would get rid of them. Otherwise, they are free parameters. You can set each to any number and still get a solution. Maybe that makes no sense in your case, but it does in some cases. What to do has to come from the problem itself and the related science or math. In other words, it's up to you. – Michael E2 Aug 22 '16 at 16:50