function nma_euler_heun2(f,rate,stepSize,startX,endX,initialY,tol) % Solve ODE using Euler-Heun (corrector-predictor method) % Function to solve ODE using Euler-Heun algorithm (called % als the corrector-predictor method. % %INPUT % f: the true solution. For debugging % % 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_heun2('2*x',rate,stepSize,startX,endX,initialY,.001) % % Author: Nasser Abbasi %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.trueY(nStep) = getTrueY(f,startStepX+stepSize); data.trueRate(nStep) = yder(rate,startStepX+stepSize,data.trueY(nStep)); data.startStepX(nStep) = startStepX; data.startStepY(nStep) = startStepY; data.startStepYRate(nStep) = yder( rate, startStepX, startStepY); data.y_predict(nStep) = y_predict; rateEndOfStep = yder( rate, endStepX, y_predict ); data.rateAtyp(nStep) = rateEndOfStep; averageRate = tan ( ( atan(rateStartOfStep) + atan(rateEndOfStep) )/2 ); % averageRate = ( rateStartOfStep + rateEndOfStep )/2; y_correct = startStepY + ( stepSize * averageRate ); data.rateAtyc(nStep) = yder( rate, endStepX, y_correct ); data.y_correct(nStep) = y_correct; relativeError = abs( (y_correct - y_predict)/y_correct )*100; data.relativeError(nStep) = relativeError; data.averageRate(nStep) = 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.trueY(nStep) = eval(f); data.trueRate(nStep) = yder(rate,x,data.trueY(nStep)); data.startStepX(nStep) = x; data.startStepY(nStep) = y; data.startStepYRate(nStep) = 0; data.y_predict(nStep) = y; data.y_correct(nStep) = y; data.relativeError(nStep) = 0; data.averageRate(nStep) = 0; data.rateAtyp(nStep) = 0; data.rateAtyc(nStep) = 0; analyze(data,stepSize); %%%%%%%%%%%%% % % %%%%%%%%%%%%%%% function v=getTrueY(f,x) v=eval(f); %%%%%%%%%%%%%%%%%%%% % % %%%%%%%%%%%%%%%%%%%% function dydx=yder(f,x,y) dydx=eval(f); %%%%%%%%%%%%%%%%%%%%%%% % % % %%%%%%%%%%%%%%%%%%%%%% function analyze(d,stepSize) row=length(d.startStepX); fprintf('n\t x\t\t y0\t\t y0''\t yp\t yp''\t yc\t\t y_c''\t y_av''\t yt\t\t y''_t\t err(a)\n'); for(i=1:row) fprintf('%d\t',i); fprintf('%4.2f\t',d.startStepX(i)); fprintf('%4.2f\t',d.startStepY(i)); fprintf('%4.2f\t',atan(d.startStepYRate(i))*180/pi); fprintf('%4.2f\t',d.y_predict(i)); fprintf('%4.2f\t',atan(d.rateAtyp(i))*180/pi); fprintf('%4.2f\t',d.y_correct(i)); fprintf('%4.2f\t',atan(d.rateAtyc(i))*180/pi); fprintf('%4.2f\t',atan(d.averageRate(i))*180/pi); fprintf('%4.2f\t',d.trueY(i)); fprintf('%4.2f\t',atan(d.trueRate(i))*180/pi ); fprintf('%4.2f\t',d.relativeError(i)); fprintf('\n'); if(i