[?] Solve two ODEs coupled together using NDSolve?

Questions
, ,
1

Hi All,

I would like to solve two ODEs coupled together using NDSolve, please see below for my code. Basically, I have two variables, y’[x] and g’[x] (not g’‘[x]). The reason why I formulate my ODE using g’‘[x] is that I have integral boundary condition on g’[x], so I use the trick from here to reformulate my equation. I will be happy to provide an original form of the equation if there is a question that is related to this part.

My problem is Mathematica complaints about “Infinite expression encountered” for my second equation(g’'[x]) part. My guess is that in the denominator of this equation, It has y[x] term which equals to 0 at left boundary. Even though I tried to avoid this by calculating y[-1+eps]=eps, It appears to still have this issue.

ClearAll["Global`*"];
w = 0.7;
q = 0.5*w*(1 - x^2);
kappa = 1.469;
sol = With[{eps = 10^-5}, NDSolve[{
y'[x] == (560*y[x]^0.5*q - 64*y[x]^4 + 6*w^(7/8)*(72*kappa - 77)*x*q*(-w*y[x]))/(3*w^(7/8)*(96*kappa - 77)*q^2),
g''[x] == (2695*y[x]^0.5*q - (864*kappa - 385)*y[x]^4 - 693*w^(7/8)*kappa*y[x]*q*(-w*x))/(9*(96*kappa - 77)*y[x]^3),
y[-1 + eps] == eps, 
g[1 - eps] - g[-1 + eps] == 0, 
g'[0] == 0}, {y[x], g'[x]},
{x, -1 + 0.001, 1 - 0.001}]];

Methods that I have tried is to give initial condition to Mathematica like below. (As suggested by here). However, I encountered error " Initial conditions should be specified at a single point."

ClearAll["Global`*"];
w = 0.7;
q = 0.5*w*(1 - x^2);
kappa = 1.469;
sol = With[{eps = 10^-5}, NDSolve[{
     y'[x] == (560*y[x]^0.5*q - 64*y[x]^4 + 6*w^(7/8)*(72*kappa - 77)*x*q*(-w*y[x]))/(3*w^(7/8)*(96*kappa - 77)*q^2),
     g''[x] == (2695*y[x]^0.5*q - (864*kappa - 385)*y[x]^4 - 693*w^(7/8)*kappa*y[x]*q*(-w*x))/(9*(96*kappa - 77)*y[x]^3),
     y[-1 + eps] == eps, g[1 - eps] - g[-1 + eps] == 0, 
     g'[0] == 0}, {y[x], g'[x]},
    {x, -1 + 0.001, 1 - 0.001}, 
    Method -> {"Shooting", 
      "StartingInitialConditions" -> {y[-1 + eps] == eps, 
        g'[eps] == 0}}]];

Any help will be greatly appreciated! Also, I have a related question: Right now, I have 2 ODEs, but I am planning to add transient terms d/dt for each variable I am solving at here. Let’s assume Mathematica could solve these ODEs. Is it possible to solve the transient PDEs? It will be a transient 1D problem, and from the NDSolve documentation, it seems that Mathematica should have capabilities to solve it.

------------------------Update for the original equation-----------------

My x range is from -1 to 1. Here are equations that I would like to solve

CodeCogsEqn.png
CodeCogsEqn.png389×333 5.31 KB

I also noticed here that this problem can be formulated as an optimization problem. However, I am having difficulty in this line of code on that link:

sol2[bc2_, {xmin_, xmax_}] := 
 NDSolveValue[{y''[x] - y[x] == (x^2), y'[0] == bc2, y[0] == 1}, y, {x, xmin, xmax}];
int[bc2_?NumericQ] := NIntegrate[sol2[bc2, {0, 2}][x], {x, 0, 2}];
**y2 = sol2[NMinimize[(int[bc2V] - 5)^2, bc2V][[-1, -1, -1]], {-3, 3}]**
Plot[{y1[x], y2[x]}, {x, -3, 3}]

I don’t know what does [-1,-1,-1] has to do with the original problem or formulation.

2

Kai,

The problem is with the constraint g[1 - eps] - g[-1 + eps] == 0. I am not sure that NDSolve (or any numerical solver) can handle a constraint relating both ends of the numerical integration. Numerical algorithms start at one end and integrate up. Furthermore, y[x] is independent of g[x] so you can integrate that separately:

sol = With[{eps = 10^-5}, 
   NDSolve[{y'[
       x] == (560*y[x]^0.5*q - 64*y[x]^4 + 
         6*w^(7/8)*(72*kappa - 77)*x*q*(-w*y[x]))/(3*
         w^(7/8)*(96*kappa - 77)*q^2), 
     y[-1 + eps] == eps}, {y[x]}, {x, -1 + eps, 1 - eps}]];

followed by

sol2 = With[{eps = 10^-5}, 
   NDSolve[{g''[
        x] == (2695*y[x]^0.5*q - (864*kappa - 385)*y[x]^4 - 
          693*w^(7/8)*kappa*y[x]*q*(-w*x))/(9*(96*kappa - 77)*
          y[x]^3) /. sol[[1, 1]], g[1 - eps] == 0, 
     g'[0] == 0}, {g'[x]}, {x, -1 + eps, 1 - eps}]];

but I had to change the constraint on the ends. I suggest you look at your physics of the problem and see if you can reformulate that one constraint.

Regards

Neil

3

You might be able to work around that issue with the boundary conditions by changing variables. I suggest you post the original equations. One approach might be to assign a new variable h[x] to the difference of g so that h[0] is your current boundary condition on g.

4

Thank you, Neil. I use your method, and the integration from 0 to 1 gives me 0.0868, which is “close” to 0. I have a feeling that depends on where I start with my initial integration, I should have an optimization problem to weakly enforce this integral constraint. I have updated my OP and includes my equations and link related to this optimization approach. I really appreciated your help! Thanks!

5

As to your second question about the line

y2 = sol2[NMinimize[(int[bc2V] - 5)^2, bc2V][[-1, -1, -1]], {-3, 3}]

the [[-1,-1,-1]] is taking the last part of an expression 3 times. NMinimize returns a list in the form

 {functionMin,{variable->value}}

they wanted to extract only the value and nothing else so they used Part ([[ ]]) to get the last part [[-1]] 3 times.

This approach numerically solves the integral and then numerically solves for the boundary condition that gives a certain value for the integral of the function over a particular range.

6

Thank you for your thorough explanation. I think I have implemented the optimization approach, and I am able to solve this set of the equations right now.

Best,
Kai