原创作品,出自 “晓风残月xj” 博客,欢迎转载,转载时请务必注明出处(http://blog.csdn.net/xiaofengcanyuexj)。
由于各种原因,可能存在诸多不足,欢迎斧正!
二、引言
计算机不具有人的思维,那么如果想要求得一个方程的精确解,需要将这类方程的求根公式输入计算机,然而多数方程不存在求根公式,因此求精确根非常困难,甚至不可能,从而寻找方程的近似根就显得特别重要。
下面总结求得方程近似解的算法,虽然每一步计算过程的具体运算复杂,但每一步的方式相同,所以可以通过编写程序,利用循环和迭代让计算机完成。
三、算法描述
1、二分法(Bisection)
原函数:f(x)=cosx-0.5
(1)初始化x1和x2分别为0和π/2,f(x1)和f(x2)显然异号。说明(x1,x2)内一定存在0点。
(2)求出f[(x1+x2)/2]的值,并求出真实解π/3与(x1+x2)/2的误差。
如果f(x1)*f[(x1+x2)/2]大于0,则0点在(x1+x2)/2的右侧,这种情况下令x1=(x1+x2)/2。
如果f(x1)*f[(x1+x2)/2]大于0,则0点在(x1+x2)/2的右侧,这种情况下令x1=(x1+x2)/2。
如果f(x1)*f[(x1+x2)/2]=0,则0点就是(x1+x2)/2。
(3)重复(2)—(3)步直到满意为止。
2、非线性内插法(Linear-interpolation)
原函数:f(x)=cosx-0.5
(1)初始化x1和x2分别为为0和π/2。
(2)求出(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1))的值。计算真实解与x0的误差。
如果f(x1)*(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1))大于0,则0点在(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1))的右侧,在这种情况下令(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1))。
如果f(x1)*(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1))小于0,则0点在(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1))的左侧,在这种情况下令(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1))。
如果f(x1)*(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1))=0,则0点就是(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1))。
(3)重复(2)—(3)步直到满意为止。
3、牛顿迭代法(Nowton raphson)
原函数:f(x)=cosx-0.5
(1)初始化x0为π/2,x0显然在0点的右侧。
(2)求出x0-(f(x0)/(-sin(x0)))的值,并赋给x0。计算真实解与x0的误差。
(3)重复(2)—(3)步直到满意为止。
4、截距法(Secant)
原函数:f(x)=cosx-0.5
(2)初始化x0和x1分别为为0和π/2,f(x1)和f(x2)显然异号。说明(x0,x2)内一定存在0点。
(2)求出x1-(x1-x0)*f(x1)/(f(x1)-f(x0))的值,并赋给x0。计算真实解与x0的误差。
(3)求出x1-(x1-x0)*f(x1)/(f(x1)-f(x0))的值,并赋给x1。计算真实解与x1的误差。
(4)重复(2)—(4)步直到满意为止。
5、直接迭代法I(Direct iteration1)
原函数:f(x)=cosx-0.5
(1)令g(x)=x+cosx-0.5
(2)初始化x0=0。
(3)求出g(x0)的值,并赋给x0。计算真实解与x0的误差。
(4)重复(2)—(4)步直到满意为止。
6、直接迭代法II(Direct iteration2)
原函数:f(x)=cosx-0.5
(1)令g(x)=2x*cosx
(2)初始化x0=0。
(3)求出g(x0)的值,并赋给x0。计算真实解与x0的误差。
(4)重复(2)—(4)步直到满意为止。
7、直接迭代法III(Direct iteration3)
原函数:f(x)=cosx-0.5
(1)令g(x)=4*x*(x-π)/(π*π)
(2)初始化x0=π/2。
(3)求出x0-f(x0)/g(x0)的值,并赋给x0。计算真实解与x0的误差。
(4)重复(2)—(4)步直到满意为止。
五、实验及分析
数据分析:
上述实验数据在我们的个人计算机数据范围,精度的因素影响下基本与教材一致。每次迭代的前几次结果几乎完全一样,足以说明程序的正确性。后几次往往由于数据太小,在有限字节内在存储,运算等方面存在一定误差,属系统误差,排除算法实现上的硬伤。
有必要提出,Direct Interration2 的初值不是书上所说的PI/2,而是0.0635232+PI/3,这点我们在实现算法时深受其害,代初值PI/2是得不到书上的测试数据的。0.0635232+PI/3是由x-x*=e这个数学关系得到的。
而且Direct Interration3的初值0很明显是有问题的,迭代式
X[n+1]=g(X[n])=X[n]-f(X[n])/h(X[n])
其中h(x)=4*x*(x-PI)/PI^2(由泰勒展开式得)
在第一次迭代过程中h(x)=h(0)=0,分母为0,这是不合法的。然而通过测试x[0]=PI/2正好与数据吻合。
六、我们的特色
1.界面:
我觉得我们的界面比较美观,如下
提供menu菜单供选择,其中1-7为单个算法演示,8为所有 7个程序的演示,9为退出,实现了人机互动,操作方便,便于门 外汉使用。
2.健壮性
我们的程序还有输入出错处理机制,在强调正确的同时,保证 了程序的健壮性,如:
有出错提示。
代码如下,通过相应输入调用相关函数
while(1)
{
tag=0;
do{
if(tag!=0)
cout<<"Your choice is illegal!Please input again!"<<endl;
tag++;
cout<<" =============Menu============="<<endl;
cout<<" || 1.Bisection|| ||"<<endl;
cout<<" || 2.linearinterpolation ||"<<endl;
cout<<" || 3.Nowton Raphson ||"<<endl;
cout<<" || 4.Secant ||"<<endl;
cout<<" || 5.Direct iteration1 ||"<<endl;
cout<<" || 6.Direct iteration2 ||"<<endl;
cout<<" || 7.Direct iteration3 ||"<<endl;
cout<<" || 8.all of seven ||"<<endl;
cout<<" || 9.return ||"<<endl;
cout<<" ==============================="<<endl;
cout<<" Please input your choice[ ]\b\b";
cin>>choice;
cout<<endl<<endl;
}while(choice<1||choice>9);
if(1==choice||8==choice)Bisection();
if(2==choice||8==choice)linearinterpolation();
if(3==choice||8==choice)NowtonRaphson();
if(4==choice||8==choice)Secant();
if(5==choice||8==choice)Directiteration1();
if(6==choice||8==choice)Directiteration2();
if(7==choice||8==choice)Directiteration3();
if(9==choice)exit(0);
}
简洁直观。
3.变量设置
我们将程序中经常用到的变量设为全局变量,避免过多的耗掉存储空间,同时也避免了频繁传参带来的不必要麻烦,double x0,x1,x2,x[20],e[20];当然这一切都是在避免冲突,保证正确性、健壮性的前提下完成的。
4.模块化
虽说算法比较简单,程序比较小,但我们小组还是采取了模块化编程,具体来说也就是少写复杂语句,多调用函数,一来程序看起来比较清晰(无论思路还是视觉上,也无论是对读者还是读编程者),再者有利于合作,让小组成员都参与进来,明确分工,紧密合作。如
我们的程序就尽量回避复杂语句,多调用函数,并把多次用到的输出功能单独写成print函数。
我们的函数清单:
double f(double x) //计算cos(x)-0.5
void print1(int n) //打印输出1
void print2(int n) //打印输出2
void print3(int n) //打印输出3
以下7个函数就么必要一一注释了,我们的函数还有变量名是尽量做到见名知意,这也是基本的编程规范吧!
void Bisection()
void linearinterpolation()
void linearinterpolation()
void Secant()
void Directiteration1()
void Directiteration2()
void Directiteration3()
5.程序的高效性(低时间消耗+低空间占用)
对于时间,由于算法固定了,就是这几个经典的,很难再有所减少了,但是编程习惯、代码风格却能在一定程度上减少空间的过度损耗,在算法中要存储每次迭代的x[n]和误差e[n],(当然可以不存储直接输出,但为了输出形式更美观我们组存储了),我们是设了全局变量,因为7个算法的实现都得用到它。另外,我们没用递归调用函数,而是在1个函数中用循环语句实现迭代,由于算法比较简单,代码量较少,用循环语句并未降低程序的可读性,还相应的减少了递归调用带来的额外计算机资源占用。我们知道,递归的话函数调用是有开销的,而且递归的次数受堆栈大小的限制,递归太深系统可能支撑不了。
下面是我们的一个Bisection函数,我们就是回避递归,转而用循环语句。不影响程序的可读性又有效地释放了空间,何乐而不为呢?
void Bisection()
{
cout<<"******************result of Bisection**********************"<<endl;
cout<<"初始化x1=0,x2=PI/2"<<endl;
x1=0;x2=PI/2;
int i=0;
while(1)
{
x[i]=(x1+x2)/2;
e[i]=x[i]-PI/3;
if(f(x1)*f(x[i])>0)
x1=x[i];
else if(f(x2)*f(x[i])>0)
x2=x[i];
else
x1=x2=x[i];
i++;
if(16==i)break;
}
print1(i);
}
6.我们在写的过程中,尝试着用字符串的形式输入函数,再通过字符串的比较由计算机自动识别函数,如cos、sin等C++标准模板库中有的函数就在直接包含头文件调用,成功是实现了部分,只是在和其他小组商量了之后才统一用书上的函数,再者是有的函数在某些点处不收敛,用NowtonRaphson();Secant();Directiteration3()可能不能快速地得到解。
六、总结
在知识层面感觉没什么好总结的,留给其他组吧。在这次以小组为单位的学习过程中,我感觉最大的收获就是大家都积极参与,无论基本功如何,我们组的每个成员都参与其中,都有事可做,并主动的做。
#include<iostream>
#include<cmath>
#include<iomanip>
using namespace std;
#define PI 3.141592653
double x0,x1,x2,x[20];
double e[20];
double f(double x)
{
return cos(x)-0.5;
}
void print1(int n)
{
int j;
cout<<"Iteration Error e[n+1]/e[n]"<<endl;
for(j=0;j<n;j++)
if(j<n-1)cout<<setw(8)<<j<<" "<<setw(20)<<x[j]-PI/3<<" "<<setw(20)<<e[j+1]/e[j]<<endl;
else cout<<setw(8)<<j<<" "<<setw(20)<<x[j]-PI/3<<" "<<endl<<endl;
}
void print2(int n)
{
int j;
cout<<"Iteration Error e[n+1]/e[n] e[n+1]/e[n]^2"<<endl;
for(j=0;j<n;j++)
if(j<n-1)cout<<setw(8)<<j<<" "<<setw(20)<<x[j]-PI/3<<" "<<setw(20)<<e[j+1]/e[j]<<" "<<setw(20)<<e[j+1]/pow(e[j],2)<<endl;
else cout<<setw(8)<<j<<" "<<setw(20)<<x[j]-PI/3<<" "<<endl<<endl;
}
void print3(int n)
{
int j;
cout<<"Iteration Error e[n+1]/e[n] |e[n+1]|/|e[n]|^1.618"<<endl;
for(j=0;j<n;j++)
if(j<n-1)cout<<setw(8)<<j<<" "<<setw(20)<<x[j]-PI/3<<" "<<setw(20)<<e[j+1]/e[j]<<" "<<setw(20)<<fabs(e[j+1])/pow(fabs(e[j]),1.618)<<endl;
else cout<<setw(8)<<j<<" "<<setw(20)<<x[j]-PI/3<<" "<<endl<<endl;
}
void Bisection()
{
cout<<"******************result of Bisection**********************"<<endl;
cout<<"初始化x1=0,x2=PI/2"<<endl;
x1=0;x2=PI/2;
int i=0;
while(1)
{
x[i]=(x1+x2)/2;
e[i]=x[i]-PI/3;
if(f(x1)*f(x[i])>0)
x1=x[i];
else if(f(x2)*f(x[i])>0)
x2=x[i];
else
x1=x2=x[i];
i++;
if(16==i)break;
}
print1(i);
}
void linearinterpolation()
{
cout<<"*****************result of linearinterpolation******************"<<endl;
cout<<"初始化x1=0,x2=PI/2"<<endl;
x1=0;x2=PI/2;
int i=0,j;
while(1)
{
x[i]=(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1));
e[i]=x[i]-PI/3;
if(f(x1)*f(x[i])>0)
x1=x[i];
else if(f(x2)*f(x[i])>0)
x2=x[i];
else
x1=x2=x[i];
i++;
if(16==i)break;
}
print1(i);
}
void NowtonRaphson()
{
int i=0;
cout<<"**********************************result of Nowton Raphson*********************************"<<endl;
cout<<"初始化x0=PI/2"<<endl;
x0=PI/2;
while(1)
{
x0=x0-(f(x0)/(-sin(x0)));
x[i]=x0;
e[i]=x0-PI/3;
i++;
if(5==i)break;
}
print2(i);
}
void Secant()
{
cout<<"***********************result of Secant**********************"<<endl;
cout<<"初始化x0=0,x1=PI/2"<<endl;
x0=0;x1=PI/2;
int i=0;
while(1)
{
if(0==i%2)
{
x0=x1-(x1-x0)*f(x1)/(f(x1)-f(x0));
x[i]=x0;
}
else
{
x1=x0-(x0-x1)*f(x0)/(f(x0)-f(x1));
x[i]=x1;
}
e[i]=x[i]-PI/3;
i++;
if(i==7)break;
}
print3(i);
}
double g(double x)
{
return x+cos(x)-0.5;
}
void Directiteration1()
{
cout<<"*********************result of Direct iteration1****************"<<endl;
cout<<"初始化x0=0"<<endl;
x0=0;
int i=0;
x0=g(x0);
e[i]=x0-PI/3;
while(1)
{
x[i]=x0;
x0=g(x0);
i++;
e[i]=x0-PI/3;
if(i==16)break;
}
print1(i);
}
double g2(double x)
{
return 2*x*cos(x);
}
void Directiteration2()
{
cout<<"*******************result of Direct iteration2******************"<<endl;
cout<<"初始化x0=PI/3+0.0635232"<<endl;
x0=PI/3+0.0635232;
int i=0;
while(1)
{
x[i]=x0;
e[i]=x0-PI/3;
x0=g2(x0);
i++;
if(i==16)break;
}
print1(i);
}
double h(double x)
{
return 4*x*(x-PI)/pow(PI,2);
}
void Directiteration3()
{
cout<<"*******************result of Direct iteration3******************"<<endl;
cout<<"初始化x0=PI/2"<<endl;
x0=PI/2;
int i=0;
while(1)
{
x0=x0-f(x0)/h(x0);
e[i]=x0-PI/3;
x[i]=x0;
i++;
if(i==10)break;
}
print1(i);
}
int main()
{
int choice;
int tag;
while(1)
{
tag=0;
do{
if(tag!=0)
cout<<"Your choice is illegal!Please input again!"<<endl;
tag++;
cout<<" =============Menu============="<<endl;
cout<<" || 1.Bisection|| ||"<<endl;
cout<<" || 2.linearinterpolation ||"<<endl;
cout<<" || 3.Nowton Raphson ||"<<endl;
cout<<" || 4.Secant ||"<<endl;
cout<<" || 5.Direct iteration1 ||"<<endl;
cout<<" || 6.Direct iteration2 ||"<<endl;
cout<<" || 7.Direct iteration3 ||"<<endl;
cout<<" || 8.all of seven ||"<<endl;
cout<<" || 9.return ||"<<endl;
cout<<" ==============================="<<endl;
cout<<" Please input your choice[ ]\b\b";
cin>>choice;
cout<<endl<<endl;
}while(choice<1||choice>9);
if(1==choice||8==choice)Bisection();
if(2==choice||8==choice)linearinterpolation();
if(3==choice||8==choice)NowtonRaphson();
if(4==choice||8==choice)Secant();
if(5==choice||8==choice)Directiteration1();
if(6==choice||8==choice)Directiteration2();
if(7==choice||8==choice)Directiteration3();
if(9==choice)exit(0);
}
return 0;
}
以上代码是大二写的,有很多错误及不合理之处,望斧正!