Matlab 4th order PDF

Preview

Summary

lecture note...

Description

clc clf clear all % Mfile name % mtl_int_sim_RK4thmethod.m % Revised: % March 7, 2008 % % Authors % Nathan Collier, Autar Kaw % Department of Mechanical Engineering % University of South Florida % Tampa, Fl 33620 % Contact: [email protected] | http://numericalmethods.eng.usf.edu/contact.html % Purpose % To illustrate the 4th 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 ; % n, number of steps to take n=5 ; %********************************************************************** % Displays title information disp(sprintf('\n\nThe 4th Order Runge-Kutta Method of Solving Ordinary Differential Equations')) disp(sprintf('University of South Florida'))

disp(sprintf('United States of America')) disp(sprintf('[email protected]\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 title information disp('The 4th 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).') % Displays inputs that will be 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(' 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(' n = %g',n)) disp(sprintf('\n----------------------------------------------------------------')) disp(sprintf('For this simulation, the following parameter is constant.\n')) h=(xf-x0)/n ; disp(sprintf(' h = ( xf - x0 ) / n ')) disp(sprintf(' = ( %g - %g ) / %g ',xf,x0,n)) disp(sprintf(' = %g',h)) xa(1)=x0 ; ya(1)=y0 ; % The simulation of the method begins here. 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, k2, k3, and k4 k1 = f(xa(i),ya(i)) ; k2 = f(xa(i)+0.5*h,ya(i)+0.5*k1*h) ; k3 = f(xa(i)+0.5*h,ya(i)+0.5*k2*h) ;

k4 = f(xa(i)+h,ya(i)+k3*h) ; % Using 4th Order Runge-Kutta formula ya(i+1)=ya(i)+1/6*(k1+2*k2+2*k3+k4)*h ; disp('1) Find k1 and k2 disp(sprintf(' k1 = disp(sprintf(' = disp(sprintf(' =

using the previous step information.') f( x%g , y%g )',i-1,i-1)) f( %g , %g )',xa(i),ya(i))) %g\n',k1))

disp(sprintf(' k2 = 1)) disp(sprintf(' = %g)',xa(i),h,ya(i),k1,h)) disp(sprintf(' = disp(sprintf(' =

f( x%g + 0.5 * h , y%g + 0.5 * k1 * h )',i-1,i-

disp(sprintf(' k3 = 1)) disp(sprintf(' = %g)',xa(i),h,ya(i),k2,h)) disp(sprintf(' = disp(sprintf(' =

f( x%g + 0.5 * h , y%g + 0.5 * k2 * h )',i-1,i-

disp(sprintf(' disp(sprintf(' disp(sprintf('

f( %g + 0.5 * %g , %g + 0.5 * %g * f( %g , %g )',xa(i)+0.5*h,ya(i)+0.5*k1*h)) %g\n',k2))

f( %g + 0.5 * %g , %g + 0.5 * %g * f( %g , %g )',xa(i)+0.5*h,ya(i)+0.5*k2*h)) %g\n',k3))

k4 = f( x%g + h, y%g + k3*h)',i-1,i-1)) = f( %g , %g )',xa(i)+h,ya(i)+k3*h)) = %g\n',k1))

disp(sprintf('2) Apply the Runge-Kutta 4th Order method to estimate y %g',i)) disp(sprintf(' y%g = y%g + 1/6*( k1 + 2*k2 + 2*k3 +k4) * h',i,i-1)) disp(sprintf(' = %g + %g * %g * %g',ya(i),1/6, (k1+2*k2+2*k3+k4),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 4th Order RungeKutta 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'))...