Change boundary conditions for an ODE system

Question

I'm trying to solve numerically a problem like the next one: \begin{cases} y'(t) = x(t)+y(t) \\ x'(t) = x(t)-y(t) \\ x(0) = x_0 \\ y(T) = y_T \end{cases} with $t\in[0,T]$

I want to use ode45 in MATLAB, but in order to use it, the boundary conditions must be initial boundary conditons, so I can't use ode45 because $y(T) = y_T$.

I defined $\tilde{y}(t)=y(T-t)$, so $\tilde{y}'(t)=-y'(T-t)$ and $y(t)=\tilde{y}(T-t)$, then I got the following system

\begin{cases} \tilde{y}'(t) = -x(T-t)+\tilde{y}(T-t) \\ x'(t) = x(t)-\tilde{y}(T-t) \\ x(0) = x_0 \\ \tilde{y}(0) = y_T \end{cases}

and now I have initial boundary conditions. Now, the problem is that I have to evaluate $x$ and $\tilde{y}$ at $T-t$ and $t$, but that doesn't seem to be an option when I use ode45.

Any idea? I would appreciate any help.

Answer 1

Instead of trying to directly use ode45, do it as a root-finding problem: finding $y(0)=y_0$ so that $y(T)=y_T$ when you solve the ODE.

Answer 2

What you are trying to do will never work. You are given conditions on the values of the solution at two different points, a boundary value problem. There is, trivially, no bijective map of the time domain that will map both points on the same one to get an initial value problem. "Two points map to one" is a contradiction to "bijective".

You will have to solve the ODE system with a free initial condition at the first point and select among all those possible solution the one(s) that also satisfy the second boundary condition. In Matlab you can realise this single-shooting method as (untested)

function defect = shoot(y0)
    [T,U] = ode45(@odefun, [0,T], [ x0 y0 ])
    return U(end,2)-yT
end

y0 = fzero(@defect,yT)

---

Using the full Matlab facilities, there is a more versatile multiple-shooting method already implemented in the Boundary-Value-Problem solvers bvp4c or bvp5c. Use them like demonstrated in the documentation.

• Code the Equation

  function dudt = bvpfun(t,u)
      x=u(1); y=u(2);
      dudt = [ x+y x-y ];
  end
  

• Code the Boundary Conditions

  function res = bcfun(u0,uT)
      res = [ u0(1)-x0  uT(2)-yT ];
  end
  

• Form an Initial Guess

  xmesh = linspace(0,T,5);
  solinit = bvpinit(xmesh, [x0 yT]);
  

• Solve the Equation

  sol1 = bvp4c(@bvpfun, @bcfun, solinit);
  plot(sol1.x,sol1.y(1,:),'-o',sol1.x,sol1.y(2,:),'-o')
  legend( 'x', 'y')