Upload New File

mysolver.m

+%
+% Copyright 2019-2020 LIST (Luxembourg Institute of Science and Technology), all right reserved.
+%
+% Author: Nicolas Rodriguez (nicolas.bjorn.rodriguez@gmail.com)
+%
+function [result,residue,iter,conv,exitflag]=mysolver(eq,a,b,tol,prec,maxiter,display)
+% Solver by dichotomy (bisection method)
+
+% INPUT:
+%   eq: function handle to be used in feval, represents the equation to be solved and must equal zero
+%   a,b, respectively lower and upper interval limits in which the root is searched
+%   tol: approximation of zero
+%   prec: convergence criteria on the solution / precision of the result
+%   maxiter: max number of iterations
+%   display: triggers the display of notifications (on or off)
+
+% OUTPUT:
+%   result: value of the root
+%   residue: value of the function at result
+%   iter: number of iterations used
+%   conv: c value difference between last iteration and the one just before
+%   exitflag:   0 = initialization failed
+%               1 = no sign change found, no calculation done
+%               2 = solution is a or b
+%               3 = result given at max iteration number
+%               4 = solution was found by tolerance
+%               5 = solution was found by precision
+
+result=NaN;
+residue=NaN;
+iter=NaN;
+conv=NaN;
+exitflag=0;
+
+if (isnan(eq(a)) || isnan(eq(b)))
+    if strcmp(display,'on')
+        warning('mysolver: wrong initial search interval');
+    end
+    return;
+end
+
+if strcmp(display,'on')
+    if prec<tol
+        warning('mysolver: the required precision on the root is lower than the tolerance on zero');
+    end
+end
+
+if abs(eq(a))<=tol
+    exitflag=2;
+    residue=eq(a);
+    result=a;
+    iter=0;
+    conv=0;
+    return;
+end
+if abs(eq(b))<=tol
+    exitflag=2;
+    residue=eq(b);
+    result=b;
+    iter=0;
+    conv=0;
+    return;
+end
+
+if eq(a)*eq(b)>0
+    exitflag=1;
+    if strcmp(display,'on')
+        warning('mysolver: no solution in the proposed search interval');
+    end
+    result=NaN;
+    residue=NaN;
+    iter=NaN;
+    conv=NaN;
+    return;
+end
+
+iter=0;
+conv=2*prec+tol;
+while conv>=prec
+    if iter>0
+        conv=abs((a+b)/2-c);
+    end
+    c=(a+b)/2;
+    if abs(eq(c))<=tol
+        exitflag=4;
+        residue=eq(c);
+        result=c;
+        iter=iter+1;
+        return;
+    end
+    if eq(a)*eq(c)<0
+        b=c;
+    else
+        a=c;
+    end
+        
+    iter=iter+1;
+    if iter+1>maxiter
+        exitflag=3;
+        residue=eq(c);
+        result=c;
+        if strcmp(display,'on')
+            warning('mysolver: maximum number of iterations was reached');
+        end
+        return;
+    end
+end
+
+if conv<prec
+    exitflag=5;
+    result=c;
+    residue=eq(result);
+    iter=iter+1;
+    return;
+end
+
+end