nonlinear.cpp

本人的计算方法课程作业

//nonlinear.cpp--二分法和割线法求非线性方程组的解
#include<iostream>
#include<cmath>
using namespace std;
const double Pi=3.141593;
const float ERROR=0.0001;            //设定误差
const int Nmax=100;                  //设定最大迭代次数
float func(float t,float alpha);
float realf(float x);
float Romberg(float t);
float SECANT1(float xx[],float e1,float e2,float x);
void BISECTION(float xt[],float e1,float e2,float x,float root[]);
void main()
{
	float x1,x2,zone[2];
    float x[2]={2,3};
	float err1=ERROR,err2=ERROR;
    BISECTION(x,err1,err2,x2,zone);
    x1=SECANT1(zone,err1,err2,x1);
    cout<<"用二分法确定一个包含解的小区间为: ["<<zone[0]<<","<<zone[1]<<"]"<<endl;
	cout<<"用割线法求得该方程在上述区间内的根为:"<<x1<<endl;
}
//func()定义被积函数
float func(float t,float alpha)
{
	float f;
	f=cos(t*cos(alpha));
	return f;
}
//realf()为非线性方程左端的函数
float realf(float x)
{
	float r;
    r=Romberg(x);
	r/=Pi;
	return r;
}
//SECANT1()用割线法求realf(x)=0的根，x0,x1是初始近似，xx[]为迭代区间
float SECANT1(float xx[],float e1,float e2,float x)
{
	int k;
	float x0,x1,f0,f1,f,t1,t2,s;
	x0=xx[0];
	x1=xx[1];
	f0=realf(x0);
	f1=realf(x1);
	if(fabs(f1)>fabs(f0))
	{
		t1=x0;x0=x1;x1=t1;
		t2=f0;f0=f1;f1=t2;
	}
	for(k=0;k<=Nmax;k++)
	{
		if(fabs(f1)<e2) 
			return x1;
		s=f1/f0;
		x=x1-((x0-x1)*s)/(1-s);
		f=realf(x);
		if(fabs(x1-x)<(e1*fabs(x)))
			return x;
		if(fabs(f)>fabs(f1))
		{x0=x;f0=f;}
		else
		{x0=x1;f0=f1;x1=x;f1=f;}
	}
	cout<<"Nmax次迭代后仍不收敛!"<<endl;
}
//BISECTION()用二分法确定realf(x)=0的根的所在区间，结果存放于数组root[]中，
//其中的参数与SECANT1()定义相同
void BISECTION(float xt[],float e1,float e2,float x,float root[])
{
	float x0,x1,f0,f1,f;
	x0=xt[0];
	x1=xt[1];
	f0=realf(x0);
	f1=realf(x1);
	if((f0*f1)>0||(fabs(f0)<e2)||(fabs(f1)<e2))
		goto stop;
A5:x=(x0+x1)/2.0;
	f=realf(x);
	if((fabs(x1-x0)<(e1*fabs(x1)))||(fabs(f)<e2))
		goto stop;
	if(f1*f<0) 
	{x0=x;f0=f;}
	else 
	{x1=x;f1=f;}
	goto A5;
stop: root[0]=x0;
	  root[1]=x1;
}
//Romberg()计算给定函数func()的积分，积分区间为[a,b],误差不超过ERROR
//结果返回积分值 
float Romberg(float t)
{
	float a=0,b=Pi,h,x,S;
	float S1,S2,T1,T2,C1,C2,R1,R2,R;
	int k=2;
	h=b-a;
	T1=(func(t,a)+func(t,b))*h/2;
A2:	S=0;
	x=a+h/2;
A3:	S=S+func(t,x);
	x=x+h;
	if(x<b)
	    goto A3;
	T2=(T1+h*S)/2;
	S2=T2+(T2-T1)/3;
	if(k>2)
	    goto A8;
A7:	k=k+1;
	h=h/2;
	T1=T2;
	S1=S2;
	goto A2;
A8:	C2=S2+(S2-S1)/15;
	if(k>3)
	    goto A11;
A10:	C1=C2;
	goto A7;
A11:	R2=C2+(C2-C1)/63;
	if(k>4)
	    goto A14;
A13:	R1=R2;
	goto A10;
A14:	if(fabs(R2-R1)>ERROR)
	     goto A13;
	R=R2;
	return R;
}