# A problem for Matlab solver

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.

-

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');

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

-
I am a young researcher in public policy. So, the word "homework" does not really apply to me. (It however reminds me of nice and young old days...) Thank you very much for your help. I will analyze your code and get back to you soon. Thanks again! –  Bill TP Aug 15 '12 at 18:27
Dear Gorcenski, I believe you are a big figure here in math forum. Would you please take a look at the following problem? math.stackexchange.com/questions/140265/…. Someone already gave the Maple code. But I wonder if you could solve the problem in Matlab, not in a symbolic way. It is an old problem in city economics. Old but difficult to me. I heard that 20,30 years ago people used to solve it in Fortran -- OMG! I discussed this problem with one of my friends in engineering, but it was not satisfactory. Sad enough. Thank you very much in advance. –  Bill TP Aug 15 '12 at 18:44
I am far from a big figure in any forum, I am simply someone trying to keep his mind fresh. I will take a look at that other question, but I cannot guarantee I will have the time to fully analyze it. –  Arkamis Aug 15 '12 at 21:27
I appreciate it. I think the problem should interest you. It looks easy, but it seems not actually. I've known this problem for two years, and am still trying solving it via Matlab. –  Bill TP Aug 15 '12 at 23:05