From 3b669a2979b4aedad9a326b5b8bdd0553adf3023 Mon Sep 17 00:00:00 2001 From: Mattrixwv Date: Sat, 22 Sep 2018 15:27:54 -0400 Subject: [PATCH] Added new files for Numerical Analysis --- Bisection.m | 35 +++++++++++++++++++++++++++++++++++ FalsePosition.m | 31 +++++++++++++++++++++++++++++++ Mullers.m | 44 ++++++++++++++++++++++++++++++++++++++++++++ NewNewton.m | 25 +++++++++++++++++++++++++ Newton.m | 24 ++++++++++++++++++++++++ Secant.m | 30 ++++++++++++++++++++++++++++++ swap.m | 8 ++++++++ 7 files changed, 197 insertions(+) create mode 100644 Bisection.m create mode 100644 FalsePosition.m create mode 100644 Mullers.m create mode 100644 NewNewton.m create mode 100644 Newton.m create mode 100644 Secant.m create mode 100644 swap.m diff --git a/Bisection.m b/Bisection.m new file mode 100644 index 0000000..2dc1500 --- /dev/null +++ b/Bisection.m @@ -0,0 +1,35 @@ +function [xList,errorList] = Bisection (f, lowerValue, upperValue, allowError) +% +%Uses the bisection method to find the possible answers to the root of the function f +% + + pkg load symbolic; + warning('off','OctSymPy:sym:rationalapprox'); + %Setting necesary values for the function + cnt = 1; + maxItterations = 50; + errorValue = 1; + currentValue = 0.0; + + + %If the lower and upper bounds are mixed up swap them + if(double(subs(f,lowerValue)) > 0) + [lowerValue,upperValue] = swap(lowerValue, upperValue); + end + + %Loop until the error is within bounds or the Maximum number of iterations is reached + while((abs(errorValue) > allowError) && (cnt < maxItterations)) + currentValue = (lowerValue + upperValue)/ 2; + errorValue = double(subs(f,currentValue)); + %Replace the correct value with the new value + %if error == 0 then the value has been found + if(errorValue < 0) + lowerValue = currentValue; + else + upperValue = currentValue; + end + xList(cnt) = currentValue; + errorList(cnt) = errorValue; + ++cnt; + end +end diff --git a/FalsePosition.m b/FalsePosition.m new file mode 100644 index 0000000..1e20e19 --- /dev/null +++ b/FalsePosition.m @@ -0,0 +1,31 @@ +function [xList, errorList] = FalsePosition(f, p0, p1, errorAllowed) +% +%FalsePosition(f, p0, p1, errorAllowed) +%This function finds the root of a function using the method of False Position +% + + pkg load symbolic; + warning('off','OctSymPy:sym:rationalapprox'); + maxIt = 50; + cnt = 2; + q0 = double(subs(f,p0)); + q1 = double(subs(f,p1)); + p = 0; + q = 0; + currentError = errorAllowed + 1; + while((cnt <= maxIt) && (currentError >= errorAllowed)) + p = p1 - (q1 * (p1 - p0))/(q1 - q0); + currentError = abs(p - p1); + %Add Values to lists + xList(end+1) = p; + errorList(end+1) = currentError; + ++cnt; + q = double(subs(f,p)); + if((q * q1) < 0) + p0 = p1; + q0 = q1; + end + p1 = p; + q1 = q; + end +end diff --git a/Mullers.m b/Mullers.m new file mode 100644 index 0000000..a4b0e07 --- /dev/null +++ b/Mullers.m @@ -0,0 +1,44 @@ +function [xList, functionValueList] = Mullers(f, p0, p1, p2, errorAllowed) +% +%Mullers(f, p0, p1, p2, errorAllowed) +%This function finds the root of a function using Muller's Method +% + + pkg load symbolic; + warning('off','OctSymPy:sym:rationalapprox'); + + h1 = p1 - p0; + h2 = p2 - p1; + s1 = (double(subs(f,p1)) - double(subs(f,p0))) / h1; + s2 = (double(subs(f,p2)) - double(subs(f,p2))) / h2; + d = (s2 - s1) / (h2 + h1); + cnt = 2; + maxIt = 50; + + while(cnt < maxIt) + b = s2 + h2 * d; + D = (b^2 - (4 * double(subs(f,p2)) * d))^(1/2); + if(abs(b - D) < abs(b + D)) + E = b + D; + else + E = b - D; + end + h = (02 * double(subs(f,p2)))/E; + p = p2 + h; + xList(end+1) = p; + functionValueList(end+1) = double(subs(f,p)); + if(abs(h) < errorAllowed) + return; + else + p0 = p1; + p1 = p2; + p2 = p; + h1 = p1 - p0; + h2 = p2 - p1; + s1 = (double(subs(f,p1)) - double(subs(f,p0)))/h1; + s2 = (double(subs(f,p2)) - double(subs(f,p1)))/h2; + d = (s2 - s1)/(h2 + h1); + ++cnt; + end + end +end diff --git a/NewNewton.m b/NewNewton.m new file mode 100644 index 0000000..e239746 --- /dev/null +++ b/NewNewton.m @@ -0,0 +1,25 @@ +function [xList, errorList] = NewNewton (f, startingValue, errorAllowed) +% +%NewNewton (f, startingValue, errorAllowed) +%This function computes the root of a function using the modified Newton's method +% + + pkg load symbolic; + warning('off','OctSymPy:sym:rationalapprox'); + oldAnswer = startingValue; + newAnswer = 0; + currentError = errorAllowed + 1; + cnt = 1; + maxIt = 50; + fp = diff(f); + fpp = diff(fp); + + + while((currentError >= errorAllowed) && (cnt < maxIt)) + newAnswer = oldAnswer - ((double(subs(f,oldAnswer)) * double(subs(fp,oldAnswer)))/(double(subs(fp,oldAnswer))^2 - (double(subs(f,oldAnswer)) * double(subs(fpp,oldAnswer))))); + currentError = abs(newAnswer - oldAnswer); + xList(end+1) = newAnswer; + errorList(end+1) = currentError; + oldAnswer = newAnswer; + end +end diff --git a/Newton.m b/Newton.m new file mode 100644 index 0000000..1d87775 --- /dev/null +++ b/Newton.m @@ -0,0 +1,24 @@ +function [xList,errorList] = Newton (f, startingValue, errorAllow) +% +%This function uses Newtons method to find a solution to the root of f +% + + pkg load symbolic; + warning('off','OctSymPy:sym:rationalapprox'); + maxIt = 50; + fp = diff(f); + oldAnswer = startingValue; + newAnswer = 0; + cnt = 1; + currentError = errorAllow + 1; + + %Loop until the error becomes small enough or the maximum number of itterations is met + while((cnt < maxIt) && (currentError > errorAllow)) + newAnswer = oldAnswer - (double(subs(f,oldAnswer))/double(subs(fp,oldAnswer))); + currentError = abs(newAnswer - oldAnswer); + xList(end+1) = newAnswer; + errorList(end+1) = currentError; + oldAnswer = newAnswer; + ++cnt; + end +end diff --git a/Secant.m b/Secant.m new file mode 100644 index 0000000..5662b4c --- /dev/null +++ b/Secant.m @@ -0,0 +1,30 @@ +function [xList, errorList] = Secant(f, p0, p1, errorAllowed) +% +%Secant(f, p0, p1, errorAllowed) +%This function find the root of a function using the Secant Method +% + + pkg load symbolic; + warning('off','OctSymPy:sym:rationalapprox'); + maxIt = 50; + cnt = 2; + q0 = double(subs(f,p0)); + q1 = double(subs(f,p1)); + currentError = errorAllowed + 1; + p = 0; + + + while((cnt <= maxIt) && (currentError >= errorAllowed)) + p = p1 - (q1 * (p1 - p0))/(q1 - q0); + currentError = abs(p - p1); + %Add the x and error values to memory + xList(end+1) = p; + errorList(end+1) = currentError; + %Setup for the next run + ++cnt; + p0 = p1; + q0 = q1; + p1 = p; + q1 = double(subs(f,p)); + end +end diff --git a/swap.m b/swap.m new file mode 100644 index 0000000..b667a62 --- /dev/null +++ b/swap.m @@ -0,0 +1,8 @@ +function [x,y] = swap(first, second) +% +%Swap the first and second values +% + + x = second; + y = first; +end