function nma_euler_heun(rate,stepSize,startX,endX,initialY,tol) % Solve ODE using Euler-Heun (corrector-predictor method) %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. % tol: In percentage. The value to stop the corrector-predictor % loop at each step. % % EXAMPLE % rate='x^2'; stepSize=.1; startX=0; endX=1; initialY=1; % nma_euler_heun(rate,stepSize,startX,endX,initialY,.001) % % copyright: Nasser M. Abbasi % 2003 % This is protect again predictor-correctpr loop not converging %to with epsilon. MAX_STEPS_PER_ONE_ITERATION=1000; y = initialY; x = startX; nStep = 0; nPoints = (endX-startX)/stepSize +1; currentPoint = 1; while(currentPoint < nPoints) converged = false; startStepX = x; endStepX = x + stepSize; startStepY = y; rateStartOfStep = yder( rate, startStepX, startStepY ); y_predict = startStepY + ( stepSize * rateStartOfStep ); thisIterStep = 0; while(~converged) nStep = nStep+1; data(nStep,1) = nStep; data(nStep,2) = startStepX; data(nStep,3) = startStepY; data(nStep,4) = y_predict; rateEndOfStep = yder( rate, endStepX, y_predict ); averageRate = ( rateStartOfStep + rateEndOfStep )/2; y_correct = startStepY + ( stepSize * averageRate ); data(nStep,5) = y_correct; relativeError = abs( (y_correct - y_predict)/y_correct )*100; data(nStep,6) = relativeError; data(nStep,7) = averageRate; y_predict = y_correct; thisIterStep = thisIterStep+1; if(relativeErrorMAX_STEPS_PER_ONE_ITERATION) converged=true; end end % Update for next point x = x+stepSize; y = y_correct; currentPoint = currentPoint+1; end nStep = nStep+1; data(nStep,1) = nStep; data(nStep,2) = x; data(nStep,3) = y; data(nStep,4) = y; data(nStep,5) = y; data(nStep,6) = 0; data(nStep,7) = 0; analyze(data,stepSize); %%%%%%%%%%%%%%%%%%%% % % %%%%%%%%%%%%%%%%%%%% function dydx=yder(f,x,y) dydx=eval(f); %%%%%%%%%%%%%%%%%%%%%%% % % % %%%%%%%%%%%%%%%%%%%%%% function analyze(d,stepSize) [row,~]=size(d); fprintf('step\t x\t\t y(x) \t predict_y(x+h) correct_y(x+h) \tepsilon(a)\n'); N=0; for i=1:row if(i>1) if(d(i-1,2) == d(i,2)) fprintf('%d\t %7.4f\t %7.4f\t\t %7.4f\t%7.4f\n',... d(i,1),d(i,4),d(i,5),d(i,6),d(i,7)); else N=N+1; cleanData(N,1)=d(i-1,2); cleanData(N,2)=d(i-1,3); fprintf('%d\t %7.4f\t %7.4f\t %7.4f\t %7.4f\t\t %7.4f\t%7.4f\n',... d(i,1),d(i,2),d(i,3),d(i,4),d(i,5),d(i,6),d(i,7)); end else fprintf('%d\t %7.4f\t %7.4f\t %7.4f\t %7.4f\t\t %7.4f\t%7.4f\n',... d(i,1),d(i,2),d(i,3),d(i,4),d(i,5),d(i,6),d(i,7)); end if(i