数值计算课程设计

1、 经典四阶龙格库塔法解一阶微分方程
1.1、 算法说明
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。

由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。

该算法是构建在数学支持的基础之上的。

4阶龙格-库塔方法(RK4)可模拟N=4的泰勒方法的精度。

这种算法可以描述为，
自初始点00(,)t y 开始，利用下面的计算方法生成近似序列
(1-1)
1.2、 经典四阶龙格库塔法解一阶微分方程算法流程图
图1-1 经典四阶龙格库塔法解一阶微分方程算法流程图
1.3、经典四阶龙格库塔法解一阶微分方程程序调试
图1-2 经典四阶龙格库塔法解一阶微分方程程序调试
1.4、经典四阶龙格库塔法解一阶微分方程代码
#include <iostream>
#include <iomanip>
using namespace std;
//f为函数的入口地址，x0、y0为初值，xn为所求点，step为计算次数
double Runge_Kuta( double (*f)(double x, double y), double x0, double y0, double xn, int step )
{ double k1,k2,k3,k4,result;
double h=(xn-x0)/step;
if(step<=0)
return(y0);
if(step==1)
{ k1=f(x0,y0);
k2=f(x0+h/2, y0+h*k1/2);
k3=f(x0+h/2, y0+h*k2/2);
k4=f(x0+h, y0+h*k3);
result=y0+h*(k1+2*k2+2*k3+k4)/6;}
else
{ double x1,y1;
x1=xn-h;
y1=Runge_Kuta(f, x0, y0, xn-h,step-1); k1=f(x1,y1);
k2=f(x1+h/2, y1+h*k1/2);
k3=f(x1+h/2, y1+h*k2/2);
k4=f(x1+h, y1+h*k3);
result=y1+h*(k1+2*k2+2*k3+k4)/6;}
return(result);}
int main()
{double f(double x, double y);
double x0,y0;
double a,b;
// int step;
cout<<"请输入初值x0,y0:";
cin>>x0>>y0;
cout<<"请输入区间：";
cin>>a>>b;
//double x0=0,y0=1;
double x,y,step;
int i;
cout<<"请输入步长：";
cin>>step;
//step=0.1;
cout.precision(10);
for(i=0;i<=(b-a)/step;i++)
{x=x0+i*step;
cout<<setw(8)<<x<<setw(18)<<Runge_Kuta(f,x0,y0,x,i)<<endl;} return(0); }
double f(double x, double y)
{ double r;
r=(x-y)/2;
return(r);
2、高斯列主元法解线性方程组
2.1、算法说明
首先将线性方程组做成增光矩阵，对增广矩阵进行行变换。

a，在第i列中，第i行及以下的元素选取绝对值最大的元素，将该对第元素ii
a消去第i行以下的元素最大的行与第i行交换，然后采用高斯消元法将新得到的ii
a。

从而得到上三角矩阵。

元素。

一次进行直到nn
再对得到的上三角矩阵进行回代操作，即可以得到方程组的解。

2.2、高斯列主元算法流程图
2.3、高斯列主元程序调试
对所编写的高斯列主元程序进行编译和链接，然后执行得如下所示的窗口，我们按命令输入增广矩阵的行数为3，输入3行5列的增广矩阵,运行界面为：
图2-2 高斯列主元程序调试
2.4、高斯列主元算法代码
#include<iostream>
#include<cmath>
using namespace std;
void load();
const N=20;
float a[N][N];
int m;
int main()
{ int i,j;
int c,k,n,p,r;
float x[N],l[N][N],s,d;
cout<<"下面请输入未知数的个数m=";
cin>>m;
cout<<endl;
cout<<"请按顺序输入增广矩阵a:"<<endl;
load();
for(i=0;i<m;i++)
{ for(j=i;j<m;j++)
c=(fabs(a[j][i])>fabs(a[i][i]))?j:i; /*找列最大元素*/ for(n=0;n<m+1;n++)
{s=a[i][n]; a[i][n]=a[c][n]; a[c][n]=s;} /*将列最大数防在对角线上*/
for(p=0;p<m+1;p++)
cout<<a[i][p]<<"\t";
cout<<endl;
for(k=i+1;k<m;k++)
{ l[k][i]=a[k][i]/a[i][i];
for(r=i;r<m+1;r++) /*化成三角阵*/
a[k][r]=a[k][r]-l[k][i]*a[i][r]; } }
x[m-1]=a[m-1][m]/a[m-1][m-1];
for(i=m-2;i>=0;i--)
{ d=0;
for(j=i+1;j<m;j++)
d=d+a[i][j]*x[j];
x[i]=(a[i][m]-d)/a[i][i];}
cout<<"该方程组的解为:"<<endl;
for(i=0;i<m;i++)
cout<<"x["<<i<<"]="<<x[i]<<"\t";
//system("pause"); return 0; }
void load()
{int i,j;
for(i=0;i<m;i++)
for(j=0;j<m+1;j++)
cin>>a[i][j];}
1122(,)(,)()(,)(,)k k k k k k k k k f p q f p q x y J P f p q f p q x y ∂∂⎡⎤
⎢⎥∂∂⎢⎥=∂∂⎢⎥
⎢⎥∂∂⎣⎦
12(,)()(,)k k k k k f p q F P f p q ⎡⎤=⎢⎥
⎣⎦
1k k P P P
+=+∆3、 牛顿法解非线性方程组
3.1、 算法说明
设k P 已知。

第1步：计算函数
（3-1）
第2步：计算雅可比矩阵
（3-2）
第3步：求线性方程组 ()()k k J P P F P ∆=-的解P ∆。

第4步：计算下一点
（3-3）
重复上述过程。

3.2、牛顿法解非线性方程组算法流程图
图3-1 算法流程图
222
20.50440
x x y x y ⎧--+=⎨+-=⎩3.3、 牛顿法解非线性方程组算法程序调试
图3-2牛顿法解非线性方程组算法程序调试
应用本程序解方程组， 初始近似值x0,y0分别为2.00和0.25，经过3次迭代求出X(1)=1.900691和X(2)=0.311213。

图3-2牛顿法解非线性方程组算法程序运行结果
3.4、牛顿法解非线性方程组算法程序代码
#include<iostream>
#include<cmath>
#define N 2 // 非线性方程组中方程个数、未知量个数
#define epsilon 0.0001 // 差向量1范数的上限
#define max 10 //最大迭代次数
using namespace std;
const int N2=2*N;
int main()
{
void ff(float xx[N],float yy[N]);
void ffjacobian(float xx[N],float yy[N][N]);
void inv_jacobian(float yy[N][N],float inv[N][N]);
void newdim(float x0[N], float inv[N][N],float y0[N],float x1[N]);
float
x0[N]={2.0,0.25},y0[N],jacobian[N][N],invjacobian[N][N],x1[N],errornorm; int i,iter=0;
cout<<"初始解向量："<<endl;
for (i=0;i<N;i++)
cout<<x0[i]<<" ";
cout<<endl;
do
{ iter=iter+1;
cout<<"第 "<<iter<<" 次迭代开始："<<endl;
//jis uan xiang liang han shu zhi---yin bian liang xiang liang y0 ff(x0,y0);
//jis uan jacobian ju zhen jacobian
ffjacobian(x0,jacobian);
//jis uan jacobian ju zhen de ni juzhen invjacobian
inv_jacobian(jacobian,invjacobian);
//you jie xiang liang x0 ji suan jie xiang liang x1
newdim(x0, invjacobian,y0,x1);
//ji suan cha xiang liang de 1 fan shu
errornorm=0;
for (i=0;i<N;i++)
errornorm=errornorm+fabs(x1[i]-x0[i]);
if (errornorm<epsilon) break;
for (i=0;i<N;i++)
x0[i]=x1[i];
} while (iter<max);
return 0;
}
void ff(float xx[N],float yy[N])
{ float x,y;
int i;
x=xx[0];
y=xx[1];
//非线性方程组
yy[0]=x*x-2*x-y+0.5;
yy[1]=x*x+4*y*y-4;
cout<<"因变量向量："<<endl;
for( i=0;i<N;i++)
cout<<yy[i]<<" ";
cout<<endl;
cout<<endl;
}
void ffjacobian(float xx[N],float yy[N][N]) {
float x,y;
int i,j;
x=xx[0];
y=xx[1];
yy[0][0]=2*x-2;
yy[0][1]=-1;
yy[1][0]=2*x;
yy[1][1]=8*y;
cout<<"雅克比矩阵： "<<endl;
for( i=0;i<N;i++)
{for(j=0;j<N;j++)
cout<<yy[i][j]<<" ";
cout<<endl;
}
cout<<endl;
}
void inv_jacobian(float yy[N][N],float inv[N][N])
{float aug[N][N2],L;
int i,j,k;
cout<<"计算雅克比矩阵的逆: "<<endl;
for (i=0;i<N;i++)
{for(j=0;j<N;j++)
aug[i][j]=yy[i][j];
for(j=N;j<N2;j++)
if(j==i+N) aug[i][j]=1;
else aug[i][j]=0;
}
for (i=0;i<N;i++)
{for(j=0;j<N2;j++)
cout<<aug[i][j]<<" ";
cout<<endl;
}
cout<<endl;
for (i=0;i<N;i++)
{
for (k=i+1;k<N;k++)
{L=-aug[k][i]/aug[i][i];
for(j=i;j<N2;j++)
aug[k][j]=aug[k][j]+L*aug[i][j];
}
}
for (i=0;i<N;i++)
{for(j=0;j<N2;j++)
cout<<aug[i][j]<<" ";
cout<<endl;
}
cout<<endl;
for (i=N-1;i>0;i--)
{ for (k=i-1;k>=0;k--)
{L=-aug[k][i] /aug[i][i];
for(j=N2-1;j>=0;j--)
aug[k][j]=aug[k][j]+L*aug[i][j];
}
}
for (i=0;i<N;i++)
{for(j=0;j<N2;j++)
cout<<aug[i][j]<<" ";
cout<<endl;
}
cout<<endl;
for (i=N-1;i>=0;i--)
for(j=N2-1;j>=0;j--)
aug[i][j]=aug[i][j]/aug[i][i];
for (i=0;i<N;i++)
{for(j=0;j<N2;j++)
cout<<aug[i][j]<<" ";
cout<<endl;
for(j=N;j<N2;j++)
inv[i][j-N]=aug[i][j];
}
cout<<endl;
cout<<"雅克比矩阵的逆： "<<endl;
for (i=0;i<N;i++)
{ for(j=0;j<N;j++)
cout<<inv[i][j]<<" ";
cout<<endl;
}
cout<<endl;
}
void newdim(float x0[N], float inv[N][N],float y0[N],float x1[N]) {
int i,j;
//计算非线性方程组的近似解向量
float sum=0;
for(i=0;i<N;i++)
{ sum=0;
for(j=0;j<N;j++)
sum=sum+inv[i][j]*y0[j];
x1[i]=x0[i]-sum;
}
cout<<"近似解向量："<<endl;
for (i=0;i<N;i++)
cout<<x1[i]<<" "<<endl; }
4、 龙贝格求积分算法
4.1、 算法说明 生成J K ≥的逼近表(,)R J K ，并以(1,1)R J J ++为最终解来逼近积分
()(,)b
a
f x dx R J J ≈⎰
（4-1）
逼近(,)R J K 存在于一个特别的下三角矩阵中，第0列元素(,0)R J 用基于2J 个[a,b]子区间的连续梯形方法计算，然后利用龙贝格公式计算(,)R J K 。

当1K J ≤≤时，第J 行的元素为
(,1)(1,1)
(,)(,1)41
K
R J K R J K R J K R J K ----=-+
- （4-2） 当|(,)(1,1)|R J J R J J tol -++<时，程序在第(1)J +行结束。

4.2、 龙贝格求积分算法流程图
图4-1 算法流程图
4.3、 龙贝格求积分算法程序调试
我们以求解积分())sin(2x x f =，精度为0.0001，最高迭代10次为例，对所编写的龙贝格求积分算法程序进行编译和链接，经执行后得如下所示的窗口
图4-2龙贝格求积分算法程序调试
说明：应用Romberg 算法求())sin(2x x f =在[]21区间上的精度为0.0001的积分为0.494508。

4.4、 龙贝格求积分算法代码 #include<iostream> #include<cmath> using namespace std;
#define f(x) sin(x*x) //举例函数 #define epsilon 0.0001 //精度
#define MAXREPT 10 //迭代次数,到最后仍达不到精度要求,则输出T(m=10). double Romberg(double aa, double bb) { //aa,bb 积分上下限
int m, n;//m 控制迭代次数, 而n 控制复化梯形积分的分点数. n=2^m double h, x; double s, q;
double ep; //精度要求
double *y = new double[MAXREPT];//为节省空间,只需一维数组
//每次循环依次存储Romberg 计算表的每行元素,以供计算下一行,算完后更新 double p;//p 总是指示待计算元素的前一个元素(同一行)
//迭代初值
h = bb - aa;
y[0] = h*(f(aa) + f(bb))/2.0;
m = 1;
n = 1;
ep = epsilon + 1.0;
//迭代计算
while ((ep >= epsilon) && (m < MAXREPT))
{
//复化积分公式求T2n(Romberg计算表中的第一列),n初始为1,以后倍增 p = 0.0;
for (int i=0; i<n; i++)//求Hn
{
x = aa + (i+0.5)*h;
p = p + f(x);
}
p = (y[0] + h*p)/2.0;//求T2n = 1/2(Tn+Hn),用p指示
//求第m行元素,根据Romberg计算表本行的前一个元素(p指示),
//和上一行左上角元素(y[k-1]指示)求得.
s = 1.0;
for (int k=1; k<=m; k++)
{
s = 4.0*s;
q = (s*p - y[k-1])/(s - 1.0);
y[k-1] = p;
p = q;
}
p = fabs(q - y[m-1]);
m = m + 1;
y[m-1] = q;
n = n + n; h = h/2.0;
}
return (q);
}
int main()
{
double a,b;
cout<<"Romberg积分,请输入积分范围a,b:"<<endl; cin>>a>>b;
cout<<"积分结果:"<<Romberg(a, b)<<endl;
return 0;
}
5、 三次样条插值算法
5.1 三次样条插值算法说明表
表5-1 三次样条插值算法说明表
策略描述
包含0m 和N m 的方程
(i)三次紧压样条，确定0()S x '，()n S x '（如果导数已知，这是“最佳选择”）
100003
(())2
m m d S x h '=
-- 1
11
3(())2
N N N N N m m S x d h ---'=
--
(ii)natural 三次样条（一条“松弛曲线”）
00m =，0N m =
(iii)外挂()S x ''到端点
021011
()h m m m m h -=-
11212
()
N N N N N N h m m m m h ------=-
(iv) ()S x ''是靠近端点的常量 01m m =，1N N m m -=
(v)在每个端点处指定()S x ''
00()m S x ''=，()N N m S x ''=
5.2、 三次样条插值算法（压紧样条）程序调试
我们将所编写的程序三次样条插值算法（压紧样条）程序进行调试
图5-1三次样条插值算法（压紧样条）程序输入界面、运行结果
图5-2三次样条插值算法程序运行结果（a ） 图5-2三次样条插值算法程序运行结果 （b ）
运行结果分析：
()
()()()
()
()
⎪
⎪
⎩
⎪⎪
⎨
⎧
≤
≤
-
+
-
-
-
+
-
≤
≤
-
+
-
-
-
-
-
-
≤
≤
-
+
-
-
-
+
-
=
3
2
)2
(
44
.1
3
62
.2
)2
(
06
.0
)3
(
62
.0
2
1
)1
(
62
.2
)2
(
08
.0
)1
(
62
.0
2
42
.0
1
)0
(
08
.0
1
06
.0
42
.0
1
06
.0
3
3
3
3
3
3
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
f（5-1）
作图程序（Matlab）：
x1=0:0.01:1;
y1=0.06*(x1 - 1).^3 + 0.42*(x1 - 0).^3 - 0.06*(x1 - 1) + 0.08*(x1 - 0);
x2=1:0.01:2;
y2=-0.42*(x2 - 2).^3 - 0.62*(x2 - 1).^3 - 0.08*(x2 - 2) + 2.62*(x2 - 1);
x3=2:0.01:3;
y3=0.62*(x3 - 3).^3 + 0.06*(x3 - 2).^3 - 2.62*(x3 - 3) + 1.44*(x3 - 2);
X=[0 1 2 3];
Y=[0 0.5 2 1.5];
plot(x1,y1,x2,y2,x3,y3,X,Y,'*')
gtext('S1')
gtext('S2')
gtext('S3')
图形为：
图5-3 三次样条插值算法Matlab 作图分析
5.3、三次样条插值算法（压紧样条）代码
#include<iostream>
#include<iomanip>
using namespace std;
const int max = 50;
float x[max], y[max], h[max];
float c[max], a[max], fxym[max];
float f(int x1, int x2, int x3)
{
float a = (y[x3] - y[x2]) / (x[x3] - x[x2]);
float b = (y[x2] - y[x1]) / (x[x2] - x[x1]);
return (a - b)/(x[x3] - x[x1]);
} //求差分
void cal_m(int n)
{ //用追赶法求解出弯矩向量M……
float B[max];
B[0] = c[0] / 2;
for(int i = 1; i < n; i++)
B[i] = c[i] / (2 - a[i]*B[i-1]);
fxym[0] = fxym[0] / 2;
for(i = 1; i <= n; i++)
fxym[i] = (fxym[i] - a[i]*fxym[i-1]) / (2 - a[i]*B[i-1]);
for(i = n-1; i >= 0; i--)
fxym[i] = fxym[i] - B[i]*fxym[i+1];
}
void printout(int n);
int main()
{
int n,i; char ch;
do{ cout<<"输入x的最大下标:";
cin>>n;
for(i = 0; i <= n; i++)
{
cout<<"Please put in X"<<i<<':';
cin>>x[i];
cout<<"Please put in Y"<<i<<':';
cin>>y[i];
}
for(i = 0; i < n; i++) //求步长
h[i] = x[i+1] - x[i];
cout<<"Please 输入边界条件\n 1 :已知两端的一阶导数\n 2 :两端的二阶导数已知\n 默认:自然边界条件\n";
int t;
float f0, f1;
cin>>t;
switch(t)
{
case 1:cout<<"输入 Y0\' Y"<<n<<"\'\n";
cin>>f0>>f1;
c[0] = 1; a[n] = 1;
fxym[0] = 6*((y[1] - y[0]) / (x[1] - x[0]) - f0) / h[0];
fxym[n] = 6*(f1 - (y[n] - y[n-1]) / (x[n] - x[n-1])) / h[n-1]; break;
case 2:cout<<"输入 Y0\" Y"<<n<<"\"\n";
cin>>f0>>f1;
c[0] = a[n] = 0;
fxym[0] = 2*f0; fxym[n] = 2*f1;
break;
default:cout<<"不可用\n";//待定
};
for(i = 1; i < n; i++)
fxym[i] = 6 * f(i-1, i, i+1);
for(i = 1; i < n; i++)
{
a[i] = h[i-1] / (h[i] + h[i-1]);
c[i] = 1 - a[i];
}
a[n] = h[n-1] / (h[n-1] + h[n]);
cal_m(n);
cout<<"\n输出三次样条插值函数：\n";
printout(n);
cout<<"Do you have another try ? y/n :";
cin>>ch;
}while(ch == 'y' || ch == 'Y');
return 0;}
void printout(int n)
{cout<<setprecision(6);
for(int i = 0; i < n; i++)
{cout<<i+1<<": ["<<x[i]<<" , "<<x[i+1]<<"]\n"<<"\t";
cout<<"S"<<i+1<<"=";
float t = fxym[i]/(6*h[i]);
if(t > 0) cout<<-t<<"*(x - "<<x[i+1]<<")^3";
else cout<<-t<<"*(x - "<<x[i+1]<<")^3";
t = fxym[i+1]/(6*h[i]);
if(t > 0) cout<<" + "<<t<<"*(x - "<<x[i]<<")^3"; else cout<<" - "<<t<<"*(x - "<<x[i]<<")^3";
cout<<"\n\t";
t = (y[i] - fxym[i]*h[i]*h[i]/6)/h[i];
if(t > 0) cout<<"- "<<t<<"*(x - "<<x[i+1]<<")"; else cout<<"- "<<-t<<"*(x - "<<x[i+1]<<")";
t = (y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i];
if(t > 0) cout<<" + "<<t<<"*(x - "<<x[i]<<")"; else cout<<" - "<<-t<<"*(x - "<<x[i]<<")";
cout<<endl; }
cout<<endl;
}
432211113211112111N N N N
k k k k k
k k k k N N N N
k k k k k k k k k N N N
k k k
k k k x A x B x C y x x A x B x C y x x A x B NC y ===========⎛⎫⎛⎫⎛⎫++= ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭
⎛⎫⎛⎫⎛⎫++= ⎪ ⎪
⎪⎝⎭⎝⎭⎝⎭⎛⎫⎛⎫
++= ⎪ ⎪⎝⎭⎝⎭
∑∑∑∑∑∑∑∑∑∑∑6、M 次多项式曲线拟合
6.1、算法说明
设{}1(,)N
k k k x y =有N 个点，横坐标是确定的。

最小二乘抛物线的系数表示为
2()y f x Ax Bx C ==++ （6-1）
求解A ，B 和C 的线性方程组为
（6-2）
6.2、 M 次多项式曲线拟合算法流程图
6.3、 M 次多项式曲线拟合算法程序调试
我们按命令依次输入命令如下命令后，得程序执行结果如下
图6-1 算法流程图
图6-2 算法调试图
作图程序：图形为：
X=[-3 -1 1 3];
Y=[15 5 1 5];
x=-3:0.01:3;
y=0.875*x.^2-1.7*x+2.145;
plot(X,Y,'go',x,y,'r-')
gtext('拟合曲线')
图6-3 图形分析
6.4、 M次多项式曲线拟合算法代码
#include<iostream>
#include<cmath>
#define MAX 20
using namespace std;
//求解任意可逆矩阵的逆，X为待求解矩阵，E为全零矩阵，非单位矩阵，也可以是单位矩阵
void inv(double X[MAX][MAX],int n,double E[MAX][MAX])
{
int i,j,k;
double temp=0;
for(i=0;i<MAX;i++)
{
for(j=0;j<MAX;j++)
if(i==j)
E[i][j]=1;
}
for(i=0;i<n-1;i++)
{
temp=X[i][i];
for(j=0;j<n;j++)
{
X[i][j]=X[i][j]/temp;
E[i][j]=E[i][j]/temp;
}
for(k=0;k<n;k++)
{
if(k==i)
continue;
temp=-X[i][i]*X[k][i];
for(j=0;j<n;j++)
{
X[k][j]=X[k][j]+temp*X[i][j];
E[k][j]=E[k][j]+temp*E[i][j];
}
}
}
}
int main()
{int n,M,i,j,k;
double X[MAX]={0},Y[MAX]={0},F[MAX][MAX]={0},B[MAX]={0};
double A[MAX][MAX]={0},BF[MAX][MAX]={0},C[MAX]={0},E[MAX][MAX]={0};
cout<<"M次多项式曲线拟合\n请先输入待拟合的点的个数:";
cin>>n;
cout<<"\n请输入"<<n<<"个点的X坐标序列:\n";
for(i=0;i<n;i++)
cin>>X[i];
cout<<"\n请输入"<<n<<"个点的Y坐标序列:\n";
for(i=0;i<n;i++)
cin>>Y[i];
cout<<"\n请输入需要拟合的次数:";
cin>>M;
for(i=0;i<n;i++)
for(k=1;k<=M+1;k++)
F[i][k-1]=pow(X[i],k-1);
//求F的转置
for(i=0;i<n;i++)
{
for(j=0;j<M+1;j++)
{
BF[j][i]=F[i][j];
}
}
//计算其转置的BF与F的乘
for(i=0;i<M+1;i++)
for(j=0;j<M+1;j++)
for(k=0;k<n;k++)
A[i][j]+=BF[i][k]*F[k][j];
//计算F的转置BF与Y的乘
for(i=0;i<M+1;i++)
for(j=0;j<n;j++)
B[i]+=BF[i][j]*Y[j];
//调用inv函数求解矩阵A的逆矩阵E inv(A,n,E);
//计算A的逆BF与B的乘
for(i=0;i<M+1;i++)
for(j=0;j<n;j++)
C[i]+=E[i][j]*B[j];
cout<<"\n拟合后的"<<M<<"次多项式系数为，幂次由高到低：\n";
for(i=M;i>=0;i--)
cout<<C[i]<<"\t";
cout<<"\n拟合后的"<<M<<"次多项式为:\n";
cout<<"\nP(x)=";
for(i=M;i>=0;i--)
{
if(i==0)
cout<<"+"<<C[i];
else
cout<<"+"<<C[i]<<"*x^"<<i;
}
cout<<endl;
return 0;
}
7、 牛顿-拉弗森迭代解非线性方程
7.1、 牛顿-拉弗森迭代解非线性方程算法概要
使用初始近似值0P ，利用迭代111()
()
k k k k f P P P f P ---=-' 其中1,2,...k =计算函数()0f x =的根的
近似值。

7.2、 牛顿-拉弗森迭代解非线性方程算法程序调试
图7-1牛顿-拉弗森迭代解非线性方程算法程序调试
7.3、 牛顿-拉弗森迭代解非线性方程算法代码 #include<iostream> using namespace std; #include<cmath> double f(double x) {
return x*x-2*x-1
}
double df(double x) {
return 2*x-2;
}
int main()
{ int k=0,max1=0;
double p0,delta,epsilon,p1,err,relerr,y;
cout<<"牛顿拉弗森法解非线性方程f(x)=x^2-2x-1"<<endl;
cout<<"初始值为p0=4.8"<<endl;
p0=4.8,delta=0.0001,epsilon=0.0001,max1=100;
for(k=1;k<=max1;k++)
{
p1=p0-f(p0)/df(p0);
err=fabs(p1-p0);
relerr=2*err/(fabs(p1)+delta);
p0=p1;
y=f(p0);
if((err<delta)||(relerr<delta)||(fabs(y)<epsilon))
break;
}
cout<<"方程的根的近似值为:"<<p0<<endl;
cout<<"方程的根的误差估计为:"<<err<<endl;
cout<<"迭代次数为:"<<k<<endl;
cout<<"方程在p0点的函数值f(p0):"<<y<<endl; return 0;
}
8、 不动点法解非线性方程
8.1、 算法说明 先将()f x 改写成()x g x =
然后对()x g x =进行迭代，即1()k k x g x -= 其中1,2,...k =
然后判断
1||k k x x ε--<是否成立，成立则返回k x ，不成立就重复以上步骤
8.2、 不动点法解非线性方程算法程序调试
我们将编写好的不动点迭代法解非线性方程算法程序进行编译，链接和执行后得如下所示结果。

