# How to solve this boundary value problem (2 coupled odes) with NDSolve

I am solving a boundary value problem which is consisted of two odes.

Mathematica prompts that SolveDelayed -> True should be considered. However, the problem can not be solved even if this option was used.

I have read this post and this post, but none of them helps.

Thanks for any suggestions.

Below is my code.

g = 15/10000;
Ri = 38/1000;
Ro = Ri + g;
m = 1/1000;
ηpos = 1/10;
τy = 2;
ϵ = 0.1;
λ = 1;
dpdz = 10^3;

η = τy/(Sqrt[u'[r]^2] + 10^-10) (1 - Exp[-m Sqrt[u'[r]^2]]) + ηpos;
τzz = τrz[r]^2 (2 λ)/η;
ode = {τrz[r] == η/(1 + (ϵ λ/η) τzz) u'[r], 1/r D[r τrz[r], r]== dpdz,
u[Ri] == 0, u[Ro] == 0};
sol = NDSolve[ode, {u, τrz}, {r, Ri, Ro}]

• I have removed your commented code, since you should endeavor to only include code relevant to the problem at hand (and since it's commented, it's likely irrelevant). Feel free to roll back those changes if you feel like they're important, though. Dec 11, 2019 at 4:57
• Thanks, it looks really better now. :) Dec 11, 2019 at 5:01
• Why do you say this is a simple boundary value problem? Dec 11, 2019 at 5:48
• Thanks for your attention. By 'simple', I mean that the physical problem behind the two odes is really common, and only two odes are loosely coupled. Sorry for this misleading word and I think I'd better to delete the word 'simple' in the title. Dec 11, 2019 at 10:42

Hint.

Solve it as an initial value problem, and then use a shooting method. The following formulation

g = 15/10000;
Ri = 38/1000;
Ro = Ri + g;
m = 1/1000;
ηpos = 1/10;
τy = 2;
ϵ = 1/10;
λ = 1;
dpdz = 10^3;

η = τy/(u'[r] + 10^-10) (1 - Exp[-m u'[r]]) +ηpos;
τzz = τrz[r]^2 (2 λ)/η;
ode = {τrz[r] == η/(1 + (ϵ λ/η) τzz) u'[r], 1/r D[r τrz[r], r] == dpdz,
u[Ri] == 0, τrz[Ri] == -0.75};
sol = NDSolve[ode, {u, τrz}, {r, Ri, Ro}, Method->{"EquationSimplification" -> "Residual"}];
Plot[Evaluate[u[r] /. sol], {r, Ri, Ro}, PlotStyle -> {Thick, Blue}]


Practically solves the boundary issue. Use then the shooting method.

• Thanks for the valuable hint and it took me some time to understand what you mean, and finally I solved the problem according to your great hint. This is the first time that I realized the relationship between ic and bc problems :) Dec 11, 2019 at 10:47

You can regard the trz separetly.

First I rationalize all parameters and rewrite equation for simple testing later.

ode = {0 == -τrz[
r] + η/(1 + (ϵ λ/η) τzz) u'[r],
1/r D[r τrz[r], r] == dpdz, u[Ri] == 0, u[Ro] == 0}


Since trz does not depend on u, treat it separatly.

dsol = DSolve[τrz[r]/r + Derivative[1][τrz][r] ==
1000, τrz, r]

{*  {{τrz -> Function[{r}, 500 r + C[1]/r]}    *}


Insert the result and apply a starting condition instead of boundary condition, since NDSolve can not find an explicit equation for u'[r] and therefore switches to differential-algbraic equation.

ode2 = Join[Most[ode] /. First@dsol /. C[1] -> c, {u'[Ri] == aa}]


Omitt the trz part

ode3 = ode2[[{1, 3, 4}]]


Find a solution for the c parameter in dependance of aa with a little trick, since Solve does not find a solution for undefined aa.

cc2[aa_] =
c /. First@
Solve[ode3[[1]] /. r -> Ri /. u'[Ri] -> -Pi, c, Reals] /.
Pi -> -aa;


Since shooting method did not work with Version 8.0, I guessed a value for aa. Find the exact aa with FindRoot.

ndsol = NDSolve[ode3 /. c -> cc2[aa] /. aa -> -88.1, u, {r, Ri, Ro}]


The solution satisfies ode.

u[Ri] /. ndsol
u[Ro] /. ndsol
Plot[u[r] /. ndsol, {r, Ri, Ro},
Epilog -> {Red, Point[{Ro, u[Ro] /. First@ndsol}]}]


and the error (would be lower with aa found by FindRoot)

Plot[Evaluate[
ode3[[1, 2]] /. c -> cc2[aa] /. aa -> -88.1 /. First@ndsol], {r, Ri,
Ro}, PlotRange -> All]


• I really like your semi-analytical way to attack the problem, and it is very useful to my problem (which is larger than the mini problem in this post). I am sorry, I cannot check both answers :) Dec 11, 2019 at 10:54
• Thanks. No problem. Dec 11, 2019 at 14:53