1
$\begingroup$

I am trying to solve the following three coupled ODEs with six boundary conditions:

deq = {pr[x] + a y'[x] + b pr'[x] - c pr''[x] == 0, 
e pt[x] + e2 pr'[x] + l pt'[x] + b2 pr[x]/x - dd pt''[x] == 0, 
v0 pr'[x] + v0 pr[x]/x - dd y'[x]/x - dd y''[x] == 0, y[L] == y0, 
pt[L] == pt0, pr[L] == pr0, y'[L] == v0 pr[L]/dd, 
pr'[L] == (pr[L] + v0 y[L])/(2 dd), 
pt'[L] == v0 (b2 pr[L] + e pt[L])/(2 dd)} /. 
Thread[Rule[{a, b, c, e, e2, l, b2, dd, L, pt0, pr0, v0}, ConstantArray[1, 12]]];

DSolve[deq, {pr, pt, y}, x]

but there is no meaningful output. Could anyone help me to solve the equations? What is wrong in my code?

I don't get solutions to the equations even when I put constants equal to 1:

L = 5;
deq = {pr[x] +  y'[x] + pr'[x] -  pr''[x] == 0, 
pt[x] + pr'[x] + l pt'[x] + pr[x]/x - pt''[x] == 0, 
pr'[x] + pr[x]/x - y'[x]/x - y''[x] == 0, y[L] == y0, 
pt[L] == pt0, pr[L] == pr0, y'[L] == pr[L], 
pr'[L] == (pr[L] + y[L]), pt'[L] == pr[L] + pt[L]} /. 
Thread[Rule[{a, b, c, e, e2, l, b2, dd, L, pt0, pr0, v0}, ConstantArray[1, 12]]];

DSolve[deq, {pr, pt, y}, x]

I tried the first answer as below:

L = 5;
DI = 1;
f = 1;
g = 1;
ff = 1;
fff = 1;
fo = 1;
gg = 1;
DD = 1;
Y0 = 1;
pt0 = Cos[Pi/6];
pr0 = Sin[Pi/6];
deq = {pr[x] + g y'[x] + gg pr'[x] - DD pr''[x] == 0, 
pt[x] + f pr'[x] + ff pt'[x] + fff pr[x]/x - fo pt''[x] == 0, 
v0 pr'[x] + v0 pr[x]/x - DI y'[x]/x - DI y''[x] == 0, y[L] == Y0, 
pt[L] == pt0, pr[L] == pr0, y'[L] == v0 pr[L]/DI, 
pr'[L] == (pr[L] + y[L]), pt'[L] == pr[L] + pt[L]} /. 
Thread[Rule[{a, b, c, e, e2, l, b2, dd, L, pt0, pr0, v0, y0}, 
 ConstantArray[1, 13]]];

sol = NDSolve[deq, {pr, pt, y}, {x, 10^-6, 4}];

Plot[Evaluate[{pr[x], pt[x], y[x]} /. sol], {x, 0, 5}, PlotLegends -> Placed[{"pr[x]", "pt[x]", "y[x]"}, {.75, .3}]]

but the output graph shows that pt and pr are not equal to boundary condition at x=l.

Improve this question
$\endgroup$
2
  • 3
    $\begingroup$ If a symbolic solver returns unevaluated, then Mathematica can't solve the problem. I believe anyone asking for help with such a case should demonstrate that it is reasonable to expect that the problem has a closed form solution. Mathematica is not a magic box that'll spit out a solution to any problem. $\endgroup$
    – Szabolcs
    Commented Jan 17, 2018 at 9:42
  • $\begingroup$ @Szabolcs So true... $\endgroup$
    – Henrik Schumacher
    Commented Jan 17, 2018 at 9:50

3 Answers 3

Reset to default
2
$\begingroup$

Solution by Maple 2017.3 (By magic box)

I chose these to solve:

deq = {pr[x] + a y'[x] + b pr'[x] - c pr''[x] == 0, 
e pt[x] + e2 pr'[x] + l pt'[x] + b2 pr[x]/x - dd pt''[x] == 0, 
v0 pr'[x] + v0 pr[x]/x - dd y'[x]/x - dd y''[x] == 0, y[L] == y0, 
pt[L] == pt0, pr[L] == pr0, y'[L] == v0 pr[L]/dd, 
pr'[L] == (pr[L] + v0 y[L])/(2 dd), 
pt'[L] == v0 (b2 pr[L] + e pt[L])/(2 dd)} /. 
Thread[Rule[{a, b, c, e, e2, l, b2, dd, L, pt0, pr0, v0}, ConstantArray[1,12]]];

enter image description here

Translated from Maple to MMA:

  {pr[x] = (-(1/6))*E^(1 - x)*(-3 + y0) + (1/6)*E^(-2 + 2*x)*(3 + y0), 
  pt[x] = -((E^(-(5/2) - Sqrt[5]/2 + (1/2)*(1 + Sqrt[5])*x)*(-5*E^2 - 
  5*Sqrt[5]*E^2))/(10*Sqrt[5])) + 
  (E^(-(5/2) + Sqrt[5]/2 - (1/2)*(-1 + Sqrt[5])*x)*(-5*E^2 + 
  5*Sqrt[5]*E^2))/(10*Sqrt[5]) + 
  (E^(Sqrt[5]*x - (1/2)*(-1 + Sqrt[5])*x)*Integrate[(-(1/6))*E*(-3 + y0)*(-
  E^(-((3*z1)/2) - (Sqrt[5]*z1)/2) + E^(-((3*z1)/2) - (Sqrt[5]*z1)/2)/z1) + 
  ((3 + y0)*(2*E^((3*z1)/2 - (Sqrt[5]*z1)/2) + E^((3*z1)/2 - 
  (Sqrt[5]*z1)/2)/z1))/(6*E^2), {z1, 1, x}])/Sqrt[5] - 
  Integrate[(-(1/6))*E*(-3 + y0)*(-E^(-((3*z1)/2) + (Sqrt[5]*z1)/2) + E^(-
  ((3*z1)/2) + (Sqrt[5]*z1)/2)/z1) + 
  ((3 + y0)*(2*E^((3*z1)/2 + (Sqrt[5]*z1)/2) + E^((3*z1)/2 + 
  (Sqrt[5]*z1)/2)/z1))/(6*E^2), {z1, 1, x}]/(E^((1/2)*(-1 + 
  Sqrt[5])*x)*Sqrt[5]), 
  y[x] == y0 + Integrate[-(((1/6)*E^(1 - z1)*(-3 + y0)*z1 - (1/6)*E^(-2 + 
  2*z1)*(3 + y0)*z1)/z1), {z1, 1, x}]}

And Integrated for x>1:

 solA = {pr[x] = -(1/6) E^(1 - x) (-3 + y0) + 1/6 E^(-2 + 2 x) (3 + y0), 
 pt[x] = -((
 E^(-(5/2) - Sqrt[5]/2 + 
 1/2 (1 + Sqrt[5]) x) (-5 E^2 - 5 Sqrt[5] E^2))/(10 Sqrt[5])) + (
 E^(-(5/2) + Sqrt[5]/2 - 
 1/2 (-1 + Sqrt[5]) x) (-5 E^2 + 5 Sqrt[5] E^2))/(10 Sqrt[5]) + (
 1/(6 Sqrt[5]))
 E^(-2 + Sqrt[5] x - 
 1/2 (-1 + Sqrt[5]) x) ((
 2 E^(3 - 
 1/2 (3 + Sqrt[5]) x) (-1 + E^(
 1/2 (3 + Sqrt[5]) (-1 + x))) (-3 + y0))/(3 + Sqrt[5]) - (
 4 E^(-(Sqrt[5]/
 2)) (-E^(3/2) + E^(1/2 (Sqrt[5] - (-3 + Sqrt[5]) x))) (3 + 
 y0))/(-3 + Sqrt[
 5]) + (3 + y0) (-ExpIntegralEi[1/2 (3 - Sqrt[5])] + 
 ExpIntegralEi[-(1/2) (-3 + Sqrt[5]) x]) + 
 E^3 (-3 + y0) (ExpIntegralEi[1/2 (-3 - Sqrt[5])] - 
 ExpIntegralEi[-(1/2) (3 + Sqrt[5]) x])) - 
 1/(6 Sqrt[5])
 E^(-2 - 1/2 (-1 + Sqrt[5]) x) ((
 2 E^(3/2) (-E^((Sqrt[5]/2)) + E^(
 1/2 (3 + (-3 + Sqrt[5]) x))) (-3 + y0))/(-3 + Sqrt[5]) + (
 4 (-E^(1/2 (3 + Sqrt[5])) + E^(1/2 (3 + Sqrt[5]) x)) (3 + y0))/(
 3 + Sqrt[5]) + 
 E^3 (-3 + y0) (ExpIntegralEi[1/2 (-3 + Sqrt[5])] - 
 ExpIntegralEi[1/2 (-3 + Sqrt[5]) x]) + (3 + 
 y0) (-ExpIntegralEi[1/2 (3 + Sqrt[5])] + 
 ExpIntegralEi[1/2 (3 + Sqrt[5]) x])), 
 y[x] = y0 + 
 1/12 (3 + 2 E^(1 - x) (-3 + y0) - 3 y0 + E^(-2 + 2 x) (3 + y0))}

Numerical check:

(* for analitic solution: *)

  y0 = 1;
  solA /. x -> 4 // N
  (* {268.969, 503.393, 135.46} *)

(* for numeric solution:*)

  L = 1;
  sol = NDSolve[deq, {pr, pt, y}, {x, 1, 4}];
  {pr[x], pt[x], y[x]} /. sol /. x -> 4 // N
  (*{{268.969, 503.393, 135.46}}*)
Improve this answer
$\endgroup$
5
  • $\begingroup$ Thank you for your great answer :). Do you think below equation can also be solved by the same method? $\endgroup$
    – Wat Watson
    Commented Jan 17, 2018 at 13:13
  • $\begingroup$ EQ = {y[x] + x y'[x] == v0 pr[x]/Dd + x v0 pr'[x]/Dd, (4 v0 pr0 x + Dd) pr[x] + v0 x^2 y[x]/2 + (2 v0 x^2 pr0 - Dd x) pr'[x] - x^2 Dd pr''[x] == 0, y[L] == 1, pr[L] == pr0, y'[L] == v0 pr[L]/Dd, pr'[L] == v0 (pr0 pr[L] + y[L]/2)/Dd} $\endgroup$
    – Wat Watson
    Commented Jan 17, 2018 at 13:13
  • $\begingroup$ Something seems to be wrong… just try solrule = Unevaluated@(*The last code block*)/. HoldPattern@(a_[x] = b_) :> (a -> Function @@ {x, b}); Simplify[deq /. solrule, x > 1] The 1st equation isn't satisfied. $\endgroup$
    – xzczd ♦
    Commented Jan 17, 2018 at 13:21
  • $\begingroup$ @WatWatson.Maple says an Error:Error, (in dsolve) found differentiated functions with same name but depending on different arguments in the given DE system: {pr(L), pr(x), y(L), y(x)} $\endgroup$
    – Mariusz Iwaniuk
    Commented Jan 17, 2018 at 13:24
  • 1
    $\begingroup$ @xzczd.I edited my answer. Maple's simplify command something does not work properly or converting from MathML is weird. $\endgroup$
    – Mariusz Iwaniuk
    Commented Jan 17, 2018 at 13:53
2
$\begingroup$

Use NDSolve

L = 5;
deq = {pr[x] + y'[x] + pr'[x] - pr''[x] == 0, 
    pt[x] + pr'[x] + l pt'[x] + pr[x]/x - pt''[x] == 0, 
    pr'[x] + pr[x]/x - y'[x]/x - y''[x] == 0, y[L] == y0, pt[L] == pt0, 
    pr[L] == pr0 , y'[L] == pr[L], pr'[L] == (pr[L] + y[L]), 
    pt'[L] == pr[L] + pt[L] } /. 
   Thread[Rule[{a, b, c, e, e2, l, b2, dd, L, pt0, pr0, v0, y0}, 
     ConstantArray[1, 13]]];

sol = NDSolve[deq, {pr, pt, y}, {x, 10^-6, 4}];

LogPlot[Evaluate[{pr[x], pt[x], y[x]} /. sol], {x, 10^-6, 4}, 
 PlotLegends -> Placed[{"pr[x]", "pt[x]", "y[x]"}, {.75, .3}]]

enter image description here

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thank you for your answer. But I prefer to solve them without ndsolve to understand the effect of each parameter. Do you think it could be possible to do that or not? $\endgroup$
    – Wat Watson
    Commented Jan 17, 2018 at 6:11
  • $\begingroup$ Actually, I need to find a functionality for the solutions. $\endgroup$
    – Wat Watson
    Commented Jan 17, 2018 at 6:13
  • $\begingroup$ Since DSolve cannot solve the system with numeric values for all parameters, it is extremely unlikely that a solution can be found for symbolic parameters. $\endgroup$
    – Bob Hanlon
    Commented Jan 17, 2018 at 6:20
  • $\begingroup$ What do you mean by "numeric values for all parameters"? Can DSolve find the functionality if I determine some parameters and not all of them? $\endgroup$
    – Wat Watson
    Commented Jan 17, 2018 at 6:25
  • 1
    $\begingroup$ The Thread[Rule[{a, b, c, e, e2, l, b2, dd, L, pt0, pr0, v0, y0}, ConstantArray[1, 13]] sets all of the parameters to 1, i.e., assigns "numeric values for all parameters". Even with these assignments, DSolve cannot solve the simplified system. Complicating the problem by leaving any of the parameters as symbolic is not going to help DSolve. $\endgroup$
    – Bob Hanlon
    Commented Jan 17, 2018 at 6:33
1
$\begingroup$

If you try to find the general solution of your equations without any boundary conditions

DSolve[{pr[x] + Derivative[1][pr][x] +Derivative[1][y][x] - (pr^\[Prime]\[Prime])[x] == 0, 
pr[x]/x + pt[x] + Derivative[1][pr][x] +Derivative[1][pt][x] - (pt^\[Prime]\[Prime])[x] == 0,
pr[x]/x + Derivative[1][pr][x] - Derivative[1][y][x]/x - (y^\[Prime]\[Prime])[x] == 0 }
, {pr, pt, y}, x]

(*DSolve[{pr[x] + Derivative[1][pr][x] +Derivative[1][y][x] - (pr^\[Prime]\[Prime])[x] == 0, 
pr[x]/x + pt[x] + Derivative[1][pr][x] +Derivative[1][pt][x] - (pt^\[Prime]\[Prime])[x] == 0, 
pr[x]/x + Derivative[1][pr][x] - Derivative[1][y][x]/x - (y^\[Prime]\[Prime])[x] == 0}, {pr, pt, y}, x]*)

MMA can't solve it! So any further activity using DSolve for the complete system seems to be useless!!!

Let's try some simplication: The first and the last ode only depends on y[x] , pr[x] and can be separated:

Derivative[1][y][x] = -pr[x] -Derivative[1][pr][x] + (pr^\[Prime]\[Prime])[x] (* first ode *)

substituting in the third ode 

oder=2 pr[x] + (1 + 2 x) Derivative[1][pr][x] - (pr^\[Prime]\[Prime])[x] +x(pr^\Prime]\[Prime])[x] - x\!\(\*SuperscriptBox[\(pr\),TagBox[RowBox[{"(","3", ")"}],Derivative],MultilineFunction->None]\)[x]

gives ode in pr which can be solved 

ergpr = DSolve[oder == 0, pr, x][[1]]
(* {pr -> Function[{x},E^-x C[1] + E^(2 x) C[2] +3 E^-x C[3] (1/3 E^(3 x) ExpIntegralEi[-2 x] - ExpIntegralEi[x]/3)]} *)

With this result y'[x]is known and y[x] can be evaluated 

ergy= DSolve[y'[x] == -pr[x] -Derivative[1][pr][x] + (pr^\[Prime]\[Prime])[x] /.ergpr, y,x][[1]]
(*{y -> Function[{x}, -E^-x C[1] + 1/2 E^(2 x) C[2] + C[4] + 
1/2 E^(2 x) C[3] ExpIntegralEi[-2 x] + 
E^-x C[3] ExpIntegralEi[x] - 1/2 C[3] Log[-2 x] + 2 C[3] Log[x]]}*)

Unfortunately the remaining second ode

odept = -E^-x C[1] + 2 E^(2 x) C[2] +3 E^(2 x) C[3] ExpIntegralEi[-2 x] +(E^-x C[1] + E^(2 x) C[2] +3 E^-x C[3] (1/3 E^(3 x) ExpIntegralEi[-2 x] -ExpIntegralEi[x]/3))/x -3 E^-x C[3] (1/3 E^(3 x) ExpIntegralEi[-2 x] -ExpIntegralEi[x]/3) +pt[x] + Derivative[1][pt][x] - (pt^\[Prime]\[Prime])[x] == 0

can not be solved analytically by MMA!

Adaption of the boundary condition y, pr perhaps simplifies the last equation...

Improve this answer
$\endgroup$

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.