function nma_RK4(rate,stepSize,startX,endX,initialY) %solve 1st order ODE using Runge-Kutta classical 4th order %function nma_RK4(rate,stepSize,initialX,finalX) % % Function to solve 1st order ODE using Runge-Kutta classical 4th order % %INPUT % rate: A string that represents the derivative dy/dx. This % is the function f(x,y) in the ODE equation dy/dx= f(x,y) % This string MUST use the letters 'x' and 'y' for the % independent and the dependent variables respectively. % For an example, in the equation dy/dx = y*x, the rate % string is passed in as 'y*x' % % stepSize: This is 'h', the step size over x % % startX: The starting value of x % endX: The ending value of x % initialY: Initial condition. The value of y at startX. % % EXAMPLE: % rate='x^2'; stepSize=.1; startX=0; % endX=1; initialY=1; nma_RK4(rate,stepSize,startX,endX,initialY) % step x y(x) end_step_y % 1 0.0000 1.0000 1.0003 % 2 0.1000 1.0003 1.0027 % 3 0.2000 1.0027 1.0090 % 4 0.3000 1.0090 1.0213 % 5 0.4000 1.0213 1.0417 % 6 0.5000 1.0417 1.0720 % 7 0.6000 1.0720 1.1143 % 8 0.7000 1.1143 1.1707 % 9 0.8000 1.1707 1.2430 % 10 0.9000 1.2430 1.3333 % 11 1.0000 1.3333 1.3333 % % copyright: Nasser M. Abbasi y = initialY; x = startX; nStep = 0; nPoints = (endX-startX)/stepSize +1; currentPoint = 1; data = zeros(nPoints,4); while(currentPoint < nPoints) startStepX = x; startStepY = y; k1 = yder(rate,startStepX ,startStepY ); k2 = yder(rate,startStepX+stepSize/2,startStepY+(1/2)*stepSize*k1 ); k3 = yder(rate,startStepX+stepSize/2,startStepY+(1/2)*stepSize*k2 ); k4 = yder(rate,startStepX+stepSize ,startStepY+stepSize*k3 ); nStep = nStep+1; data(nStep,1) = nStep; data(nStep,2) = startStepX; data(nStep,3) = startStepY; y = startStepY + stepSize*( (1/6) * (k1+ 2*k2+ 2*k3+ k4) ); data(nStep,4) = y; currentPoint = currentPoint+1; x = x + stepSize; end nStep = nStep+1; data(nStep,1) = nStep; data(nStep,2) = x; data(nStep,3) = y; data(nStep,4) = y; analyze(data); %my_test(data,stepSize); end %%%%%%%%%%%%%%%%%%%% % % %%%%%%%%%%%%%%%%%%%% function dydx=yder(f,x,y) dydx=eval(f); end %%%%%%%%%%%%%%%%%%%%%%% % % % %%%%%%%%%%%%%%%%%%%%%% function analyze(d) [row,~]=size(d); fprintf('step\t x\t\t y(x) \t end_step_y\n'); for i=1:row fprintf('%d\t %7.4f\t %7.4f\t %7.4f\n',... d(i,1),d(i,2),d(i,3),d(i,4)); end end %%%%%%%%%%%%%%%%%% % %%%%%%%%%%%%%%%%%% function my_test(d,stepSize) plotAnalyticAndRKOnly(d,stepSize); figure; x_5=[0 0.5 1 1.5 2]; y_5=[1 0.4 0.21 0.189 0.288225]; x_25=[0 0.25 0.5 0.75 1 1.25 1.5 1.75 2]; y_25=[1 0.7 0.5009375 0.38196484375 0.321089196777344 ... 0.305034736938477 0.332678509973526 0.332678509973526 0.487581941179949]; y_analytic='exp(x*(x^2/3 - 1.2))'; ezplot(y_analytic,0,2); hold on; plot(x_5,y_5,'o') plot(x_5,y_5,'-.') plot(x_25,y_25,'x') plot(x_25,y_25,'--') plot(d(:,2),d(:,3),'.:') xlim([-0.1 2.1]); ylim([0 2.2]) legend('analytic','h=0.5 pts','h=0.5','h=0.25 pts','h=0.25','RK4 h=0.5'); ylabel('y'); title(sprintf('analytic solution of %s compared to RK4 (h=%.2f)',... y_analytic,stepSize)); end %%%%%%%%%%%%%%%%%%%% % % %%%%%%%%%%%%%%%%%%%%%%%% function plotAnalyticAndRKOnly(cleanData,stepSize) figure; y_analytic='exp(x*(x^2/3 - 1.2))'; ezplot(y_analytic,0,2); hold on; plot(cleanData(:,2),cleanData(:,3),'.:') xlim([-0.1 2.1]); ylim([0 2.2]) legend('analytic','RK4 h=0.5'); ylabel('y'); title(sprintf('analytic solution of %s compared to RK4(h=%.2f)',... y_analytic,stepSize)); end