1
$\begingroup$

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}]
Improve this question
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – march
    Dec 11, 2019 at 4:57
  • $\begingroup$ Thanks, it looks really better now. :) $\endgroup$
    – xinxin guo
    Dec 11, 2019 at 5:01
  • $\begingroup$ Why do you say this is a simple boundary value problem? $\endgroup$
    – xzczd
    Dec 11, 2019 at 5:48
  • $\begingroup$ 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. $\endgroup$
    – xinxin guo
    Dec 11, 2019 at 10:42

2 Answers 2

Reset to default
1
$\begingroup$

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.

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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 :) $\endgroup$
    – xinxin guo
    Dec 11, 2019 at 10:47
2
$\begingroup$

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}]}]

enter image description here

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]

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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 :) $\endgroup$
    – xinxin guo
    Dec 11, 2019 at 10:54
  • $\begingroup$ Thanks. No problem. $\endgroup$
    – Akku14
    Dec 11, 2019 at 14:53

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.