clc
clf
clear all

% Mfile name
%       mtl_int_sim_spehigherorder.m

% Revised:
%       March 7, 2008

% % Authors
%       Nathan Collier, Autar Kaw
%       Department of Mechanical Engineering
%       University of South Florida
%       Tampa, Fl 33620
%       Contact: kaw@eng.usf.edu | http://numericalmethods.eng.usf.edu/contact.html

% Purpose
%       To illustrate the use of the 4th order Runge-Kutta method to solve 
%       a higher order ODE as applied to a function of the user's choosing.

% Keyword
%       Runge-Kutta
%       Ordinary Differential Equations
%       Higher Order
%       Initial Value Problem

% Inputs
%       This is the only place in the program where the user makes the changes
%       based on their wishes

% Define the differential equation to be solved for of the form
% 
%      A1*d/dx(dy/dx) + A2*dy/dx + A3*y + A4 = 0
%
% The coefficients a1, a2, and a3 can be functions of x and y while
% a4 can be a function of x only. Define these functions here. 

   A1 = @(x,y) 1 ;
   A2 = @(x,y) 2 ;
   A3 = @(x,y) 1;
   A4 = @(x) -exp(-x) ;

% Create functions 1 and 2 for use in the method

   f1 = @(x,y,z) z ;
   f2 = @(x,y,z) 1/A1(x,y)*(-A2(x,y)*z-A3(x,y)*y-A4(x)) ;

% For the exact solution, the previous two functions need expressed in the 
% following for. Replace z and y with y(2) and y with y(1)

   f = @(x,y) [ y(2); 1/A1(x,y(1))*(-A2(x,y(1))*y(2)-A3(x,y(1))*y(1)-A4(x)) ] ;

% x0, x location of known initial condition

   x0=0 ;

% y0, corresponding value of y at x0

   y0=1 ;

% xf, x location at where you wish to see the solution to the ODE

   xf=0.75 ;

% dydx0, the derivative dy/dx at x = x0

   dydx0 = 2 ;

% n, number of steps to take

   n=3 ;
   
%**********************************************************************

disp(sprintf('\n\nSolving Higher Order, Ordinary Differential Equations'))
disp('University of South Florida')
disp('United States of America')
disp('kaw@eng.usf.edu')
disp(sprintf('\nNOTE: This worksheet demonstrates the use of Matlab to illustrate'))
disp('the 4th Order Runge-Kutta method to solve higher order ODE''s')

disp(sprintf('\n***************************Introduction*****************************'))
 
disp('The ordinary differential equations solved in Euler''s and 2nd and 4th')
disp('order Runge-Kutta method files were all of the first order. These methods')
disp('can be used to solve higher order ordinary differential equations by converting')
disp('them into a coupled system of first order equations. The following sheet')
disp('demonstrates how this is done.')

disp(sprintf('\n\n****************************Input Data******************************'))
disp('Below are the input parameters to begin the simulation. Although this method')
disp('can be used to solve ODE''s of any order, the ODE being approximated here is')
disp('of second order.')
disp(sprintf('\n     A1*d/dx(dy/dx) + A2*(dy/dx) + A3*y + A4 = 0\n'))
disp('This second order equation can be written in terms of 2 first order equations.')
disp('The user must be able to break up his ODE into these two functions.') ;
disp(sprintf('\n     f1 = dy/dx = z    '))
disp(sprintf('     f2 = dz/dx = 1/A1 * ( A2*z + A3*y + A4 )'))
disp(sprintf('     x0 = initial x'))
disp(sprintf('     y0 = initial y'))
disp(sprintf('     xf = final x'))
disp(sprintf('     n = number of steps to take')) 
format short g
disp(sprintf('\n-----------------------------------------------------------------'))
disp(sprintf('\n     x0 = %g',x0))
disp(sprintf('     y0 = %g',y0))
disp(sprintf('     dydx0 = %g',dydx0))
disp(sprintf('     xf = %g',xf))
disp(sprintf('     n = %g',n))
disp(sprintf('\n-----------------------------------------------------------------'))
disp(sprintf('For this simulation, the following parameters are constant.'))

% compute the spacing
h=(xf-x0)/n ;
disp(sprintf('\n     h = ( xf - x0 ) / n'))
disp(sprintf('       = ( %g - %g ) / %g',xf,x0,n))
disp(sprintf('       = %g',h))