8.3、 不动点法解非线性方程算法代码
#include<iostream> #include<cmath> #include<iomanip> #define MAX 20 #define eps 1e-10 using namespace std; double g(double x) { return 1-x*x/4+x; } main()
{ double P[MAX]={0},err=0.0,relerr=0.0,tol=0.0,p=0.0,p0=0.0; int k=0,max1=0,i=0; cout<<"
不
动
点
法
解
非
线
性方程f(x)=1-x^2/2"<<endl;
图8-1 程序调试图
cout<<"方程在[0，1]上有解，初始值为p0=0"<<endl;
//初始化
P[0]=p0=0;
max1=100;
tol=0.001;
for(k=2;k<=max1;k++)
{ P[k-1]=g(P[k-2]);
err=fabs(P[k-1]-P[k-2]);
relerr=err/(fabs(P[k-1]+eps));
p=P[k-1];
if((err<tol)||(relerr<tol))
break; }
if(k==max1)
cout<<"迭代次数超过允许的最大迭代次数！"<<endl;
cout<<"不动点的近似值为:"<<p<<endl;
cout<<"程序迭代次数为:"<<k<<endl;
cout<<"近似值的误差为:"<<err<<endl;
cout<<"求解不动点近似值的序列:"<<endl;
for(i=0;i<k;i++)
{ cout<<setw(10)<<P[i];
}
cout<<endl;
return 0;
}
9. 拉格朗日插值
9.1、 算法说明
()()()0011N N+1,,,,...,,N P n n x y x y x y x 有个点的次数最高为的多项式（)的构造方法，它具有
()()()N ,,0
P =N
k N k N k k x y L x L x =∑的形式，其中是基于节点
()()()()()()()()()()0111N 0111()k k n k k k k k k k k n x x x x x x x x x x L x x x x x x x x x x x -+-+-----=-----，的拉格朗日系数多项式
对每个固定的k,拉格朗日系数多项式(),N k L x 具有性质为：
,1,()0,N k i k i
L x k i
=⎧=⎨≠⎩. （9-1）
9.2、 拉格朗日插值算法程序调试
首先编写好程序，然后编译链接，运行程序，按照程序提示依次输入点的个数、点数对应的x 值、y 值、纵坐标、再输入x 值，输入结果如下：
图9-1拉格朗日插值算法程序调试
图9-2拉格朗日插值算法程序计算结果
9.3、 拉格朗日插值算法程序代码： #include<iostream> using namespace std;
double func(double X,int k,double x[],int n); int main()
{
double Sn=0;
int n;
cout<<"请输入点的个数n:";
cin>>n;
double*x=(double*)malloc(n*sizeof(double));
double*y=(double*)malloc(n*sizeof(double));
double X;
int i;
for(i=0;i<n;i++)
{
cout<<"请输入x"<<i+1<<",y"<<i+1<<":"<<endl;
cin>>x[i]>>y[i];
}
cout<<"请输入x";
cin>>X;
for(i=0;i<n;i++)
{
Sn=Sn+func(X,i,x,n)*y[i];
}
cout<<"通过拉格朗日插值公式所得的结果："<<"当x="<<X<<"时，y="<<Sn<<endl; return 0;
}
double func(double X,int k,double x[],int n)
{
int i;
double Pn=1;
double p;
for(i=0;i<n;i++)
h
{
if(i==k) continue;
else p=(X-x[i])/(x[k]-x[i]);
Pn=Pn*p;
}
return Pn;
10、 雅克比迭代
10.1、 雅克比迭代的基本算法说明
（1），建立判定条件来判断雅可比迭代是否收敛，因此我们定义N N *矩阵的严格对角优势：∑≠=>N
k j j kj kk a a 1
，其中N k ,...,2,1=，可知可以收敛。

（2）问题重点方程可以表示成下面的形式： 5
*2158
*4214
7y x z z
x y z y x -+=
++=-+=
（10-1）
这样就可以提出下列雅可比迭代过程： 5
*2158
*4214
7111k
k k k
k k k
k k y x z z x y z y x -+=
++=-+=
+++ （10-2）
代入初始点（x0,y0,z0）进行迭代。

10.2、 雅克比迭代的程序运行
图10-1雅克比迭代的程序运行
10.3、 雅克比迭代的程序代码 #include <iostream>
#include <cmath>
using namespace std;
#define delta 0.0000001 #define n 3
#define max1 1000
main()
{
int i,j,k;
double err,norm,A[n][n],B[n],P[n],X[n];
cout<<"Please input matrix A("<<n<<"*"<<n<<")"<<endl; for (i=0;i<n;i++)
for (j=0;j<n;j++)
cin>>A[i][j];
cout<<"Please input matrix B("<<n<<"*1)"<<endl;
for (i=0;i<n;i++)
cin>>B[i];
cout<<"Please input matrix P("<<n<<"*1)"<<endl;
for (i=0;i<n;i++)
cin>>P[i];
for (k=0;k<max1;k++)
{
for (j=0;j<n;j++)
{
X[j]=B[j];
for (i=0;i<n;i++)
X[j]=X[j]-A[j][i]*P[i];
X[j]=(X[j]+A[j][j]*P[j])/A[j][j]; }
norm=0;
for (i=0;i<n;i++)
norm=norm+pow((X[i]-P[i]),2);
norm=pow(norm,0.5);
err=fabs(norm);
for (i=0;i<n;i++)
P[i]=X[i];
if (err<delta)
goto loop; }
loop:cout<<"Solve X="<<endl;
for (i=0;i<n;i++) cout<<X[i]<<endl;}
11、设计体会及今后的改进意见
课程设计是我们专业课程知识综合应用的实践训练，着是我们迈向社会，从事职业工作前一个必不少的过程。

”千里之行始于足下”，通过这次课程设计，我深深体会到这句千古名言的真正含义。

我今天认真的进行课程设计，学会脚踏实地迈开这一步，就是为明天能稳健地在社会大潮中奔跑打下坚实的基础。

通过这次“典型数值算法的C++语言程序设计”，使得我在多方面都有所提高。

通过这次课程设计，综合运用本专业所学课程的理论和C++语言编程的实际训练从而培养和提高学生独立工作能力，巩固与扩充了数值分析所学的内容，通过对编程的锻炼，对Matlab 和C++运行的环境有了更深入的了解。

在此感谢我们的刘海峰老师.，老师严谨细致、一丝不苟的作风一直是我工作、学习中的榜样；老师循循善诱的教导和不拘一格的思路给予我无尽的启迪；这次数值分析课程设计的每个实验细节和每个数据，都离不开老师您的细心指导。

同时感谢对我帮助过的同学们，谢谢你们对我的帮助和支持，让我感受到同学的友谊。

对于独立编程这个模块，自己还是缺乏训练，导致很多思路，算法无法用C++语言实现，因此自己还得“百尺竿头更进一步”，继续努力学好C++语言和MatlAB软件的使用。

继续努力做好编程的基础，为自己的以后工作打下扎实基础。

由于我的设计能力有限，在设计过程中难免出现错误，恳请老师们多多指教,我十分乐意接受你们的批评与指正，我将万分感谢。