# Solving 4 differential equations using NDSolve

I am having a problem using NDSolve to solve the following system of 4 equations, in 4 dependent variables $$(r,t,\gamma,p_{r})$$ and one independent variable $$z$$. Here we have:

\begin{align} t &= t(z) \\ r &= r(z) \\ \gamma &= \gamma(r, t; z) \\ p_{r} &= p_{r}(r, t; z) \end{align}

We also have 3 known data sets called $$AA(r, t; z), BB(r, t; z)$$, and $$CC(r, t; z)$$ from experiment. The 4 equations are given below by total (not partial) derivatives as

\begin{align} \frac{\mathrm{d}t}{\mathrm{d}z} &= \frac{\gamma}{\sqrt{\gamma^{2}-p_{r}^{2}-1}} \\ \frac{\mathrm{d}r}{\mathrm{d}z} &= \frac{p_{r}}{\sqrt{\gamma^{2}-p_{r}^{2}-1}} \\ \frac{\mathrm{d}\gamma}{\mathrm{d}z} &= \frac{\mathrm{d}r}{\mathrm{d}z} A(r, t; z) + B(r, t; z) \\ \frac{\mathrm{d}p_r}{\mathrm{d}z} &= \frac{\mathrm{d}t}{\mathrm{d}z} A(r, t; z) - C(r, t; z) \end{align}

But when I try to solve this using NDSolve, different errors keep popping up. For example, see this case and error below:

NDSolve[{D[t[z], z] == γ[r[z], t[z], z]/
Sqrt[γ[r[z], t[z], z]^2 - pr[r[z], t[z], z]^2 - 1],
D[r[z], z] == pr[r[z], t[z], z]/
Sqrt[γ[r[z], t[z], z]^2 - pr[r[z], t[z], z]^2 - 1],
D[γ[r[z], t[z], z], z] ==
D[r[z], z] AA[r[z], t[z], z, 0.0] + BB[r[z], t[z], z, 0.0],
D[pr[r[z], t[z], z], z] ==
D[t[z], z] AA[r[z], t[z], z, 0.0] - CC[r[z], t[z], z, 0.0],
r[0] == 0, t[0] == 0, γ[r[0], t[0], 0] == 1.1174,
pr[r[0], t[0], 0] == 0}, {r[z], t[z],
pr[r[z], t[z], z], γ[r[z], t[z], z]}, {z, 0.0, 1.76}]

NDSolve::overdet: There are fewer dependent variables, {r[z],t[z]}, than equations, so the system is overdetermined.

Originally I tried to put Dt in each left-hand side on each line, but I read in the Mathematica help about NDSolve that one should use D not Dt. It is a bit confusing and errors keep appearing.

Not sure what is it that I am doing wrong? Any suggestions on how to solve this issue?

UPDATE: user "Nasser" has asked that I clarify the defintions of AA, BB, CC functions. Since I am not able to retrieve the definition for these functions, which I call from another code file that includes C++, please use these approximate/equivalent definitions instead if needed:

AA[r_, t_, z_, k_] :=
N[-BesselJ[1, 240 r] Cos[2 Pi 11 10^9 t] Sin[(z Pi)/1.76]];
BB[r_, t_, z_, k_] :=
N[BesselJ[0, 240 r] Sin[2 Pi 11 10^9 t] Sin[(z Pi)/1.76]];
CC[r_, t_, z_, k_] :=
N[1/370 BesselJ[1, 240 r] Cos[2 Pi 11 10^9] Sin[(z Pi)/1.76]];

Please note that k here is a fixed constant (usually 0), so the variables are really just (r, t; z) here.

