# Solving a system of differential and integral equations in Mathematica

Here is a system of differential equation I am trying to solve:

$$u''(x) - a(x) u(x) + b(x) v(x) = 0 \quad (1)$$ $$\gamma v''(x) + a(x) u(x) - b(x) v(x) = 0 \quad (2)$$

where $$a(x)$$, $$b(x)$$ are known functions in $$x$$ and $$\gamma$$ is a constant. The following are the boundary conditions:

$$u'(x=0) = u'(x=1) = v'(x=0) = v'(x=1) =0 \quad (3)$$

In order to avoid trivial solutions there is mass conservation constraint as follows:

$$\int_{0}^{1} (u(x)+v(x))dx = 1 \quad (4)$$

One can apply the constraint of equation 4 using finite differences, by applying the conservation law instead of expanding one of the equations for a particular cell.

1. Is there a built-in way to include this constraint in NDSolve instead of writing a code for the discretization myself?

2. Can finite element or shooting method be used to solve this problem in mathematica in case I wanted to avoid finite differences?

I tried this first:

eqs = {u''[x] - 15*(1 - x)*u[x] + 5*v[x] == NeumannValue[0, x == 0] + NeumannValue[0, x == 1], 1.2* v''[x] + 15*(1 - x)*u[x] - 5*v[x] == NeumannValue[0, x == 0] + NeumannValue[0, x == 1], Integrate[(u[x] + v[x]), {x, 0, 1}] == 1};
sol = NDSolve[eqs, {u, v}, {x, 0, 1}]

NDSolve::overdet: There are fewer dependent variables, {u[x],v[x]}, than equations, so the system is overdetermined.

• We can introduce new variable w'[x]== u[x] + v[x] with DirichletCondition[w[x] == 1, x == 1], DirichletCondition[w[x] == 0 x == 0], but problem has no solution, so we have message NDSolve::fempsf: PDESolve could not find a solution. Nov 27, 2022 at 18:30
• All your equations are linear: the differential equations as well as the integral constraint. So you can just solve without regards to the integral equation, and then scale the solution so that the mass constraint is satisfied. Nov 27, 2022 at 18:32
• @AlexTrounev How can a first order differential equation have two Dirichlet BCs? Nov 28, 2022 at 1:40
• @Roman can you please elaborate? If I try to solve the Eq. (1) and (2) only, NDSolve returns a trivial solution. Nov 28, 2022 at 1:42
• @Brownian_Motion Actually this is coupled system of ODEs, so we can use several appropriate boundary conditions, for example, use DirichletCondition[w[x] == 1,x==1] and DirichletCondition[w[x] == 0, x == 0], but remove one NeumannValue. Nov 28, 2022 at 4:33

eqs = {u''[x] - 15*(1 - x)*u[x] + 5*v[x] == 0,
1.2*v''[x] + 15*(1 - x)*u[x] - 5*v[x] == 0, u'[0] == 0, u'[1] == 0,
v'[0] == 0, v'[1] == 0};


If you have a solution {u,v} you may multiply it and the new function still fulfill the equations. Therefore, we may get a solution without regard to the integral and subsequently scale the solution to give the correct integral:

sol = {u[x], v[x]} /. NDSolve[eqs, {u, v}, {x, 0, 1}][[1]]
Plot[sol, {x, 0, 1}]


Now the integral:

int = NIntegrate[sol[[1]] + sol[[2]], {x, 0, 1}]
(* 2.4553 *)


has not the correct value. We therefore scale the solution:

sol1 = sol/int;
NIntegrate[sol1[[1]] + sol1[[2]], {x, 0, 1}]
(* 1. *)

Plot[sol, {x, 0, 1}]


• This makes a lot of sense. Why does not the code work with NeumannValues? I tried removing the integral and solve the equations with NeumannValue and it complains that since there are no Dirichlet BC therefore the solution might not be unique. Nov 28, 2022 at 14:12
• Neumann gives derivatives and Dirichlet values. Derivatives without values are not unique. Nov 28, 2022 at 15:24
• Thanks for the comment. I still am slightly confused why directly defining the derivative condition (e.g. u'[1] == 0 worked and NeumannValue[0,x==1] did not work. For my system, don't they mean the same thing? Nov 28, 2022 at 17:16
• NeumannValue[0,x==1] seems to imply the trivial solution. Nov 28, 2022 at 17:28

I don't claim this is optimal, but you can go one order further and then have an ODE system with a boundary condition in U==Integrate[u] and likewise for v. In order for this to work you need to add a boundary condition though, so I choose to have U[0]==0. Here is code for this (I kept the lower-case letters and we'll just remember to take a derivative at the end).

eqs = {u'''[x] - 15*(1 - x)*u'[x] + 5*v'[x] == 0,
1.2*v'''[x] + 15*(1 - x)*u'[x] - 5*v'[x] == 0, u''[0] == 0,
u''[1] == 0, v''[0] == 0, v''[1] == 0,
u[1] - u[0] + v[1] - v[0] == 1, u[0] == 0};
sol = NDSolveValue[eqs, {u[x], v[x]}, {x, 0, 1}];

Plot[Evaluate@D[sol, x], {x, 0, 1}]