# Solving a system of coupled ODEs

I am trying to solve these 9 coupled ODEs

(v1[t] - t) n1'[t] + n1[t] v1'[t] == 0,

3 sigma1 n1[t] n1'[t] + (v1[t] - t) v1'[t] +phi'[t] == 0,

(v2[t] - t) n2'[t] + n2[t] v2'[t] == 0,

3 Q sigma2 n2[t] n2'[t] + (v2[t] - t) v2'[t] + Q phi'[t] == 0,

phi'[t] - (1/ne[t]) ne'[t] == 0,

phi'[t] - (sigmas1/ns1[t]) ns1'[t] == 0,

(vs20 + vs2[t] - t) ns2'[t] + ns2[t] vs2'[t] == 0,

3 Qs2 sigmas2 ns2[t] ns2'[t] + (vs20 + vs2[t] - t) vs2'[t] + Qs2 phi'[t] == 0,

a n1'[t] + b n2'[t] - ne'[t] - c ns1'[t] + d ns2'[t] == 0,


My problem is that I have a precision error that I can't fix probably due to the initial values.

Please, is there a numerical method to evaluate the initial conditions for each variable before solving the ODE?

sol = NDSolve[{(v1[t] - t) n1'[t] + n1[t] n1'[t] == 0, 3 sigma1 n1[t] n1'[t]
+ (v1[t] - t) v1'[t] + phi'[t] == 0, (v2[t] - t) n2'[t] + n2[t] v2'[t] ==
0, 3 Q sigma2 n2[t] n2'[t] + (v2[t] - t) v2'[t] + Q phi'[t] == 0, phi'[t] -
(1/ne[t]) ne'[t] == 0, phi'[t] - (sigmas1/ns1[t]) ns1'[t] == 0, (vs20 +
vs2[t]- t) ns2'[t] + ns2[t] vs2'[t] == 0, 3 Qs2 sigmas2 ns2[t] ns2'[t] +
(vs20+vs2[t]- t) vs2'[t] + Qs2 phi'[t] == 0, a n1'[t] + b n2'[t] - ne'[t] - c
ns1'[t] + d ns2'[t] == 0, n1[0] == 3, v1[0] == 1, phi[0] == 0, n2[0] == 0.52,
v2[0] == 0.25,ne[0] == ns1[0] == ns2[0] == 0, vs2[0] == 0.1}, {n1, v1, phi,
n2, v2, ne, ns1, ns2, vs2}, {t, 0, 0.9}, WorkingPrecision -> 50, PrecisionGoal
-> 30]


NDSolve::precw: The precision of the differential equation (.......) is less than WorkingPrecision (50.).

the plots give

 Block[{vs20 = 6, Q = 1/8, Qs2 = 1, sigma1 = 0.25, sigma2 = 0.25,sigmas1 = 1,
sigmas2 = 1, a = 0.335, b = 0.564, c = 1.9, d =3},Plot[Evaluate[{n1[t],v1[t],
phi[t]} /. sol], {t, 0, 0.9},PlotRange -> {-1, 4}, PlotRangePadding -> 0,
PlotLegends -> {"n1", "v1", "phi"}]]


Please, someone for help

## 1 Answer

Try DSolve (without timerange):

DSolve[{(v1[t] - t) n1'[t] + n1[t] n1'[t] == 0,
3 sigma1 n1[t] n1'[t] + (v1[t] - t) v1'[t] + phi'[t] ==
0, (v2[t] - t) n2'[t] + n2[t] v2'[t] == 0,
3 Q sigma2 n2[t] n2'[t] + (v2[t] - t) v2'[t] + Q phi'[t] == 0,
phi'[t] - (1/ne[t]) ne'[t] == 0,
phi'[t] - (sigmas1/ns1[t]) ns1'[t] ==
0, (vs20 + vs2[t] - t) ns2'[t] + ns2[t] vs2'[t] == 0, 3 Qs2 sigmas2 ns2[t] ns2'[t] + (vs20 + vs2[t] - t) vs2'[t] +
Qs2 phi'[t] == 0, a  n1'[t] + b n2'[t] - ne'[t] - c ns1'[t] + d ns2'[t] == 0,
n1[0] == 3, v1[0] == 1, phi[0] == 0, n2[0] == 0.52, v2[0] == 0.25,
ne[0] == ns1[0] == ns2[0] == 0, vs2[0] == 0.1}
, {n1, v1, phi, n2,v2, ne, ns1, ns2, vs2}, t]

(*{{n1 -> Function[{t}, 3], n2 -> Function[{t}, 0.52],
ne -> Function[{t}, 0], ns1 -> Function[{t}, 0],
ns2 -> Function[{t}, 0], phi -> Function[{t}, 0],
v1 -> Function[{t}, 1], v2 -> Function[{t}, 0.25],
vs2 -> Function[{t}, 0.1]}}*)


alternatively Rationalize your equation before applying NDSolve!

• But I get the same plots @Ulrich Neumann Jan 18, 2020 at 16:21
• without error message and no precision issue Jan 18, 2020 at 16:25
• @gallagher precw` is merely a warning, it doesn't necessarily mean the result is wrong. Jan 19, 2020 at 6:34
• Now when I try to plot the output, I always get the plot above, i.e. constant values for all solutions. I made a mistak? Jan 19, 2020 at 10:17
• @Gallagher What should be the solution? What does this model describe? Jan 19, 2020 at 17:07