# [✓] Solve two ODEs coupled together using NDSolve?

GROUPS:
 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 solveI 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.
27 days ago
5 Replies
 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.RegardsNeil
27 days ago
 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.
 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.