mtl_ode_sim_rk2ndmethod.m

CFD（计算流体力学）中ODE下RK（RungeKutta）2的计算编程方法。

clc
clf
clear all

% Mfile name
%       mtl_int_sim_RK2ndmethod.m

% Revised:
%       March 7, 2008

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

% Purpose
%       To illustrate the 2nd order Runge-Kutta method applied
%       to a function of the user's choosing.

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

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

% dy/dx in form of f(x,y). In general it can be a function of both 
% variables x and y. If your function is only a function of x then
% you will need to add a 0*y to your function.

   fcnstr='y*x^2-1.2*y' ;
   f=inline(fcnstr) ;

% 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=2 ;

% a2, parameter which must be between 0 and 1. Certain names are associated
% with different parameters.
% 
% a2 = 0.5 Heun's Method
%    = 2/3 Ralston's Method
%    = 1.0 Midpoint Method

   a2=0.5 ;

% n, number of steps to take

   n=5 ;

%**********************************************************************

% Displays title information
disp(sprintf('\n\nThe 2nd Order Runge-Kutta Method of Solving Ordinary Differential Equations'))
disp(sprintf('University of South Florida'))
disp(sprintf('United States of America'))
disp(sprintf('[e-mail]\n'))
disp('NOTE: This worksheet demonstrates the use of Matlab to illustrate ')
disp('the Runge-Kutta method, a numerical technique in solving ordinary ')
disp('differential equations.') 

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

% Displays introduction text 
disp('The 2nd Order Runge-Kutta method approximates the solution to an ')
disp('ordinary differential equation by using the equation expressed in')
disp('the form dy/dx = f(x,y) to approximate the slope. This slope is used')
disp('to project the solution to the ODE a fixed distance away.')

% displays what inputs are used
disp(sprintf('\n\n***************************Input Data*******************************'))
disp('Below are the input parameters to begin the simulation which can be')
disp('changed in the m-file input section')
disp(sprintf('\n     f = dy/dx ')) 
disp(sprintf('     x0 = initial x ')) 
disp(sprintf('     y0 = initial y ')) 
disp(sprintf('     xf = final x '))
disp(sprintf('     a2 = constant value between 0 and 1.'))
disp(sprintf('        = 0.5, Heun Method'))
disp(sprintf('        = 1.0, Midpoint Method'))
disp(sprintf('        = 0.66667, Ralston''s Method')) 
disp(sprintf('      n = number of steps to take')) 
format short g
disp(sprintf('\n-----------------------------------------------------------------\n'))
disp(sprintf(['     f(x,y) = dy/dx = ' fcnstr]))
disp(sprintf('     x0 = %g',x0))
disp(sprintf('     y0 = %g',y0))
disp(sprintf('     xf = %g',xf))
disp(sprintf('     a2 = %g',a2))
disp(sprintf('     n = %g',n))
disp(sprintf('\n-----------------------------------------------------------------'))
disp(sprintf('For this simulation, the following parameter is constant.\n'))

% Calculates constants used in the method
h=(xf-x0)/n ;
disp(sprintf('     h = ( xf - x0 ) / n '))
disp(sprintf('       = ( %g - %g ) / %g ',xf,x0,n))
disp(sprintf('       = %g',h))
a1=1-a2 ;
disp(sprintf('\n    a1 = 1 - a2'))
disp(sprintf('       = 1 - %g',a2))
disp(sprintf('       = %g',a1))
p1=1/2/a2 ;
disp(sprintf('\n    p1 = 1 / ( 2 * a2 )'))
disp(sprintf('       = 1 / ( 2 * %g )',a2))
disp(sprintf('       = %g',p1))
q11=p1 ;
disp(sprintf('\n   q11 = p1'))
disp(sprintf('       = %g',q11))

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

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

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

% Adding Step Size
  xa(i+1)=xa(i)+h ;

% Calculating k1 and k2
  k1 = f(xa(i),ya(i)) ;
  k2 = f(xa(i)+p1*h,ya(i)+q11*k1*h) ;

% Using 2nd Order Runge-Kutta formula
  ya(i+1)=ya(i)+(a1*k1+a2*k2)*h ;

  disp('1) Find k1 and k2 using the previous step information.')
  disp(sprintf('     k1 = f( x%g , y%g )',i-1,i-1))
  disp(sprintf('        = f( %g , %g ) )',xa(i),ya(i)))
  disp(sprintf('        = %g\n',k1))
  disp(sprintf('     k2 = f( x%g + p1 * h , y%g + q11 * k1 * h )',i-1,i-1))
  disp(sprintf('        = f( %g + %g * %g , %g + %g * %g * %g)',xa(i),p1,h,ya(i),q11,k1,h))
  disp(sprintf('        = f( %g , %g )',xa(i)+p1*h,ya(i)+q11*k1*h))
  disp(sprintf('        = %g\n',k2))

  disp(sprintf('2) Apply the Runge-Kutta 2nd Order method to estimate y%g',i))
  disp(sprintf('     y%g = y%g + ( a1 * k1 + a2 * k2 ) * h',i,i-1))
  disp(sprintf('        = %g + %g * %g',ya(i),a1*k1+a2*k2,h))
  disp(sprintf('        = %g\n',ya(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
xspan = [x0 xf];
[x,y]=ode45(f,xspan,y0);
[yfi dummy]=size(y);
yf=y(yfi);

% Plotting the Exact and Approximate solution of the ODE.
hold on
xlabel('x');ylabel('y');
title('Exact and Approximate Solution of the ODE by the 2nd Order Runge-Kutta Method');
plot(x,y,'--','LineWidth',2,'Color',[0 0 1]);            
plot(xa,ya,'-','LineWidth',2,'Color',[0 1 0]);
legend('Exact','Approximation');

disp('The figure window that now appears shows the approximate solution as ') 
disp('piecewise continuous straight lines. The blue line represents the exact')
disp('solution. In this case ''exact'' refers to the solution obtained by the')
disp(sprintf('Matlab function ode45.\n'))

disp('Note the approximate value obtained as well as the true value and ')
disp('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))

disp(sprintf('\nThe approximation can be more accurate if we made our step'))
disp(sprintf('size smaller (that is, increasing the number of steps).\n\n'))