Message Boards Message Boards

Solve ODE system with BVP?

GROUPS:

I have a system of 3 non-linear ODEs of 2nd order with 3 unknown functions and 6 boundary conditions at Pi. I need to find functions numerically, get their plots and then calculate an integral containing these 3 functions. First is the system: I tried solving it with NDSolve but it gives me an error (NDSolve::ndsz: At x == 1.9672446444378016`, step size is effectively zero; singularity or stiff system suspected.) and the plots are wrong. I solved this system in Maple earlier so I think I'm making some syntax mistake since I'm not very familiar with Wolfram... Here's my code:

a = -10^2; b = 1; c = -2*10^(-3); m = 10^(-2); Rpi = 10^(-3); mupi = 4; Ypi = 10^(-2);
system = {R0''[x] + R0'[x]*Cot[x] == Exp[2*mu[x]]/2*((R0[x])^2 - c/a - m^2*(Y0[x])^2/a), R0[x]*Exp[2*mu[x]] == 2*(1 - mu''[x] - mu'[x]*Cot[x]), Y0''[x] + Y0'[x]*Cot[x] - m^2*Exp[2*mu[x]]*Y0[x] == 0, R0[Pi] == Rpi, R0'[Pi] == 0, mu[Pi] == mupi, mu'[Pi] == 0, Y0[Pi] == Ypi, Y0'[Pi] == 0}
sol = NDSolve[system, {R0, mu, Y0}, {x, 0, Pi}]
Plot[Evaluate[R0[x] /. sol], {x, 0, Pi}, PlotRange -> All]

What am I doing wrong? I know that a cot(x) can create a singularuty, but only at 0 or Pi, not at 1.97... On Maple plots there are singularities at 0 and that's fine, but not in R(x):

enter image description here

POSTED BY: Aaeru Michi
Answer
25 days ago

Can you show the original system of equations?

POSTED BY: Alexander Trounev
Answer
25 days ago

It's x instead of theta, doesn't really matter...

enter image description here

POSTED BY: Aaeru Michi
Answer
25 days ago

I changed the sign of parameter $a$ and added options in NDSolve. Apparently, the parameters should be selected more carefully to exit the singular point.

a = 10^2; b = 1; c = -2*10^(-4); m = 10^(-2); Rpi = 
10^(-3); mupi = 4; Ypi = 10^(-2);
system = {R0''[x] + R0'[x]*Cot[x] == 
Exp[2*mu[x]]/2*((R0[x])^2 - c/a - m^2*(Y0[x])^2/a), 
R0[x]*Exp[2*mu[x]] == 2*(1 - mu''[x] - mu'[x]*Cot[x]), 
Y0''[x] + Y0'[x]*Cot[x] - m^2*Exp[2*mu[x]]*Y0[x] == 0, 
R0[Pi] == Rpi, R0'[Pi] == 0, mu[Pi] == mupi, mu'[Pi] == 0, 
Y0[Pi] == Ypi, Y0'[Pi] == 0};
sol = NDSolveValue[system, {R0, mu, Y0}, {x, 0, Pi}, 
Method -> {"StiffnessSwitching", "NonstiffTest" -> False}] // 
Quiet;
Plot[sol[[1]][x], {x, 0., Pi}, PlotRange -> All, 
AxesLabel -> {\[Theta], R}] // Quiet

fig1

POSTED BY: Alexander Trounev
Answer
24 days ago

Thanks a lot! I realized a mistake in signs in equations and got a correct answer thanks to you. enter image description here

Now I need to integrate some expression numerically, I hope you can help. Or should I create another discussion? Basically I need to integrate this:

Integrate[Sin[x]*Exp[2*mu[x]]*(a*(R0[x])^2+b*R0[x]+c-(Y0'[x])^2/2/Exp[2*mu[x]]-m^2/2*(Y0[x])^2), {x, 0, Pi}]

And in this case I don't know a correct answer because Maple doesn't seem to work with this kind of integral correctly. I'm not sure how to "get" solutions from sol to put them into the integral, like this maybe:

R1[x] = Evaluate[R0[x] /. sol]
mu1[x] = Evaluate[mu[x] /. sol]
Y1[x] = Evaluate[Y0[x] /. sol]
Y2[x] = Evaluate[Y0'[x] /. sol]
Integrate[Sin[x]*Exp[2*mu1[x]]*(a*(R1[x])^2+b*R1[x]+c-(Y2[x])^2/2/Exp[2*mu1[x]]-m^2/2*(Y1[x])^2), {x, 0, Pi}]

But it doesn't really help: Wolfram gives me errors and doesn't calculate it...

POSTED BY: Aaeru Michi
Answer
24 days ago

You solve numerically then better is use NIntegrate than Integrate. The correct syntax is:

R1[x] = Evaluate[sol[[1]][x]];
mu1[x] = Evaluate[sol[[2]][x]];
Y1[x] = Evaluate[sol[[3]][x]];
Y2[x] = Evaluate[D[sol[[3]][x], x]];
NIntegrate[Sin[x]*Exp[2*mu1[x]]*(a*(R1[x])^2 + b*R1[x] + c - (Y2[x])^2/2/Exp[2*mu1[x]] - m^2/2*(Y1[x])^2), {x, 0, Pi},WorkingPrecision -> 20]
(* 9.0569457332837941551 *)
POSTED BY: Mariusz Iwaniuk
Answer
24 days ago

Thank you, you totally helped me!

POSTED BY: Aaeru Michi
Answer
23 days ago

Now it looks like in Maple

a = -10^2; b = 1; c = -2*10^(-4); m = 10^(-2); Rpi = 
 10^(-3); mupi = 4; Ypi = 10^(-2);
system = {R0''[x] + 
     R0'[x]*Cot[x] == -Exp[2*mu[x]]/
      2*((R0[x])^2 - c/a - m^2*(Y0[x])^2/a), 
   R0[x]*Exp[2*mu[x]] == 2*(1 - mu''[x] - mu'[x]*Cot[x]), 
   Y0''[x] + Y0'[x]*Cot[x] - m^2*Exp[2*mu[x]]*Y0[x] == 0, 
   R0[Pi] == Rpi, R0'[Pi] == 0, mu[Pi] == mupi, mu'[Pi] == 0, 
   Y0[Pi] == Ypi, Y0'[Pi] == 0};
sol = NDSolveValue[system, {R0, mu, Y0}, {x, 0, Pi}, 
    Method -> {"StiffnessSwitching", "NonstiffTest" -> False}] // 
   Quiet;
Plot[sol[[1]][x], {x, 0., Pi}, PlotRange -> All, 
  AxesLabel -> {\[Theta], R}, PlotRange -> All] // Quiet

fig

POSTED BY: Alexander Trounev
Answer
24 days ago

Now we can calculate the integral (Sorry Mariusz, we changed the system of equations and parameters):

In[57]:= NIntegrate[
 Sin[x]*Exp[
   2*sol[[2]][x]]*(a*(sol[[1]][x])^2 + b*sol[[1]][x] + 
    c - (sol[[3]]'[x])^2/2/Exp[2*sol[[2]][x]] - 
    m^2/2*(sol[[3]][x])^2), {x, 0, Pi}]

Out[57]= 0.128767
POSTED BY: Alexander Trounev
Answer
24 days ago

Thank you once more!

POSTED BY: Aaeru Michi
Answer
23 days ago

Group Abstract Group Abstract