[Math] A problem for Matlab solver

Question

I can numercally solve a system of two ODEs for x and p:

dy1 = -y1 + y2
dy2 = -1 + y2

with the following Matlab code:

%main.m
[t y]=rk4sys(20,0.05);
[t' y]

%rk4sys.m
function [t,y] = rk4sys(n,h)
C = [0;-1];
M = [-1,1;0,1];
y(1, :) = [0, 0];
f =@(t,y)(C+M*y')';
t(1) = 0;
  for j = 1 : n
    k1 = f(t(j), y(j,:));
    k2 = f(t(j)+h/2, y(j,:)+h*k1/2 );
    k3 = f(t(j)+h/2, y(j,:)+h*k2/2 );
    k4 = f(t(j)+h, y(j,:)+h*k3);
    y(j+1,:) = y(j,:)+h*(k1+2*(k2+k3)+k4)/6;
    t(j+1)=t(j)+h;
  end
end

The main program produces:

ans =

     0         0         0
0.0500   -0.0013   -0.0513
0.1000   -0.0050   -0.1052
0.1500   -0.0113   -0.1618
0.2000   -0.0201   -0.2214
0.2500   -0.0314   -0.2840
0.3000   -0.0453   -0.3499
0.3500   -0.0619   -0.4191
0.4000   -0.0811   -0.4918
0.4500   -0.1030   -0.5683
0.5000   -0.1276   -0.6487
0.5500   -0.1551   -0.7333
0.6000   -0.1855   -0.8221
0.6500   -0.2188   -0.9155
0.7000   -0.2552   -1.0138
0.7500   -0.2947   -1.1170
0.8000   -0.3374   -1.2255
0.8500   -0.3835   -1.3396
0.9000   -0.4331   -1.4596
0.9500   -0.4862   -1.5857
1.0000   -0.5431   -1.7183

Now, I have to change the element of (1,3) in the matrix above (that is 0(zero) for now) in order for the last number of the third column (that is -1.7183) to be equal to 0(zero).

In other words, the element (1,3), which is 0(zero) for now, is just my initial guess. I have to change this number to meet that y2 (at t = 21), that is -1.7183 for now, be 0(zero).

I can do this in Excel, but I have failed finding the way to do it in Matlab. Any help would be greatly appreciated.

Answer 1

You are looking for a shooting method.

In other words, you have the following problem description: $$\mathbf{y}'(t) = f(\mathbf{y},t),$$ $$y_1(0) = 0,$$ $$y_2(\tau) = 0.$$

This is a boundary value problem, and you need to transform it into an initial value problem.

To do so, instead of solving the equations above, you are really solving $$\mathbf{y}'(t) = f(\mathbf{y},t),$$ $$y_1(0) = 0,$$ $$y_2(0) = b.$$

for some estimate of $b$. Your problem is to find $b$ such that the original condition is met.

One technique is known as shooting. You make an initial guess for $b$, compute the error at the right boundary for $y_2(\tau)$, and adjust $b$ based on your error. This is essentially a root-finding problem, and now you can use root-finding techniques to find the "root" $b$.

Essentially, the solution to an ODE initial value problem can be written as $y(t;y_0)$, or the solution $y(t)$ parameterized by the initial conditions $y_0$. This is because, broadly speaking, solutions to ODEs are not unique until you specify initial conditions. So the actual solution trajectory depends greatly on your choice of initial conditions $y_0$.

In your case, you have a second-order system, so you can compute $\mathbf{y}(t;\mathbf{y}_0)$, where $\mathbf{y}_0 = \begin{pmatrix} 0 \\ b\end{pmatrix}$. Since the trajectory that leads you to $y_2(\tau) = 0$ is unique, you can iterate over different values of $b$ to find the right result, and indeed, formulate the problem as solving $F(b) = 0$, where $F$ is your "function" that computes the solution to your ODE. In other words, treat the MATLAB code that you've written above as a function of the parameter $b$ that returns the value $y_2(\tau)$. Then, compute the root.

There are a number of ways to do this; the simplest method, although by far not the fastest, is through bisection. Newton's method may also work, but that is more difficult to set up. Broyden's method (the secant method) will also work.

EDIT: Below is a solution. If this is homework, you must attribute this code properly. This code is presented to show you how shooting works, not to solve the problem for you

The following code can be pasted in a single file called "odetest.m". This is not the highest-quality code; I banged it together to illustrate a point, and I tried to design it in such a way so as to be mathematical in how it was structured. Certainly it can be optimized. That's not the point.

First, we define an RK4 ODE solver. Then, we define a function $F(b)$ such that $F$ returns the value of $y_2(t=1)$. This value depends entirely on the parameter $b$; as such, we wish to find $b$ such that $F(b) = 0$.

We then use Broyden's method, which is a secant method, to compute the estimate for $b$ from an initial guess of both the value and the derivative. Finally, we call this, and then plot some things.

The resulting plot is here:

Broyden's method is not the only way to solve this problem; however, it converges quite quickly and as you can see, it is easy to implement.

function odetest
    [i,lambda,lh] = broyden(@F,0,3);

    figure;
    hold on
    colors = 'brgky';
    legend_string = {};
    for j=1:length(lh)
        F(lh(j),1,colors(mod(j-1,length(colors))+1));
        legend_string{j} = ['j  = ' num2str(j) ', b = ' num2str(lh(j))];
    end
    legend(gca,legend_string);

    plot([0 1],[0 0],'k'); % Plot the x-axis

    set(gcf,'PaperPositionMode','auto')
    print('-f1','-r360','-dpng','broydenshooting.png');

    answer = lambda
end

function x=F(b,ploton,color)
    % Extra arguments for plotting mode; when the function is called as a
    % function, just use a single argument
    ode = @(t,y)[-y(1)+y(2);-1+y(2)];
    y0 = [0; b];
    T = [0 1];
    h = 0.05;
    t = T(1):h:T(2);

    y = RK4(ode,.05,t,y0);
    if nargin > 1
        plot(t,y(2,:),color,'LineWidth',2)
    end
    x = y(2,end);
end

function y=RK4(f,h,t,y0)
    n = length(t);
    y = zeros(length(y0),n);
    y(:,1) = y0;
    for i=1:n-1
        k1 = f(t(i),y(:,i));
        k2 = f(t(i)+h/2, y(:,i)+h/2*k1);
        k3 = f(t(i)+h/2, y(:,i)+h/2*k2);
        k4 = f(t(i)+h, y(:,i)+h*k3);
        y(:,i+1) = y(:,i)+h/6*(k1+2*k2+2*k3+k4);
    end
end

function [i,lambda,lh] = broyden(func,IX,IJ)
% lambda = broyden(func,IX,IJ)
% IX -- initial guess for X
% IJ -- initial guess for J = dF/dX
%
% i -- number of steps needed to converge
% lambda -- solution
% lh -- lambda history, the computed lambda at each step
%
% Attempt to compute the solution to a system of nonlinear equations
% without explicitly requiring the computation of the derivatives of the
% function.
lh = [];
lambda = IX;
B = IJ;
    for i=1:30
        lh = [lh; lambda];
        f = func(lambda);
        s = B\(-f);
        x = s + lambda;
        y = func(x)-f;
        B = B + ((y-B*s)*s')/(s'*s);
        if norm(x-lambda)/norm(lambda) < 1e-10
            lambda = x;
            break
        end
        lambda = x;
    end
end