Matlab: Is it possible to numerically solve a system of ode's with a mixture of initial and terminal conditions?

Question

I'm trying to use ode45 to solve a system of ODE's:

[X,Y]=  ode45(@sys,[0, T],y0);

where,

function dy = sys(t,y)

        dy(1) = f_1(y)
        dy(2) = f_2(y)
        dy(3) = f_3(y)
end

The problem is that the function ode45 requires that y0 be initial values [y_1(0), y_2(0), y_3(0)], while in my system, I only have the values [y_2(0), y_3(0), y_3(T)] available.

Mathematically, this set of initial/terminal conditions should be enough to pin down the system, but is there any way I can work with that by ode45 or any other functions in MATLAB?

Thanks!

Answer 1

After digging in the Matlab documentation for a little bit, I think the more elegant way is to use the bvp4c function. bvp4c is a function specifically designed to handle boundary value problems like this, as opposed to ode**, which are really for initial value problems only. In fact, there's a whole set of other functions such as deval and bvpinit in Matlab that really facilitate the use of bvp4c. Here's the link to the Matlab documentation.

I'll post a brief (and perhaps a bit contrived) example here:

function [y1, y2, y3] = test(start,T)

solinit = bvpinit(linspace(0,3,10), [1,1,0]);
sol = bvp4c(@odefun,@bvpbc,solinit);

tspan = linspace(start,T,100);
S = deval(sol, tspan);
y1 = S(1,:);
y2 = S(2,:);
y3 = S(3,:);

plot (tspan,y1)

figure
plot (tspan,y2)

figure
plot (tspan,y3)


%% system definition & BVCs

function dydx = odefun(t,y)
    dydx(1) = y(1) + y(2) + t;

    dydx(2) = 2*y(1) + y(2);

    dydx(3) = 3 * y(1) - y(2);
end

function res = bvpbc(y0,yT)
   res= [y0(3) yT(2) yT(3)];
end

end

The test function outputs 3 vectors of solutions points for y1, y2 and y3 respectively, and also plots them.

Here are the variable paths plotted by Matlab:

Also, I found this video by professor Jake Blanchard from WMU to be very helpful.

Answer 2

One approach is to use the shooting method to iteratively solve for the unknown initial state y_1(0) such that the desired final state y_3(T) is satisfied.

The iteration proceeds by solving a nonlinear equation F = 0:

F(y_1(0)) = Y_3(T) - y_3(T)

where the uppercase function Y is the solution obtained by integrating the ODE's forward in time from a set of initial conditions. The task is to guess the value of y_1(0) to obtain F = 0.

Since we're now solving a nonlinear equation, all of the usual approaches apply. Specifically you could use either bisection or Newton's method to update your guess for the unknown initial state y_1(0). Note a couple of things:

1. The ODE's are integrated on [0,T] at each iteration of the nonlinear solve.
2. There may be multiple solutions for F = 0, depending on the structure of your ODE's.
3. Newton's method may converge faster than bisection, but may also be numerically unstable unless you can provide a good starting guess for y_1(0).

Using existing MATLAB functions, the bounded scalar nonlinear solver FMINBND might be a good choice as a nonlinear solver.

Hope this helps.