# Mathematica - using NIntegrate within NDSolve

• Mathematica
• maves
In summary, the problem is that the program is trying to solve a system of differential equations, but it can't do it analytically, so it uses new variables x1, x2, y1, y2, z1, z2, to try and solve them 1st order equations. The solutions it's looking for should be in particular planes, but when it tries to plot the solutions, it always gets a linear line. It seems that it's missing one minus sign in the expression for the second integral, which makes the trajectories for particles in the gravitational potential "hyperbolic".

#### maves

Hello,
I am computing orbits of constant energy in a gravitational potential of a complicated planet (sth like a torus).

I need to solve a system of 2nd order differential equations, that look something like
$$\ddot x=-\frac{\partial \phi}{\partial x}$$ and the same for y and z. The problem is that ∅ is a complicated integral that can not be solved analitically. I transformed the system of equations to a system of 1st order equations using new variables x1, x2, y1, y2, z1, z2 (is that even necessary?) and tried to solve them using NDSolve. The solutions I am looking for should be lying in particular planes, so I put for instance z=0 and then defined the integrals (already partially derived):

int1[x1_?NumericQ, y1_?NumericQ] :=
NIntegrate[((R (x1 - R Cos[\[Phi]]))/((x1 - R Cos[\[Phi]])^2 + (y1 - R Sin[\[Phi]])^2 + (0 -
Z)^2))^(3/2), {\[Phi], 0, 2 Pi}, {R, b, b + c}, {Z, R-(b+c), -R + (b + c)}];
int2[x1_?NumericQ, y1_?NumericQ] :=
NIntegrate[((R (y1 - R Sin[\[Phi]]))/((x1 - R Cos[\[Phi]])^2 + (y1 - R Sin[\[Phi]])^2 + (0 -
Z)^2))^(3/2), {\[Phi], 0, 2 Pi}, {R, b, b + c}, {Z, R-(b+c), -R + (b + c)}];

and then used these in NDSolve:
s = NDSolve[{
x1'[t] == x2[t],
y1'[t] == y2[t],
x2'[t] == int1[x1[t], y1[t]],
y2'[t] == int2[x1[t], y1[t]],
x1[0] == 9.6, x2[0] == 0, y1[0] == 0, y2[0] == 1.52}, {x1, x2, y1,
y2}, {t, 0, 100}][[1]],

having chosen some random initial values (well, not exactly random, I believe x1-the initial position in direction x - should be close to the planet (9.5), y1 is not needed and y2 - the initial velocity in direction y - should be "something") .

I plotted the solutions {x(t),y(t)}:
ParametricPlot[Evaluate[{x[t], y[t]} /. s], {t, 0, 100},
PlotRange -> All, AxesOrigin -> {0, 0},
AxesLabel -> {"x1[t]", "y1[t]"}]

which is the trajectory of a particle in my gravitational potential ∅. The initial instruction tells me to find the orbits of lowest (constant) energy (so, as I understand, by changing the initial conditions/values), but when I plot my solutions, I always get a linear line, now matter what i put into the initial values. I don't know, whether the problem is in my understanding of the problem, the part written in Mathematica or somewhere else. I would be very glad if anyone could help me out, and I apologize for the long post, I hope I made it clear, what the problem is.

So I still have hope, that maybe someone could help me :) What I figured out until now is, that i forgot to take into account one minus sign, which made everything wrong, as the gravitational force in my calculations was repulsive and so the trajectory was hyperbolic (looked like linear on first sight). I also left the expression in the second-order form.

