Numerically solving a system of nonlinear ODEs with boundary conditions

Question

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!

Answer 1

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.