diff --git a/hws/hw2/main.c b/hws/hw2/main.c index 26efcd4..d940ede 100644 --- a/hws/hw2/main.c +++ b/hws/hw2/main.c @@ -2,6 +2,7 @@ #include #include "nr.h" +#include "rtmuller.h" void dbessj0(float x, float *y, float *dy) { *y = bessj0(x); @@ -25,16 +26,17 @@ int main() { int i; zbrak(bessj0, a, b, n, xb1, xb2, &nb); - for (i = 1; i <= nb; i++) { + /* for (i = 1; i <= nb; i++) { printf("[%f, %f]\n", xb1[i], xb2[i]); - } + } */ float *roots_bisect = calloc(sizeof(float), nb); float *roots_linter = calloc(sizeof(float), nb); float *roots_secant = calloc(sizeof(float), nb); float *roots_newt = calloc(sizeof(float), nb); float *roots_newtsafe = calloc(sizeof(float), nb); - + float *roots_muller = calloc(sizeof(float), nb); + for (i = 0; i < nb; i++) { float l = xb1[i + 1]; float r = xb2[i + 1]; @@ -43,6 +45,7 @@ int main() { roots_secant[i] = rtsec(bessj0, l, r, xacc); roots_newt[i] = rtnewt(dbessj0, l, r, xacc); roots_newtsafe[i] = rtsafe(dbessj0, l, r, xacc); + roots_muller[i] = rtmuller(bessj0, l, r, xacc); } printf("besection roots: "); @@ -75,6 +78,12 @@ int main() { } printf("\n"); + printf("muller method roots: "); + for (i = 0; i < nb; i++) { + printf("%f ", roots_muller[i]); + } + printf("\n"); + /* free */ free(xb1); free(xb2); @@ -83,6 +92,8 @@ int main() { free(roots_linter); free(roots_secant); free(roots_newt); + free(roots_newtsafe); + free(roots_muller); return 0; } \ No newline at end of file diff --git a/hws/hw2/muller.c b/hws/hw2/muller.c deleted file mode 100644 index e69de29..0000000 diff --git a/hws/hw2/muller.h b/hws/hw2/muller.h deleted file mode 100644 index e69de29..0000000 diff --git a/hws/hw2/rtmuller.c b/hws/hw2/rtmuller.c new file mode 100644 index 0000000..3004228 --- /dev/null +++ b/hws/hw2/rtmuller.c @@ -0,0 +1,51 @@ +#include "nrutil.h" + +#include "rtmuller.h" +#include + +float rtmuller(float (*func)(float), float x1, float x2, float xacc) { + int j; + float fh, fl, fm, fnew, s, xl, r, x3, x0, x_prev, x_curr, x_next; + float d0, d1, d2, a, b, c; + + x0 = x1; + x_prev = x2; + x_curr = (x0 + x_prev) / 2.0; + + if (x_prev == x_curr) x_curr += xacc; + if (x0 == x_curr) x0 -= xacc; + + for (j = 1; j <= 200; j++) { + d0 = func(x0); + d1 = func(x_prev); + d2 = func(x_curr); + + float h1 = x_prev - x0; + float h2 = x_curr - x_prev; + float delta1 = (d1 - d0) / h1; + float delta2 = (d2 - d1) / h2; + a = (delta2 - delta1) / (h2 + h1); + b = a * h2 + delta2; + c = d2; + + float disc = b * b - 4 * a * c; + if (disc < 0) nrerror("Muller's method: complex roots encountered"); + + float denom1 = b + sqrt(disc); + float denom2 = b - sqrt(disc); + + if (fabs(denom1) < fabs(denom2)) { + x_next = x_curr - (2 * c) / denom2; + } else { + x_next = x_curr - (2 * c) / denom1; + } + + if (fabs(x_next - x_curr) < xacc) return x_next; + + x0 = x_prev; + x_prev = x_curr; + x_curr = x_next; + } + nrerror("Too many iterations in rtmuller"); + return 0.0; +} \ No newline at end of file diff --git a/hws/hw2/rtmuller.h b/hws/hw2/rtmuller.h new file mode 100644 index 0000000..e49e0fb --- /dev/null +++ b/hws/hw2/rtmuller.h @@ -0,0 +1,3 @@ +#pragma once + +float rtmuller(float (*func)(float), float x1, float x2, float xacc);