xa(1)=x0 ;
ya(1)=y0 ;
za(1)=dydx0 ;

disp(sprintf('\n\n********************************Simulation**********************************'))

for i=1:n
    disp(sprintf('\nStep %g',i))
    disp(sprintf('-----------------------------------------------------------------'))

% Use the 4th Order Runge-Kutta Method to determine the next step estimate.
    xa(i+1)=xa(i)+h ;
    k1y = f1(xa(i),ya(i),za(i)) ;
    k1z = f2(xa(i),ya(i),za(i)) ;
    k2y = f1(xa(i)+0.5*h,ya(i)+0.5*k1y*h,za(i)+0.5*k1z*h) ;
    k2z = f2(xa(i)+0.5*h,ya(i)+0.5*k1y*h,za(i)+0.5*k1z*h) ;
    k3y = f1(xa(i)+0.5*h,ya(i)+0.5*k2y*h,za(i)+0.5*k2z*h) ;
    k3z = f2(xa(i)+0.5*h,ya(i)+0.5*k2y*h,za(i)+0.5*k2z*h) ;
    k4y = f1(xa(i)+h,ya(i)+k3y*h,za(i)+k3z*h) ;
    k4z = f2(xa(i)+h,ya(i)+k3y*h,za(i)+k3z*h) ;
    ya(i+1)=ya(i)+1.0/6.0*(k1y+2*k2y+2*k3y+k4y)*h ;
    za(i+1)=za(i)+1.0/6.0*(k1z+2*k2z+2*k3z+k4z)*h ;

    disp(sprintf('1) Find k1y, k1z, k2y, k2z, k3y, k3z, k4y, and k4z using the previous'))
    disp(sprintf('step information and the constant parameters.'))
    disp(sprintf('     k1y = f1( x(%g) , y(%g), z(%g) )',i-1,i-1,i-1))
    disp(sprintf('         = f1( %g , %g, %g ) ',xa(i),ya(i),za(i)))
    disp(sprintf('         = %g',k1y))
    disp(sprintf('     k1z = f2( x(%g) , y(%g), z(%g) )',i-1,i-1,i-1))
    disp(sprintf('         = f2( %g , %g, %g ) ',xa(i),ya(i),za(i)))
    disp(sprintf('         = %g',k1z))
    disp(sprintf('     k2y = f1( x(%g) + 0.5 * h , y(%g) + 0.5 * k1y * h , z(%g) + 0.5 * k1z * h )',i-1,i-1,i-1))
    disp(sprintf('         = f1( %g + 0.5 * %g , %g + 0.5 * %g * %g , %g + 0.5 * %g * %g ) ',xa(i),h,ya(i),k1y,h,za(i),k1z,h))
    disp(sprintf('         = f1( %g , %g , %g )',xa(i)+0.5*h,ya(i)+0.5*k1y*h,za(i)+0.5*k1z*h))
    disp(sprintf('         = %g',k2y))
    disp(sprintf('     k2z = f2( x(%g) + 0.5 * h , y(%g) + 0.5 * k1y * h , z(%g) + 0.5 * k1z * h )',i-1,i-1,i-1))
    disp(sprintf('         = f2( %g + 0.5 * %g , %g + 0.5 * %g * %g , %g + 0.5 * %g * %g ) ',xa(i),h,ya(i),k1y,h,za(i),k1z,h)) 
    disp(sprintf('         = f2( %g , %g , %g )',xa(i)+0.5*h,ya(i)+0.5*k1y*h,za(i)+0.5*k1z*h))
    disp(sprintf('         = %g',k2z))
    disp(sprintf('     k3y = f1( x(%g) + 0.5 * h , y(%g) + 0.5 * k2y * h , z(%g) + 0.5 * k2z * h )',i-1,i-1,i-1))
    disp(sprintf('         = f1( %g + 0.5 * %g , %g + 0.5 * %g * %g , %g + 0.5 * %g * %g ) ',xa(i),h,ya(i),k2y,h,za(i),k2z,h))
    disp(sprintf('         = f1( %g , %g , %g )',xa(i)+0.5*h,ya(i)+0.5*k2y*h,za(i)+0.5*k2z*h))
    disp(sprintf('         = %g',k3y))
    disp(sprintf('     k3z = f2( x(%g) + 0.5 * h , y(%g) + 0.5 * k2y * h , z(%g) + 0.5 * k2z * h )',i-1,i-1,i-1))
    disp(sprintf('         = f2( %g + 0.5 * %g , %g + 0.5 * %g * %g , %g + 0.5 * %g * %g ) ',xa(i),h,ya(i),k2y,h,za(i),k2z,h))
    disp(sprintf('         = f2( %g , %g , %g )',xa(i)+0.5*h,ya(i)+0.5*k2y*h,za(i)+0.5*k2z*h))
    disp(sprintf('         = %g',k3z))
    disp(sprintf('     k4y = f1( x(%g) + h , y(%g) + k3y * h , z(%g) + k3z * h )',i-1,i-1,i-1))
    disp(sprintf('         = f1( %g + %g , %g + %g * %g , %g + %g * %g ) ',xa(i),h,ya(i),k3y,h,za(i),k3y,h))
    disp(sprintf('         = f1( %g , %g , %g ) ',xa(i)+h,ya(i)+h*k3y,za(i)+h*k3z))
    disp(sprintf('         = %g',k4y))
    disp(sprintf('     k4z = f2( x(%g) + h , y(%g) + k3y * h , z(%g) + k3z * h )',i-1,i-1,i-1)) 
    disp(sprintf('         = f2( %g + %g , %g + %g * %g , %g + %g * %g ) ',xa(i),h,ya(i),k3z,h,za(i),k3z,h))
    disp(sprintf('         = f2( %g , %g , %g ) ',xa(i)+h,ya(i)+h*k3y,za(i)+h*k3z))
    disp(sprintf('         = %g',k4z))
    disp(sprintf('2) Apply the Runge-Kutta 4th Order method to estimate y(%g) and z(%g)',i,i))
    disp(sprintf('     y(%g) = y(%g) + 1/6*( k1y + 2 * k2y + 2 * k3y + k4y ) * h ',i,i-1))
    disp(sprintf('           = %g + %g * %g ',ya(i),1/6*(k1y+2*k2y+3*k3y+k4y),h))
    disp(sprintf('           = %g',ya(i+1)))
    disp(sprintf('     z(%g) = z(%g) + 1/6*( k1z + 2 * k2z + 2 * k3z + k4z ) * h ',i,i-1))
    disp(sprintf('           = %g + %g * %g ',ya(i),1/6*(k1z+2*k2z+3*k3z+k4z),h))
    disp(sprintf('           = %g',za(i+1)))
    disp(sprintf('   at x%g = %g',i,xa(i+1)))

