数值计算课设.docx
《数值计算课设.docx》由会员分享，可在线阅读，更多相关《数值计算课设.docx（48页珍藏版）》请在冰点文库上搜索。
数值计算课设
1.经典四阶龙格库塔法解一阶微分方程
1.1算法说明
龙格-库塔（Runge-Kutta）方法是一种在工程上应用广泛的高精度单步算法。
由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。
该算法是构建在数学支持的基础之上的。
4阶龙格-库塔方法（RK4）可模拟N=4的泰勒方法的精度。
这种算法可以描述为，自初始点
开始，利用下面的计算方法生成近似序列
1.2经典四阶龙格库塔法解一阶微分方程算法流程图
图1-1算法流程图
1.3经典四阶龙格库塔法解一阶微分方程程序调试
将编写好的代码放在VC6.0环境中编译，直接执行程序便可以得到求解微分方程
，并且
的结果。
如下图：
图1-2运行结果
然后将这些点进行插值或者拟合后就可以得到微分方程的解。
1.4经典四阶龙格库塔法解一阶微分方程程序代码：
#include
#include
usingnamespacestd;
doubleRunge_Kuta（double（*f）（doublex,doubley）,doublex0,doubley0,doublexn,longstep）
{doublek1,k2,k3,k4,result;
doubleh=（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
{
doublex1,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）;}
intmain（）
{
doublef（doublex,doubley）;
doublex0,y0;
doublea,b;
cout<<"请输入初值x0,y0:
";
cin>>x0>>y0;
cout<<"请输入区间：
";
cin>>a>>b;
//doublex0=0,y0=1;
doublex,y,step;
longi;
cout<<"请输入步长：
";
cin>>step;
//step=0.1;
cout.precision（10）;
for（i=0;i<=（b-a）/step;i++）
{
x=x0+i*step;
cout<return（0）;}
doublef（doublex,doubley）
{doubler;
r=（x-3*y）/5;
return（r）;
}
2.高斯列主元法解线性方程组
2.1算法说明
首先将线性方程组做成增光矩阵，对增广矩阵进行行变换。
对第元素
，在第i列中，第i行及以下的元素选取绝对值最大的元素，将该元素最大的行与第i行交换，然后采用高斯消元法将新得到的
消去第i行以下的元素，一次进行直到
，从而得到上三角矩阵。
再对得到的上三角矩阵进行回代操作，即可以得到方程组的解。
2.2高斯列主元算法流程图
开始
输入
（增广矩阵）
是
否
交换
中
两行
是
否
输出
结束
图2-1高斯列主元法解线性方程组算法流程图
2.3高斯列主元程序调试
对所编写的高斯列主元程序进行编译和链接，然后执行得如下所示的窗口，我们按命令输入增广矩阵的行数为3，输入3行4列的增广矩阵，输入时界面图如下，按回车键后，程序执行结果如下图所示：
图2-2运行结果图
2.4高斯列主元算法程序代码：
#include
#include
usingnamespacestd;
voidload（）;
constN=20;
floata[N][N];
intm;
intmain（）
{inti,j;
intc,k,n,p,r;
floatx[N],l[N][N],s,d;
cout<<"下面请输入未知数的个数m=";
cin>>m;
cout<cout<<"请按顺序输入增广矩阵a:
"<load（）;
for（i=0;i{for（j=i;jc=（fabs（a[j][i]）>fabs（a[i][i]））?
j:
i;
for（n=0;n{s=a[i][n];a[i][n]=a[c][n];a[c][n]=s;}
for（p=0;pcout<cout<for（k=i+1;k{l[k][i]=a[k][i]/a[i][i];
for（r=i;ra[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;jd=d+a[i][j]*x[j];
x[i]=（a[i][m]-d）/a[i][i];
}
cout<<"该方程组的解为:
"<for（i=0;icout<<"x["<}
voidload（）
{
inti,j;
for（i=0;ifor（j=0;jcin>>a[i][j];
}
3.牛顿法解非线性方程组
3.1算法说明
设
已知。
第1步：
计算函数
第2步：
计算雅可比矩阵
第3步：
求线性方程组
的解
。
第4步：
计算下一点
重复上述过程。
3.2牛顿法解非线性方程组算法流程图
图3-1算法流程图
3.3牛顿法解非线性方程组算法程序调试
以方程组
为例，初始近似值x0,y0分别为2.00和0.25，经过10次迭代求出X
（1）=4.30017和X
（2）=17.6736。
运行程序，可以得到结果如下：
图3-2运行结果图
3.4牛顿法解非线性方程组算法程序代码
#include
#include
#defineN2
#defineepsilon0.0001
#definemax10
usingnamespacestd;
constintN2=2*N;
intmain（）
{voidff（floatxx[N],floatyy[N]）;
voidffjacobian（floatxx[N],floatyy[N][N]）;
voidinv_jacobian（floatyy[N][N],floatinv[N][N]）;
voidnewdim（floatx0[N],floatinv[N][N],floaty0[N],floatx1[N]）;
floatx0[N]={2.0,0.25},y0[N],jacobian[N][N],invjacobian[N][N],x1[N],errornorm;
inti,iter=0;
cout<<"初始解向量：
"<for（i=0;icout<cout<do
{iter=iter+1;
cout<<"第"<"<ff（x0,y0）;
ffjacobian（x0,jacobian）;
inv_jacobian（jacobian,invjacobian）;
newdim（x0,invjacobian,y0,x1）;
errornorm=0;
for（i=0;ierrornorm=errornorm+fabs（x1[i]-x0[i]）;
if（errornormfor（i=0;ix0[i]=x1[i];
}while（iterreturn0;}
voidff（floatxx[N],floatyy[N]）
{floatx,y;
inti;
x=xx[0];
y=xx[1];
yy[0]=3*x*x+2*y-4;
yy[1]=2*x+2*y-3;
cout<<"因变量向量：
"<for（i=0;icout<cout<cout<}
voidffjacobian（floatxx[N],floatyy[N][N]）
{floatx,y;
inti,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<<"雅克比矩阵：
"<for（i=0;i{for（j=0;jcout<cout<cout<voidinv_jacobian（floatyy[N][N],floatinv[N][N]）
{floataug[N][N2],L;
inti,j,k;
cout<<"计算雅克比矩阵的逆:
"<for（i=0;i{for（j=0;jaug[i][j]=yy[i][j];
for（j=N;jif（j==i+N）aug[i][j]=1;
elseaug[i][j]=0;}
for（i=0;i{for（j=0;jcout<cout<cout<for（i=0;i{for（k=i+1;k{L=-aug[k][i]/aug[i][i];
for（j=i;jaug[k][j]=aug[k][j]+L*aug[i][j];}}
for（i=0;i{for（j=0;jcout<cout<cout<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{for（j=0;jcout<cout<cout<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{for（j=0;jcout<cout<for（j=N;jinv[i][j-N]=aug[i][j];}
cout<cout<<"雅克比矩阵的逆：
"<for（i=0;i{for（j=0;jcout<cout<cout<voidnewdim（floatx0[N],floatinv[N][N],floaty0[N],floatx1[N]）
{inti,j;
floatsum=0;
for（i=0;i{sum=0;
for（j=0;jsum=sum+inv[i][j]*y0[j];
x1[i]=x0[i]-sum;}
cout<<"近似解向量：
"<for（i=0;icout<4.龙贝格求积分算法
4.1算法说明
生成
的逼近表
，并以
为最终解来逼近积分
逼近
存在于一个特别的下三角矩阵中，第0列元素
用基于
个[a,b]子区间的连续梯形方法计算，然后利用龙贝格公式计算
。
当
时，第
行的元素为
当
时，程序在第
行结束。
4.2龙贝格求积分算法流程图
图4-1算法流程图
4.3龙贝格求积分法程序调试
我们以求解积分方程
，精度为0.0001，最高迭代10次为例，对所编写的龙贝格求积分算法程序进行编译和链接，经执行后得如下所示的窗口：
图4-2运行结果
4.4龙贝格求积分算法程序代码
#include
#include
#definef（x）sin（x*x）
#defineepsilon0.0001
#defineMAXREPT10
usingnamespacestd;
doubleRomberg（doubleaa,doublebb）
{
intm,n;
doubleh,x;
doubles,q;
doubleep;
double*y=newdouble[MAXREPT];
doublep;
h=bb-aa;
y[0]=h*（f（aa）+f（bb））/2.0;
m=1;
n=1;
ep=epsilon+1.0;
while（（ep>=epsilon）&&（m{
p=0.0;
for（inti=0;i{
x=aa+（i+0.5）*h;
p=p+f（x）;
}
p=（y[0]+h*p）/2.0;
s=1.0;
for（intk=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）;
}
intmain（）
{
doublea,b;
cout<<"Romberg积分,请输入积分范围a,b:
"<cin>>a>>b;
cout<<"积分结果:
"<return0;
}
5.三次样条插值算法
5.1算法说明
表5-1算法说明
策略描述
包含
和
的方程
（i）三次紧压样条，确定
，
（如果导数已知，这是“最佳选择”）
（ii）natural三次样条（一条“松弛曲线”）
，
（iii）外挂
到端点
（iv）
是靠近端点的常量
，
（v）在每个端点处指定
，
5.2三次样条插值算法（压紧样条）程序调试
求解三次紧压样条曲线，以过点（0，0.0），（1，1.5），（2，2.5）和（3，3.5），而且一阶导数边界条件S'（0）=0.5和S'
（2）=-5的三次压紧样条曲线。
将所编写的程序进行编译，链接和执行后得如下所示的结果：
图5-1运行结果图
5.3三次样条插值算法（压紧样条）代码
#include
#include
usingnamespacestd;
constintmax=50;
floatx[max],y[max],h[max];
floatc[max],a[max],fxym[max];
floatf（intx1,intx2,intx3）
{
floata=（y[x3]-y[x2]）/（x[x3]-x[x2]）;
floatb=（y[x2]-y[x1]）/（x[x2]-x[x1]）;
return（a-b）/（x[x3]-x[x1]）;
}
voidcal_m（intn）
{
floatB[max];
B[0]=c[0]/2;
for（inti=1;iB[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];
}
voidprintout（intn）;
intmain（）
{
intn,i;charch;
do{cout<<"输入x的最大下标:
";
cin>>n;
for（i=0;i<=n;i++）
{
cout<<"PleaseputinX"<';
cin>>x[i];
cout<<"PleaseputinY"<';
cin>>y[i];
}
for（i=0;ih[i]=x[i+1]-x[i];
cout<<"Please输入边界条件\n1:
已知两端的一阶导数\n2:
两端的二阶导数已知\n默认:
自然边界条件\n";
intt;
floatf0,f1;
cin>>t;
switch（t）
{
case1:
cout<<"输入Y0\'Y"<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;
case2:
cout<<"输入Y0\"Y"<cin>>f0>>f1;
c[0]=a[n]=0;
fxym[0]=2*f0;fxym[n]=2*f1;
break;
default:
cout<<"不可用\n";
};
for（i=1;ifxym[i]=6*f（i-1,i,i+1）;
for（i=1;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<<"Doyouhaveanothertry?
y/n:
";
cin>>ch;
}while（ch=='y'||ch=='Y'）;
return0;
}
voidprintout（intn）
{
cout<for（inti=0;i{
cout<["<cout<<"S"<floatt=fxym[i]/（6*h[i]）;
if（t>0）cout<<-t<<"*（x-"<elsecout<<-t<<"*（x-"<t=fxym[i+1]/（6*h[i]）;
if（t>0）cout<<"+"<elsecout<<"-"<cout<<"\n\t";
t=（y[i]-fxym[i]*h[i]*h[i]/6）/h[i];
if（t>0）cout<<"-"<elsecout<<"-"<<-t<<"*（x-"<t=（y[i+1]-fxym[i+1]*h[i]*h[i]/6）/h[i];
if（t>0）cout<<"+"<elsecout<<"-"<<-t<<"*（x-"<cout<}
cout<}
6.M次多项式曲线拟合
6.1算法说明
设
有N个点，横坐标是确定的。
最小二乘抛物线的系数表示为
求解A，B和C的线性方程组为
6.2M
升级会员 