数值分析源代码.docx
《数值分析源代码.docx》由会员分享，可在线阅读，更多相关《数值分析源代码.docx（13页珍藏版）》请在冰豆网上搜索。
数值分析源代码
一,三次样条插值
#include
#include
voidready（double*&A,double*&B,double*h,double*y,intm,intn,doublem0=0,doublemn=0）
{A=newdouble[n];B=newdouble[n];
if（m==2）{
A[0]=0;B[0]=2*m0;
A[n]=1;B[n]=2*mn;
}
else{A[0]=1;B[0]=3*（y[1]-y[0]）/h[0];
A[n]=0;B[n]=3*（y[n]-y[n-1]）/h[n-1];
}
for（inti=1;i{A[i]=h[i-1]/（h[i-1]+h[i]）;
B[i]=3*（（1-A[i]）*（y[i]-y[i-1]）/h[i-1]+A[i]*（y[i+1]-y[i]）/h[i]）;
}
}
voidcalculate1（double*&a,double*&b,double*A,double*B,intn）
{a=newdouble[n];b=newdouble[n];
a[0]=-A[0]/2;b[0]=B[0]/2;
for（inti=1;i<=n;i++）
{a[i]=-A[i]/（2+（1-A[i]）*a[i-1]）;
b[i]=（B[i]-（1-A[i]）*b[i-1]）/（2+（1-A[i]）*a[i-1]）;
}
}
voidcalculate2（double*&m,double*a,double*b,intn）
{m=newdouble[n+1];
m[n+1]=0;
for（inti=n;i>=0;i--）m[i]=m[i+1]*a[i]+b[i];
}
voidmain（）
{inte,t,w,v;
out:
cout<<"求三次样条插值函数:
"<"<cin>>e;
double*x,*y,*h,*A,*B,*a,*b,*m,s0,sn,u,s;
cout<<"请输入x的一系列值:
"<x=newdouble[e];
for（inti=0;i>x[i];
cout<<"请输入y的一系列值:
"<y=newdouble[e];
for（intf=0;f>y[f];
h=newdouble[e];
for（intd=0;dloop:
cout<<"边界条件选择:
1.提供了s\"（x1）=0,s\"（x2）=0;2.提供了s'（x0）=m0,s'（xn）=mn（m0,mn!
=0）"
<cin>>t;
if（t!
=1&&t!
=2）{cout<<"输入信息错误!
请重新输入!
"<if（t==1）ready（A,B,h,y,t,e）;
else{cout<<"请输入s'（x0）,s'（xn）的值:
"<cin>>s0>>sn;
ready（A,B,h,y,t,e,s0,sn）;}
calculate1（a,b,A,B,e）;
calculate2（m,a,b,e）;
in:
cout<<"请输入所需要计算的x的值:
（"<cin>>u;
if（ux[e-1]）{cout<<"输入数据错误!
"<for（intc=0;c{if（u>x[c]&&us=（1+2*（u-x[w]）/（x[w+1]-x[w]））*（（u-x[w+1]）*（u-x[w+1]）/（（x[w]-x[w+1]）*（x[w]-x[w+1]）））*y[w]+
（1+2*（u-x[w+1]）/（x[w]-x[w+1]））*（（u-x[w]）*（u-x[w]）/（（x[w+1]-x[w]）*（x[w+1]-x[w]）））*y[w+1]+
（u-x[w]）*（（u-x[w+1]）*（u-x[w+1]）/（（x[w]-x[w+1]）*（x[w]-x[w+1]）））*m[w]+
（u-x[w+1]）*（（u-x[w]）*（u-x[w]）/（（x[w+1]-x[w]）*（x[w+1]-x[w]）））*m[w+1];
cout<<"f（"<"<In:
cout<<"请选择:
1,继续进行新的操作计算2,退出不计算"<cin>>v;
if（v!
=1&&v!
=2）{cout<<"输入错误!
请重新输入!
"<if（v==1）{system（"cls"）;gotoout;}
}
运行结果:
根据课本52页的数据输入,得:
第9题:
f（78.3）=3.00304
第10题:
f（0.35）=0.591591
二,自动选取步长梯形法
#include
#include
#include
doublef（doublex）
{return2/（1+x*x）;}
voidmain（）
{doublea,b,e,h,T0,T1,S;
out:
cout<<"请输入区间上下界:
a和b"<cin>>a>>b;
cout<<"请输入给定的精度e:
"<cin>>e;
h=（b-a）/2;
T1=（f（a）+f（b））*h;
intn=1,t,v;
do
{T0=T1;
S=0;
for（intk=1;kS=S+f（a+（2*k-1）*h/n）;
T1=T0/2+S*h/n;
t=n;
if（fabs（T1-T0）>=3*e）n=2*n;
}while（t==n）;
cout<<"该积分的近似值为:
"<In:
cout<<"请选择:
1,继续进行新的操作计算2,退出不计算"<cin>>v;
if（v!
=1&&v!
=2）{cout<<"输入错误!
请重新输入!
"<if（v==1）{system（"cls"）;gotoout;}
}
运行结果:
根据课本97页第7题的数据输入,得结果为:
0.75
三,Romberg求积分法
#include
#include
#include
usingnamespacestd;
doublef（doublex）
{returnsqrt（x）;}
voidmain（）
{doublea,b,e,T[20][20],S=0,s;intk=1,v;
out:
cout<<"请输入区间上下界:
a和b"<cin>>a>>b;
cout<<"请输入给定的精度e:
"<cin>>e;
T[0][0]=（b-a）*（f（a）+f（b））/2;
do
{for（inti=1;double（i）<=pow（2,double（k））-1;i++）
{S=S+f（a+（2*i-1）*（b-a）/pow（2,double（k）））;}
T[0][k]=0.5*（T[0][k-1]+（b-a）*S/pow（2,double（k）-1））;
for（intm=1;m<=k;m++）
{T[m][k-m]=（T[m-1][k-m+1]*pow（4,double（m））-T[m-1][k-m]）/（pow（4,double（m））-1）;}
s=k;
if（fabs（T[k][0]-T[k-1][0]）>=e）k=k+1;
}while（s!
=k）;
cout<<"该积分的近似值为:
"<In:
cout<<"请选择:
1,继续进行新的操作计算2,退出不计算"<cin>>v;
if（v!
=1&&v!
=2）{cout<<"输入错误!
请重新输入!
"<if（v==1）{system（"cls"）;gotoout;}
}
运行结果:
根据课本97页第8题的数据输入（输入a=1,b=9,e=0.000001）,得结果为:
92.121
四,列主元高斯消去法
#include
#include
#include
usingnamespacestd;
voidmain（）
{intn,v;
out:
cout<<"请输入系数矩阵A的行列数n:
"<cin>>n;
doubleA[10][10],*b=newdouble[n],*x=newdouble[n],e,T;inti1;
cout<<"请输入系数矩阵A的值:
"<for（inti=1;i<=n;i++）
{
for（intj=1;j<=n;j++）
cin>>A[i][j];
}
cout<<"系数矩阵A为:
"<for（inti=1;i<=n;i++）
{
for（intj=1;j<=n;j++）
cout<cout<}
cout<<"请输入方程组的右端项b:
"<for（inti=1;i<=n;i++）cin>>b[i];
cout<<"右端项b为:
"<for（inti=1;i<=n;i++）cout<cout<cout<<"请输入精度值e:
"<cin>>e;
for（intk=1;k{doublemax=fabs（A[k][k]）;
for（inti=k;i<=n;i++）
{
if（fabs（A[i][k]）>max）{max=fabs（A[i][k]）;i1=i;}
}
T=max;
if（T"<if（i1!
=k）
{doubletemp1;
for（intk1=1;k1<=n;k1++）
{doubletemp;
temp=A[i1][k1];
A[i1][k1]=A[k][k1];
A[k][k1]=temp;
}
temp1=b[i1];
b[i1]=b[k];
b[k]=temp1;
}
for（inti=k+1;i<=n;i++）
{
T=A[i][k]/A[k][k];
b[i]=b[i]-T*b[k];
for（intj=k+1;j<=n;j++）
A[i][j]=A[i][j]-T*A[k][j];
}
}
if（fabs（A[n][n]）<=e）{cout<<"求解失败!
"<x[n]=b[n]/A[n][n];
for（inti=n-1;i>=1;i--）
{doublesum=0;
for（intj=i+1;j<=n;j++）sum+=A[i][j]*x[j];
x[i]=（b[i]-sum）/A[i][i];
}
for（inti=1;i<=n;i++）cout<cout<stop:
;
In:
cout<<"请选择:
1,继续进行新的操作计算2,退出不计算"<cin>>v;
if（v!
=1&&v!
=2）{cout<<"输入错误!
请重新输入!
"<if（v==1）{system（"cls"）;gotoout;}}
运行结果:
根据课本126页第1题的数据输入,
系数矩阵为235右端项为5精度为0.00001;
3476
1335
结果为-4
1
2
五,迭代法--Jacobi迭代与Seidel迭代
#include
#include
#include
usingnamespacestd;
voidJacobi（）;
voidSeidel（）;
voidmain（）
{intw,v;
out:
loop:
cout<<"请选择所要的迭代方式:
1,Jacobi迭代2,Seidel迭代"<cin>>w;
switch（w）
{
case1:
Jacobi（）;break;
case2:
Jacobi（）;break;
default:
{cout<<"输入数据错误!
请重新输入!
"<}
In:
cout<<"请选择:
1,继续进行新的操作计算2,退出不计算"<cin>>v;
if（v!
=1&&v!
=2）{cout<<"输入错误!
请重新输入!
"<if（v==1）{system（"cls"）;gotoout;}
}
voidJacobi（）
{intn,k;
cout<<"请输入系数矩阵A的行数n:
"<cin>>n;
doubleA[10][10],*b=newdouble[n],*x=newdouble[n],*y=newdouble[n],*g=newdouble[n],e,T;intM;
cout<<"请输入系数矩阵A的值:
"<for（inti=1;i<=n;i++）
{
for（intj=1;j<=n;j++）
cin>>A[i][j];
}
cout<<"系数矩阵A为:
"<for（inti=1;i<=n;i++）
{
for（intj=1;j<=n;j++）
cout<cout<}
cout<<"请输入方程组的右端项b:
"<for（inti=1;i<=n;i++）cin>>b[i];
cout<<"右端项b为:
"<for（inti=1;i<=n;i++）cout<cout<cout<<"请输入容许误差e:
"<cin>>e;
cout<<"请输入容许最大迭代次数M:
"<cin>>M;
cout<<"请输入初始向量Y:
"<for（inti=1;i<=n;i++）cin>>y[i];
k=1;
//--------形成迭代矩阵-------
for（inti=1;i<=n;i++）
{
if（fabs（A[i][i]）"<T=A[i][i];
for（intj=1;j<=n;j++）A[i][j]=-A[i][j]/T;
A[i][i]=0;
g[i]=b[i]/T;
}
//--------迭代--------
loop:
doublesum=0;
for（inti=1;i<=n;i++）
{
for（intj=1;j<=n&&j!
=i;j++）sum+=（A[i][j]*y[j]）;
x[i]=g[i]+sum;
}
doublefanshu=0;
for（inti=1;i<=n;i++）
{
fanshu+=fabs（x[i]-y[i]）;
}
if（fanshu{for（inti=1;i<=n;i++）cout<cout<}
elseif（kelsecout<<"求解失败!
"<stop:
;}
voidSeidel（）
{
intn,k;
cout<<"请输入系数矩阵A的行数n:
"<cin>>n;
doubleA[10][10],*b=newdouble[n],*x=newdouble[n],*y=newdouble[n],*g=newdouble[n],e,T;intM;
cout<<"请输入系数矩阵A的值:
"<for（inti=1;i<=n;i++）
{
for（intj=1;j<=n;j++）
cin>>A[i][j];
}
cout<<"系数矩阵A为:
"<for（inti=1;i<=n;i++）
{
for（intj=1;j<=n;j++）
cout<cout<}
cout<<"请输入方程组的右端项b:
"<for（inti=1;i<=n;i++）cin>>b[i];
cout<<"右端项b为:
"<for（inti=1;i<=n;i++）cout<cout<cout<<"请输入容许误差e:
"<cin>>e;
cout<<"请输入容许最大迭代次数M:
"<cin>>M;
cout<<"请输入初始向量Y:
"<for（inti=1;i<=n;i++）cin>>y[i];
k=1;
for（inti=1;i<=n;i++）x[i]=y[i];
//--------形成迭代矩阵-------
for（inti=1;i<=n;i++）
{
if（fabs（A[i][i]）"<T=A[i][i];
for（intj=1;j<=n;j++）A[i][j]=-A[i][j]/T;
A[i][i]=0;
g[i]=b[i]/T;
}
//--------迭代--------
loop:
doublesum=0;
for（inti=1;i<=n;i++）
{
for（intj=1;j<=n&&j!
=i;j++）sum+=（A[i][j]*x[j]）;
x[i]=g[i]+sum;
}
doublefanshu=0;
for（inti=1;i<=n;i++）
{
fanshu+=fabs（x[i]-y[i]）;
}
if（fanshu{for（inti=1;i<=n;i++）cout<cout<}
elseif（kelsecout<<"求解失败!
"<stop:
;
}
运行结果:
根据课本139页第1题的数据输入,
简单迭代:
x1=1.2,x2=1.4,x3=1.6,x4=0.79999
Seidel迭代:
x1=1.2,x2=1.4,x3=1.6,x4=0.79999
升级会员 