《計算方法》 李曉紅等 第二章 解非線性方程(組) 例題 程式碼
#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 = int( log10(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;
}
}
}