Matlab Application : Differential Equations and

        Boundary Value Problems

Differential Equations:

Solution of differential equations is also referred to as simulation ( specially if the independent variable is time and the equations describe system dynamics.  For non linear systems very often this is the only way to obtain the solution.

Here we consider the Blasius solution to the laminar flow over a flat plate
 

The nonlinear ordinary differential equations(ODE) are described by the following two-point boundary value problem:  The prime indicates derivative wrt  the independent variable

 Boundary Value Problems

The previous problem we considered was actually a nonlinear two-point boundary value problem.   Matlab does not have a procedure to solve such a problem (not that I am aware of) .  We are going to develop our own.
 

Let us speculate how we can solve the problem.

f f ''      + 2 f ''' = 0

x = 0 ,     f = 0 ,     f' = 0
x = inf,     f' = 1

Note: the final value of the indepenedent variable is infinity.  Normally in computation we replace this with a finite value.  We do not however know this final value -  these problems are FREE FINAL VALUE problems.  In TPBVP this final vaue is solved along with any missing initial value.  To handle these problems (as well as efficient programming style) these problems have the final value of the independent variable normalized to 1.  To do this a transformation is involved:

xnew = x * tau ;        0 <= xnew  <= 1;
          tau is the actual final value (to be solved for in free time problem)

Also all the system equations must be multiplied by  'tau'

y1' = tau * y2 ,   etc
 

To use the Runge-Kutta method we need to know    f''(0) and we don't need  f'(infinity) = 1.   One way is to GUESS  a value of f''(0)  and use Runge-Kutta and discover what the value of f'(inf).   If it is 1 then we have a solution .  If not we have to make another guess.   We can use the information at the final point to change our guesses.   This can take a lot of guesses - as the equations are non-linear.   Such methods are natural in the subject of optimization.   Matlab provides one function ( not a efficient one) to handle minimization chores - if we do not have the optimization toolbox).   In boundary value problems these methods are called shooting methods
 

If we do not complete this in class it needs to be completed as part of home work

Two Point Boundary Value Problem: There will be four files used

    The Tpbvp.m script m-file    (initializes and starts the solution)
    The MinFun.m  function m-file    (sets up the optimization problem)
    The BC.m  function m-file    (sets up the Boundary conditions)
    the state.m function m-file      (holds the set of differential equations)

The Tpbvp file:

% Solution to TBBVP ( Two Point Boundary value Problem) for ODE
%
% Dr. P.Venkataraman, RIT Mechanical Engineering
%  June 1998
%
% Using ODE45 and FMINS (non-stiff differential equations)
% Problem must be described in state space as a
% System of first order equations in a function M-file - state.m
%-----------
%  STATE.m  |
%-----------
% [ STATE.m is similar to Function F(T,Y) described in ODE45]
% except the differential system must be normalized with respect
% to final time - i.e. the final time is one (1)-final time not known before
%
% function dsys = state(x, y, flag, tf)
% where x: independent variable
%       y: dependent vector - state derivatives
%       flag: ignore
%       tf:  final time
% ----   all derivatives must be multipled by tf ---
% This is useful for unknown final time problems
%
% Boundary Conditions are described in a function M-file  - BC.m
%--------
% BC.m  |
%-------
%  function bvector = BC(t0,tf, yi, yf)
%           Y(T0) = Vector of initial conditions of state
%           Y(Tf) = Vector of Final conditions of state after call to ODE45
% The function should return returns an vector of
%  Boundary Conditions set up with the
%  Right Hand Side to be 0 (ZERO) when satisfied
%
%-----------
% MinFun.m
%----------
% Calls the Runge-Kutta method  and sets up the optimization problem
% based on the least square error in the boundary conditions
% and final time
% function ReturnVal = MinFun(yd,NS,ftime)
%
 
 

% input the starting design vector - intial conditions fo the problem

global x y        % global makes the variables available
                % in the other functions
                % otherwise variables only have function scope
 
 

ok = 1;
while (ok)  % this is true
   string1 = ['\nInput a starting guess for the initial value vector' ...
         '\nThis is mandatory as the length of the vector indicates the ' ...
         '\norder of the system. Please enter it now and hit return : \n'];

   yinit = input(string1);
   if (isempty(yinit) == 0)
      break;
   end
end
NSYS = length(yinit);
fprintf('\nThe initial state vector(guessed) is : \n');
disp(yinit);
fprintf('\n')  %  line feed in the display

% check if final time is unknown
finaltime = input('Is final time unknown ? (y/n)','s');

if finaltime == 'y' | finaltime == 'Y' | finaltime == 'yes' ...
      | finaltime == 'YES'
   tfguess  = input('final time - if unknown then this is a guess [1]: ')
   ftm = tfguess;
   if isempty(tfguess)
      tfguess = 1;
      ftm = tfguess;
   end
   fprintf('\n');
   y0 = [yinit tfguess]
else
    ftm = input('final time :  ')
   y0 = yinit;
end

% call the fmins function with additional parameters
[YD, OPT] = fmins('MinFun',y0,[0 , 1.0e-08],[], NSYS,ftm);

fprintf('\nNumber of iterations : %i5 ',OPT(10)); % this ia s format specifier
fprintf('\n\n')
fprintf('The final values for the design vector: \n')
disp(YD)

fprintf('\nStrike any key to continue \n')
fprintf('The last element above is the final time \n')
pause

z = [ x y];
disp (z);
%  you can plot if you like the output
 

The MinFun File
 

function ReturnVal = MinFun(yd, NS, ftime)
% The function which will organize the optimal adjustment
%  of the boundary conditions and final time
global x y
n= length(yd);
nstate = n-1;

if n > NS
 for i = 1:1:n-1
    ys(i) = yd(i);
   end
else
   for i = 1:1:n;
      ys(i) = yd(i);
   end
end

if n > NS
   tf = yd(n);
else
   tf = ftime;
end

tspan = [0, 1];
 

[x y] = ode45('state', tspan,ys,[],tf);

nd = length(x);

for i = 1:1:n-1
   yf(i) = y(nd,i);
end

bvec = feval('Bc',0,tf,ys,yf);

% setting up the least square error
NBC = length(bvec);
sum = 0.0;
for i = 1:NBC
   sum = sum + bvec(i)*bvec(i);
end

ReturnVal = sum;

 

The BC.m File

function bvector = BC(t0, tf, yi, yf)
%
% This is where the boundary conditions are specified
% The conditions are specified such that the right hand sides are zero
%  This will allow mixed BC's

one = yi(1) ;
two = yi(2) ;
three = yf(2) -1.0;
 

bvector = [one two three];

 

The state.m  File

function dsys = state(x, y, flag, tf)
% the state space definition of the system
dsys =tf* [ y(2), y(3), ( -y(1)*y(3))]';