end


disp(sprintf('\n\n********************************Results**********************************'))

% The following finds what is called the 'Exact' solution. This technique
% is numerical but we are using it to avoid using dsolve, as dsolve will not
% solve odes that do not have exact solutions. 
xspan = [x0 xf];
[xe,ye]=ode45(f,xspan,[y0;dydx0]);
[yfi dummy]=size(ye);
yf=ye(yfi);

hold on
xlabel('x');ylabel('y');
title('Exact and Approximate Solution of the ODE by the 4th Order Runge-Kutta Method');
plot(xe,ye(:,1),'-','LineWidth',2,'Color',[0 0 1]);            
plot(xa,ya,'-','LineWidth',2,'Color',[0 1 0]);
legend('Exact','Approximation');

disp(sprintf('The figure window that now appears shows these %g steps as a series of',n))
disp(sprintf('green lines. The blue line represents the exact solution. In this case'))
disp(sprintf(' ''exact'' refers to the solution obtained by the Matlab function ode45.'))
disp(sprintf('\nNote the approximate value obtained as well as the true value and'))
disp(sprintf('relative true error at our desired point x = xf.'))

disp(sprintf('\n   Approximate = %g',ya(n+1))) 
disp(sprintf('   Exact       = %g',yf))
disp(sprintf('\n   True Error = Exact - Approximate')) 
disp(sprintf('              = %g - %g',yf,ya(n+1)))
disp(sprintf('              = %g',yf-ya(n+1)))
disp(sprintf('\n   Absolute Relative True Error Percentage'))
disp(sprintf('              = | ( Exact - Approximate ) / Exact | * 100'))
disp(sprintf('              = | %g / %g | * 100',yf-ya(n+1),yf))
disp(sprintf('              = %g',abs( (yf-ya(n+1))/yf )*100))









