1. 程式人生 > >《計算方法》 李曉紅等 第二章 解非線性方程(組) 例題 程式碼

《計算方法》 李曉紅等 第二章 解非線性方程(組) 例題 程式碼

#include "stdafx.h"

#include <math.h>

// 解非線性方程(組)

double fvalue(double x)

{

return (x*x*x-2.5*x+1);

}

// 例1,對分法求解三次方程 x^3-2.5x+1 = 0;

//

void DuiFenMethod()

{

printf("對分法求解三次方程 x^3-2.5x+1 = 0\n");

// step 1

double a = 0.;

double b = 1.;

double epsilon = 0.5e-3;

int k = 0;

// step 2

//int N = log10(2.);

int N = int((log10(double(b-a))-log10(epsilon))/log10(2.) - 1)+1;

printf("迭代次數:%d\n",N);

double x1,f;

while(k<=N)

{

x1 = (b+a)/2;

f = fvalue(x1);

if (abs(f)<1e-6)

{

break;

}

else if(fvalue(x1)*fvalue(a) < 0)

{

b = x1;

}

else

{

a = x1;

}

printf("a%d=%6f\t",k,a);

printf("b%d=%6f\t"

,k,b);

printf("x%d=%6f\t",k,x1);

printf("f%d=%6f\n",k,f);

k++;

}

printf("root = x%d = %f\n",k-1,x1);

}

// x = sqrt(10/(4+x))

void DieDaiMethod()

{

printf("非線性方程迭代法:\n");

double x0(1.5);

double x1 = sqrt(10/(4+x0));

double epsilon = 1e-7;

double L = 0.2/sqrt(2.);

printf("x0=%.9f\n",x0);

printf

("x1=%.9f\n",x1);

printf("L=%.9f\n",L);

int n = intlog10(epsilon*(1.-L)/abs(x1-x0))/log10(L) );

printf("n=%d\n",n);

for (int k = 2; k <= 8; k++) {

x1 = sqrt(10./(4.+x1));

printf("x%d=%.9f\n",k,x1);

}

}

void AitkenMethod()

{

printf("埃特肯加速法:\n");

double a = 0.;

double b = 1.;

double epsilon = 0.5e-3;

double x0(0.5),x1;

double y,z;

int n = 0;

while (1)

{

y = (x0*x0*x0*2+2)*0.2;

z = (y*y*y*2+2)*0.2;

x1 = x0 - (y-x0)*(y-x0)/(z-2*y+x0);

if (abs(x0-x1)<epsilon)

{

printf("x%d=%.6f\n",n,x0);

break;

}

else

{

x0 = x1;

}

printf("x%d=%.6f\t",n,x0);

printf("y%d=%.6f\t",n,y);

printf("z%d=%.6f\n",n,z);

n++;

}

}

// 例4

void NewtonMethod()

{

printf("牛頓迭代法:\n");

double x0(0.);

double x1;

double epsilon = 0.5e-3;

int k = 0;

while (1)

{

printf("x%d=%.11f\n",k,x0);

x1 = (4*x0*x0*x0-2)/(6*x0*x0-5);

x1 = x0 - fvalue(x0) / (3*x0*x0-2.5);

k++;

if (abs(x0-x1)<epsilon)

{

printf("x%d=%.11f\n",k,x1); break;

}

else

x0 = x1;

}

}

void GeXianMethod()

{

printf("割線法:\n");

double x0(0.5);

double x1(0.2);

double epsilon = 0.5e-6;

int k = 0;

double x2;

while (1)

{

printf("x%d=%.6f\t",k,x0);

printf("x%d=%.6f\n",k+1,x1);

x2 = x0 - fvalue(x0)*(x1-x0)/(fvalue(x1)-fvalue(x0));

k++;

if (abs(x2-x1)<epsilon)

{

printf("x%d=%.6f\n",k,x2);  break;

}

else

{

x0 = x1; x1 = x2;

}

}

}