function nma_euler_midpoint(rate,stepSize,startX,endX,initialY) % solve ODE using Euler-mid-point algorithm % %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. % % % Author: Nasser Abbasi y = initialY; x = startX; nStep = 0; nPoints = (endX-startX)/stepSize +1; currentPoint = 1; while(currentPoint < nPoints) startStepX = x; endStepX = x + stepSize; startStepY = y; rateStartOfStep = yder( rate, startStepX, startStepY ); y_midpoint = startStepY + ( stepSize/2 * rateStartOfStep ); nStep = nStep+1; data(nStep,1) = nStep; data(nStep,2) = startStepX; data(nStep,3) = startStepY; rateMidPoint = yder( rate, startStepX + stepSize/2 , y_midpoint ); data(nStep,4) = rateMidPoint; y = startStepY + ( stepSize * rateMidPoint ); data(nStep,5) = 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; data(nStep,5) = y; analyze(data,stepSize); %%%%%%%%%%%%%%%%%%%% % % %%%%%%%%%%%%%%%%%%%% function dydx=yder(f,x,y) dydx=eval(f); %%%%%%%%%%%%%%%%%%%%%%% % % % %%%%%%%%%%%%%%%%%%%%%% function analyze(d,stepSize) [row,col]=size(d); fprintf('step\t x\t\t y(x) \t midpoint_rate end_step_y\n'); N=0; for(i=1:row) fprintf('%d\t %7.4f\t %7.4f\t %7.4f\t %7.4f\n',... d(i,1),d(i,2),d(i,3),d(i,4),d(i,5)); end plotAnalyticAndMidPointOnly(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','midpoint h=0.5'); ylabel('y'); title(sprintf('analytic solution of %s compared to Euler and midpoint(h=%.2f)',y_analytic,stepSize)); %%%%%%%%%%%%%%%%%%%% % % %%%%%%%%%%%%%%%%%%%%%%%% function plotAnalyticAndMidPointOnly(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','MidPoint h=0.5'); ylabel('y'); title(sprintf('analytic solution of %s compared to MidPoint(h=%.2f)',y_analytic,stepSize));