# Numerically solving a system of nonlinear ODEs with boundary conditions

I have a system of 6 second-order nonlinear ODEs involving 5 different functions of a variable $t$. Every function has a boundary condition at $0$. I've never taken a differential equations class and know very little about them, but unfortunately they just popped up in my research and now I really need to solve this system. I think it's impossible or very difficult to solve analytically and a numerical approximation would be good enough for me. I was hoping to do it in Matlab because I'm very familiar with Matlab. If Matlab already has a function for this that would be amazing but I couldn't find anything. My supervisor suggested the shooting method but upon looking this up (and Runge-Kutta methods in general), it appears to be only for one equation and not a system. Am I missing something? If not, can anybody recommend an algorithm and somewhere I can read about it for problem? Speed and ease-of-programming is more important than accuracy (although within 5%ish of the true solution would be nice) because I have to solve the system many times for different-valued constants. I can list the equations if that would help. Thank you!

• See MATLAB ode45 and other ode** functions as a place to begin with. These handle systems of equations as simply as single ODEs. – Daryl Jul 11 '13 at 1:05
• @Daryl I had seen that, but I didn't know you could turn a system of equations into a single equation. How exactly does that work? Also, reading this article it seems like they need all the equations to be of the form "Y' = ...", but my equations are all implicit. One last thing - isn't this an exact solver? I actually tried putting it through dsolve, but that made Matlab take two hours then crash without a solution. – Samuel Reid Jul 11 '13 at 2:47

## 1 Answer

Since you have not given the set of differential equations, this answer cannot respond exactly. However, as indicated in the comments, the system of equations is implicit. These can be written as $F(t,y,y')=0$, where $$y(t)=[y_1(t),y_2(t),\ldots,y_n(t)]^T \text{ and } y'(t)=[y'_1(t),y'_2(t),\ldots,y'_n(t)]^T.$$This system can be solved numerically in MATLAB using ode$15$i.

As an example, suppose we want to solve the implicit system of ODEs $$F(t,y,y')=\begin{bmatrix}\cos(y'_1)+y'_2-y_1\\t\sin(y'_2)+y'_1-y_2^2\end{bmatrix}$$ with $y(0)=[0,0]^T$ for $0\leq t\leq T$.

ode$15$i requires consistent initial conditions for both $y$ and $y'$. Given that we do not have initial conditions, they can be approximated using the decic command.

This can be done in MATLAB as follows.

T = 5;
y0 = [0;0];
fun = @(t,y,yp) [cos(yp(1))+yp(2)-y(1);t*sin(yp(2))+yp(1)-y(2)^2];

[y0,yp0] = decic(fun,0,y0,[1;1],[0;0],[0;0]);%get consistent initial conditions

[tsol,ysol] = ode15i(fun,[0,T],y0,yp0); % solve the system

plot(tsol,ysol(:,1),'b',tsol,ysol(:,2),'r') % plot the solution


Edit: Solution to include equations

System

The system of DEs to be solved is as follows: \begin{align} \theta '' &= \frac{\theta' \cos \theta}{R} + \frac{R'\sin \theta}{R} + \frac{\nu\sin\theta}{\kappa R} + \frac{\eta\cos\theta}{\kappa R} \tag{1}\\ \nu' &= \frac{\kappa}{2} \left( \theta' + \frac{\sin\theta}{R}-c_0 \right)\left( \theta'-\frac{\sin\theta}{R}-c_0\right) + \lambda\tag{2}\\ \eta' & =0\tag{3}\\ R' &= \cos \theta \tag{4}\\ z' &=- \sin \theta \tag{5}\end{align} No initial conditions are given, though it is stated that all zero initial conditions may be appropriate.

Solution

Noting $(3)$, the solution of $\eta'=0$ must be $\eta(t)=c$ for some constant $c$. Since $\eta(0)=0$, we must have $c=0$, giving $\eta(t)=0$. As such, $(3)$ is not needed, since $\eta(t)$ can be obtained analytically. I will assume, however, throughout the remainder of the solution below, that $\eta(t)$ need not be a zero constant, but simply any constant, as given by the initial conditions.

Furthermore, $(1)$ can be reduced to first order by introducing the auxiliary variable $A=\theta'$. Then $(1)$ becomes $$A'=\frac{A \cos \theta}{R} + \frac{R'\sin \theta}{R} + \frac{\nu\sin\theta}{\kappa R} + \frac{\eta\cos\theta}{\kappa R} \tag{6}$$ with the additional equation $$\theta'=A\tag{7}.$$ The system then required to be solved is $(2),\,(4),\,(5),\,(6)\text{ and }(7)$.

In MATLAB, an m-file for the ode function should be created similar to as follows

function F = odefun(t,y,yp,parms)
% evaluate the ode functions with the given parameters
% y and yp are ordered as [theta,nu,R,z,A]
% parms is a structure, containing eta, kappa, lambda, c0

% extract components from parms structure
eta = parms.eta;
kappa = parms.kappa;
lambda = parms.lambda;
c0 = parms.c0;

% extract components from y and yp for more clear use
theta = y(1); thetap = yp(1);
nu = y(2); nup = yp(2);
R = y(3); Rp = yp(3);
z = y(4); zp = yp(4);
A = y(5); Ap = yp(5);

% allocate memory for return variable
F = zeros(5,1);

% precompute some quantities
sint = sin(theta);
cost = cos(theta);

% evaluate equation (2)
F(1) = nup - (kappa/2*(thetap+sint/R-c0)*(thetap-sint/R-c0)+lambda); %****

% evaluate equation (4)
F(2) = Rp - cost;

% evaluate equation (5)
F(3) = zp + sint;

% evaluate equation (6)
F(4) = Ap - (A*cost + Rp*sint + (nu*sint+eta*cost)/kappa)/R; %****

% evaluate equation (7)
F(5) = thetap - A;


Note: Both equations marked by %**** $\left(\text{i.e.$(2)$and$(6)$}\right)$ will have a divide by zero if $R(0)=0$.

From my previous smaller example, you should be able to make this function work with the ode15i(...) and decic(...) commands.

• Could you possibly do a second-order example? This was great but matlab doesn't like my attempt at generalizing. – Samuel Reid Jul 11 '13 at 19:56
• I added the actual equations to my question. Any help would be very much appreciated :) – Samuel Reid Jul 11 '13 at 20:59
• @SamuelReid With only a brief peruse, it looks like your only second order derivative is on $\theta$. This is usually overcome by introducing an auxiliary variable, say $\Theta$, such that $\Theta=\theta'$. Then $\theta''=\Theta'$, and the system is reduced to first order. I will expand my answer later today to include this, but it is too difficult on a mobile device. – Daryl Jul 11 '13 at 23:03
• @SamuelReid I also want to check your third equation. Is it actually $\eta'=0$. If so, then it can be simply removed, since $\eta=constant$. – Daryl Jul 11 '13 at 23:06
• Your edit is awesome and so are you. Thank you! – Samuel Reid Jul 24 '13 at 17:09