So, the right expressions for integrals in xy plane are:
1) Firstly, I define the first integral, which is equal to $-\frac{\partial \phi}{\partial x}$: (This "?NumericQ" must be here, if I want the program to proceed without evaluating the integral for a specific value of x and y. I integrate over $\Phi, R, Z$ and leave x and y "free", so that I will be able to evaluate the integrals at the points that NDSolve will calculate.

int1[x_?NumericQ, y_?NumericQ] :=
NIntegrate[((R (x - R Cos[\[Phi]]))/((x - R Cos[\[Phi]])^2 + (y - R Sin[\[Phi]])^2 + (0 -
Z)^2))^(3/2), {\[Phi], 0, 2 Pi}, {R, b, b + c}, {Z, R-(b+c), -R + (b + c)}];

2) Similarly, I define the second integral, which is equal to $-\frac{\partial \phi}{\partial y}$:

int2[x_?NumericQ, y_?NumericQ] :=
NIntegrate[((R (y - R Sin[\[Phi]]))/((x - R Cos[\[Phi]])^2 + (y - R Sin[\[Phi]])^2 + (0 -
Z)^2))^(3/2), {\[Phi], 0, 2 Pi}, {R, b, b + c}, {Z, R-(b+c), -R + (b + c)}];

Now I can write my set of equations and set the initial values, for instance:

s = NDSolve[{
x''[t] == int1[x[t], y[t]],
y''[t] == int2[x[t], y[t]],
x[0] == 9.6, x'[0] == 0, y[0] == 0, y'[0] == 1.5}, {x, x', y,
y'}, {t, 0, 10}][[1]]

And plot my solutions:
ParametricPlot[Evaluate[{x[t], y[t]} /. s], {t, 0, 10}]

x[0] has been chosen close to the toroid (9.5), because the potential is lowest near the toroid. The orbits in xy plane should be circular and go "around" the torus (like a satellite).
For comparison, I also plotted the torus in the xy plane, as the orbits should be circular or at least elliptical, and they seem completely wrong - like if there was no torus at all, they cross the planet, which is of course wrong - something like the case shown in attachment, remember, this is the xy plane!

I would be very glad if anyone could help me here - is there anything wrong with my code? but what could it be, as the code is not exactly rocket science -; or at least tell me, if you think my thread was posted to the wrong section.

#### Attachments

• tor21.jpg
9.6 KB · Views: 630
Last edited:

## 1. How do I use NIntegrate within NDSolve in Mathematica?

In order to use NIntegrate within NDSolve, you need to specify the integration variable and limits in the NIntegrate function and then include it in the differential equation as an integral term. For example: NDSolve[{y'[x] == x + NIntegrate[f[t], {t, 0, x}], y[0] == 0}, y, {x, 0, 1}]

## 2. Can I use NIntegrate for multiple integrals within NDSolve?

Yes, you can use NIntegrate for multiple integrals by specifying the variables and limits for each integral in the NIntegrate function and including it in the differential equation. For example: NDSolve[{y'[x] == x + NIntegrate[f[t,u], {t, 0, x}, {u, 0, 1}], y[0] == 0}, y, {x, 0, 1}]

## 3. Can I use NIntegrate for improper integrals within NDSolve?

Yes, you can use NIntegrate for improper integrals by specifying the limits as Infinity or -Infinity in the NIntegrate function. For example: NDSolve[{y'[x] == x + NIntegrate[f[t], {t, 0, Infinity}], y[0] == 0}, y, {x, 0, 1}]

## 4. Is it possible to use NIntegrate with symbolic constants within NDSolve?

Yes, you can use NIntegrate with symbolic constants by defining them beforehand and then including them in the NIntegrate function. For example: NDSolve[{y'[x] == x + NIntegrate[f[t], {t, 0, a}], y[0] == 0}, y, {x, 0, 1}] where a is a symbolic constant.

## 5. How do I specify options for NIntegrate within NDSolve?

You can specify options for NIntegrate by including them in the NIntegrate function. For example: NDSolve[{y'[x] == x + NIntegrate[f[t], {t, 0, x}, Method -> "AdaptiveMonteCarlo"], y[0] == 0}, y, {x, 0, 1}] This will use the AdaptiveMonteCarlo method for integration.