• And what is AA, BB and CC there? your example does not include them. So can't really try the code fully like this. You also have r[z] but when you write pr[r, t, z] you just use r inside. Too many problems to try the code like this. – Nasser May 21 '20 at 8:39
• @Nasser AA, BB and CC are functions that I call from another code, and they return a single real number for each (t,r,z) value. So, I don't know how to share them here. Please feel free to use any other simple function as dummy to replace them. Thanks for pointing our that some points didn't have [z]. I have fixed them (see edit). However, the problem is still there. Can you please advise on how else you would write such problem? I am a beginner so any suggestions are welcome! – user135626 May 21 '20 at 8:54
• Please feel free to use any other simple function as dummy to replace them but I think it is you who should make your example complete. This is what MWE means. Others should not have to make extra code just to try your code, when they do not even know the context and what to add to replace these missing functions. – Nasser May 21 '20 at 8:59
• @Nasser I understand. I would love to provide these functions here, but I don't know how to provide such a function. I am not that skilled. I currently call them from another file which calculates them. That file calls another C++ code that is huge and I don't know how it works. I am trying, to the best of my ability, to ask a reasonable question here, even though it might not look very complete... – user135626 May 21 '20 at 9:04
• @Nasser I got an idea. Give me a few moments, and I will edit the post with functions that replace these ones. – user135626 May 21 '20 at 9:07

Try this:

AA[r_, t_, z_,
k_] := -BesselJ[1, 240 r] Cos[2 Pi 11 10^9 t] Sin[(z Pi)/1.76];
BB[r_, t_, z_, k_] :=
BesselJ[0, 240 r] Sin[2 Pi 11 10^9 t] Sin[(z Pi)/1.76];
CC[r_, t_, z_, k_] :=
1/370 BesselJ[1, 240 r] Cos[2 Pi 11 10^9] Sin[(z Pi)/1.76];

NDSolve[{
D[t[z], z] == γ[z]/Sqrt[γ[z]^2 - pr[z]^2 - 1],
D[r[z], z] == pr[z]/Sqrt[γ[z]^2 - pr[z]^2 - 1],
D[γ[z], z] ==
D[r[z], z] AA[r[z], t[z], z, 0] + BB[r[z], t[z], z, 0],
D[pr[z], z] ==
D[t[z], z] AA[r[z], t[z], z, 0] - CC[r[z], t[z], z, 0],
r[0] == 0, t[0] == 0, γ[0] == 1.1174, pr[0] == 0},
{r, t, pr, γ}, {z, 0.0, 1.76}]

The quantities you're integrating, PR[z] = pr[r[z], t[z], z], and G[z] = γ[r[z], t[z], z], are functions of a single variable z. If you're hoping to recover pr[r, t, z] and γ[r, t, z] as functions of three variables this way, you cannot. But you can obtain the composition pr[r[z], t[z], z] and γ[r[z], t[z], z] along the trajectory of {r[z], t[z], z}.

• Now this is just a brilliant answer! I really appreciate it. Not only does the code now work, but the explanation you gave is also quite clear. Many thanks! – user135626 May 21 '20 at 17:43
• @user135626 You're welcome. And thanks for the accept and appreciation. :) – Michael E2 May 21 '20 at 18:23
• May I ask you please about the intuition that led you to writing it in this format? I mean, what is the basic idea here about when to write a composite function out in z only versus when to keep in terms of the dependents like r[z] and t[z]? In dynamical equations such terms happen a lot and this tends to be quite confusing for me - any general advice on this? – user135626 May 21 '20 at 21:32
• @user135626 In terms of problem-solving, it was somewhat backwards. I know NDSolve cannot solve f'[g[z]] == ... implicitly for f[u] where u = g[z] (it can do some delay diff eqs); NDSolve wants the equation code in the form F'[z] == .... Then I thought about what you were trying to do. Your TeX equations (just the four DEs) are exactly what I coded and would have coded if that was all I had. The simple explanation is that in NDSolve the dependent variables like pr have to be functions of the independent variable z, not a composition. Not sure how to explain it better. HTH – Michael E2 May 21 '20 at 22:00
• Yes, it is very helpful, indeed. Thank you once again for being so clear. The game of tracing all dependent and independent variables in some Hamiltonian dynamics can be quite a challenge sometimes. So it is good to have some clarity here! :) – user135626 May 22 '20 at 1:55