数值分析作业.docx
《数值分析作业.docx》由会员分享，可在线阅读，更多相关《数值分析作业.docx（20页珍藏版）》请在冰豆网上搜索。
数值分析作业
数值分析作业
（鲁鹏：
2009020979岩土工程）
1．已知A和b
1.1用超松弛法求解Ax=b（取松弛因子
=1.4，
=0，迭代9次）
解：
1.1.1原理：
其基本思想是在高斯方法已求出x（m）,x（m-1）的基础上，经过重新组合的新序列，而此新序列收敛速度加快。
其算式是：
xi（m）=（1-ω）xi（m-1）+ω（
bijxi（m）+
xj（m-1）+gi）
其中ω是超松弛因子，当ω>1时，可以加快收敛速度。
1.1.2.程序：
#include
#include
#include
#definesorc1.4
#defineA9
voidmain（）
{inti,j,k;
doubleaa[A][A]={{12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743},
{2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124},
{-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103},
{1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585},
{-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317},
{0.71819,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417},
1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2.213474},
{3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782},
{-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001}};
doublebb[A]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};
doublexl[A],x[9],tem=0.0;
clrscr（）;
for（i=0;i{xl[i]=0;
x[i]=0;
}
for（k=0;k<9;k++）
{for（i=0;i{tem=0;
for（j=0;j{if（j!
=i）
tem=tem-aa[i][j]*xl[j];
}
x[i]=（tem+bb[i]）/aa[i][i];
x[i]=xl[i]+sorc*（x[i]-xl[i]）;
xl[i]=x[i];
}
}
printf（"\n\n\t\ttheresultsareshownbelow!
\n\n\n"）;
for（i=0;iprintf（"\t\tx[%d]=%lf\n",i+1,x[i]）;
printf（"\n\n\t\tpressanykeytoexit!
\n"）;
getch（）;
}
1.1.3输出结果：
1.1.4讨论：
超松弛法，对G-S方法每一步得到的结果与前一步所得到的结果进行了处理，迭代收敛更快，实践证明，选择恰当的松弛系数不仅可以加快收敛速度，而且可以得到满足精度的解答！
1.2请用列主元素消去法求解Ax=b：
1.2.1原理：
对矩阵作恰当的调整，选取绝对值尽量大的元素作为主元素。
然后进行行变换，把矩阵化为上三角阵，再进行回代，求出方程的解。
1.2.2程序：
#include
#include
#include
voidmain（）
{inti,j,k,l,m,n,r,maxi;
doubletem,li,x[9];
doubleaa[9][9]={{12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743},
{2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124},
{-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103},
{1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585},
{-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317},
{0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417},
{1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2.213474},
{3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782},
{-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001}};
doublebb[9]
={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};
clrscr（）;
for（i=0;i<8;i++）
{maxi=i;
for（j=i;j<9;j++）
if（abs（aa[j][i]）>abs（aa[maxi][i]））
for（k=i;k<9;k++）
{tem=aa[i][k];
aa[i][k]=aa[maxi][k];
aa[maxi][k]=tem;
}
for（l=i+1;l<9;l++）
{li=aa[l][i]/aa[i][i]*（-1）;
for（m=i;m<9;m++）
aa[l][m]=aa[l][m]+aa[i][m]*li;
bb[l]=bb[l]+bb[i]*li;
}
}
x[8]=bb[8]/aa[8][8];
for（i=7;i>=0;i--）
{tem=0;
for（j=8;j>i;j--）
tem=tem+aa[i][j]*x[j];
x[i]=（bb[i]-tem）/aa[i][i];
}
for（i=0;i<9;i++）
printf（"\t\tx[%d]=%lf\n",i,x[i]）;
getch（）;
}
1.2.3结果输出：
1.2.4讨论：
列主元素消去法减小了误差传递，结果更为精确！
3．已知函数如下表：
x
1
2
3
4
5
f（x）
0
0.69314718
1.0986123
1.3682944
1.6094378
x
6
7
8
9
10
f（x）
1.7917595
1.9459101
2.079445
2.1972246
2.3025851
f'（x）
f'
（1）=1
f'（10）=0.1
试用三次样条差值求f（4.563）及f’（4.563）的近似值
3.1原理：
根据对任意分划的三弯矩插值法：
把区间分为n个点，
取区间[
],在这个区间上,s（x）是三次多项式，s”（x）是一次多项式，设
（待定常数）,则s”（x）是过（
），（
）两点的直线。
则
（记
）
将s”（x）积分两次，利用
确定任意常数，得：
记
边界条件：
则
由于x的划分为均匀的，
则
（其中对于
来说j=1~8）则矩阵可化为：
3.2程序：
#include
#include
#include
#defineN10
voidmain（）
{inti;
doublex1,y1,mj1;
doublex[N]={1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0},h[N-1],mj[N],yd0=1,ydn=0.1,d[N],a[N-1],b[N],c[N-1],ld[N-1],bd[N],r[N-1],dt[N-1];
doubley[N]={0.0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1972246,2.3025851};
clrscr（）;
for（i=0;ib[i]=2;
for（i=0;ih[i]=x[i+1]-x[i];
d[0]=6/h[0]*（（y[1]-y[0]）/h[0]-yd0）;
for（i=1;id[i]=6/（h[i-1]+h[i]）*（（y[i+1]-y[i]）/h[i]-（y[i]-y[i-1]）/h[i-1]）;
d[N-1]=6/h[N-2]*（ydn-（y[N-1]-y[N-2]）/h[N-2]）;
for（i=0;i{a[i]=h[i]/（h[i]+h[i+1]）;
c[i]=1-a[i];
}
bd[0]=b[0];
for（i=0;i{r[i]=c[i];
ld[i]=a[i]/b[i];
bd[i+1]=b[i+1]-r[i]*ld[i];
}
dt[0]=d[0];
for（i=1;idt[i]=d[i]-ld[i-1]*dt[i-1];
mj[N-1]=dt[N-1];
for（i=N-2;i>=0;i--）
mj[i]=（dt[i]-r[i]*mj[i+1]）/bd[i];
x1=4.563;
printf（"\n\n\ttheresultsareshownbelow!
\n\n"）;
for（i=0;iif（x1>=x[i]&&x1<=x[i+1]）
{y1=mj[i]*pow（（x[i+1]-x1）,3）/（6*h[i]）+mj[i+1]*pow（（x1-x[i]）,3）/（6*h[i]）+（y[i]-mj[i]*pow（h[i],2）/6）*（x[i+1]-x1）/h[i]+（y[i+1]-mj[i+1]*pow（h[i],2）/6）*（x1-x[i]）/h[i];
mj1=pow（（x1-x[i]）,2）/（2*h[i]）*mj[i+1]-mj[i]*pow（（x[i+1]-x1）,2）/（2*h[i]）+（mj[i]-mj[i+1]）*h[i]/6+（y[i+1]-y[i]）/h[i];
printf（"\tf（%lf）=%lf\n\tdf（%lf）=%lf\n",x1,y1,x1,mj1）;
break;
}
printf（"\n\tpleasepressanykeytoexittheprogramme!
\n"）;
getch（）;
}
3.3结果输出：
3.4讨论：
其基本思想是对均匀分划的插值函数的构造,三次样条函数空间中不取1,x,,x2,x3,（x-xj）+3为基函数,而取B样条函数Ω3（x-xj/h）为基函数.由于三次样条函数空间是N+3维的,故我们把分点扩大到X-1,XN+1,则任意三次样条函数可用Ω3（x-xj/h）线性组合来表示S（x）=
cjΩ3（x-xj/h）这样对不同插值问题,若能确定cj由解的唯一性就能求得解S（x）.
样条插值因为没有Runge现象，计算比高次插值函数简单，效果比高次插值函数好。
且光滑性比低次分段插值函数好。
它集中了高次插值函数，低次插值函数，低次分段插值函数的优点。
4.用牛顿法求方程：
在（0.1，1.9）中的近似根（初始近似值取为区间断点，迭代6次货误差小于0.00001）
4.1原理：
设函数在有限区间[a,b]上二阶导数存在,且满足条件
ⅰ）f（a）f（b）<0
ⅱ）f”（x）在区间[a,b]上不变号.
ⅲ）f’（x）≠0
ⅳ）|f（c）|/b-a≤|f’（c）|其中c是a,b中使min[|f’（a）,f’（b）]达到的一个,则对任意时近似值x0€[a,b],Newton迭代过程为
xk+1=φ（xk）=x-f（xk）/f’（xk）,k=1,2,3…
用Newton法求7次方程的根，当N=1时，有：
满足Newton法的条件,所以设其初始值为1.9，故方程为：
x7-28*x4+14=0,
Newton法迭代公式为：
如此一次一次的迭代，逼近x的真实根。
当前后两个的差<=ε时，就认为求出了近似的根。
本程序用Newton法求代数方程（最高次数不大于10）在（a,b）区间的根
4.2程序：
#include
#include
#include
doublef1（doublex）
{doubley;
y=pow（x,7）-28*pow（x,4）+14;
returny;
}
doublef2（doublex）
{doubley;
y=7*pow（x,6）-28*4*pow（x,3）;
returny;
}
voidmain（）
{intk=0;
doublex1,x2,f1,f2,tep,i;
clrscr（）;
printf（"\n\n\tthewholeprocessisshownbelow:
\n"）;
for（i=0.1;i<2.0;i=i+0.1）
{x1=i;
do
{
f1=pow（x1,7）-28*pow（x1,4）+14;
f2=7*pow（x1,6）-28*4*pow（x1,3）;
tep=f1/f2;
x2=x1-tep;
x1=x2;
}
while（fabs（tep）>0.00001）;
printf（"\tx0=%1.1lf---->x=%lf",i,x2）;
k=k+1;
if（k%2==0）
printf（"\n"）;
}
getch（）;
}
4.3结果输出：
4.4讨论：
由以上不同初值的迭代结果可见，取不同的初值，该迭代法会收敛于不同的解.
5,用Romberg算法求
（允许误差
=0.00001）
5.1原理：
当
就停止运算。
5.2程序：
#include
#include
#include
#defineA1
#defineB3
doublef（doublex）
{doubley;
y=pow（3,x）*pow（x,1.4）*（5*x+7）*sin（pow（x,2））;
returny;
}
voidmain（）
{inti,k=1,j;
doublet[10][10],tep;
clrscr（）;
t[0][0]=（B-A）/2*（f（A）+f（B））;
t[0][1]=（t[0][0]+（B-A）*f（（B+A）/2））/2;
t[1][0]=（4*t[0][1]-t[0][0]）/3;
printf（"\n\n\ttheprocessisshownbelow!
\n\tt[0][0]=%lf\n",t[0][0]）;
printf（"\tt[1][0]=%lf\n",t[1][0]）;
while（fabs（t[k][0]-t[k-1][0]）>0.00001）
{k=k+1;
tep=0;
for（j=1;j<=pow（2,k-1）;j++）
tep=tep+f（A+（2*j-1）*（B-A）/pow（2,k））;
t[0][k]=（t[0][k-1]+（B-A）/pow（2,k-1）*tep）/2;
for（i=1;i<=k;i++）
t[i][k-i]=（pow（4,i）*t[i-1][k-i+1]-t[i-1][k-i]）/（pow（4,i）-1）;
printf（"\tt[%d][%d]=%lf\n",k,0,t[k][0]）;
}
printf（"\n\tpressanykeytoexit!
\n"）;
getch（）;
}
5.3计算过程输出：
5.4讨论：
Romberge算法的优点是:
1.把积分化为代数运算,而实际上只需求T1（i）,以后由递推法可得.
2.算法简单且收敛速度快,一般4或5次即能达到要求.
3.节省存储量.
Romberge算法的缺点是:
1.对函数的光滑性要求较高.
2.计算新分点的数值时,这些数值的个数将成倍增加.
6.用定步长四阶Runge-Kutta法求解
h=0.0005,打印
（0.025）,
（0.045），
（0.085），
（0.1），（i=1，2，3）
6.1原理：
表示
6.2程序：
#include
#include
#include
#defineH0.0005
doublef（intn,doublex,doubley1,doubley2,doubley3）
{doubletem;
if（n==0）
tem=1+0*x+0*y1+0*y2+0*y3;
if（n==1）
tem=0+0*x+0*y1+0*y2+y3;
if（n==2）
tem=1000+0*x+0*y1-1000*y2-100*y3;
returntem;
}
voidmain（）
{inti,j;
doubleyy[3]={0,0,0};
doublek[3][4],xx=0.0005,tt[3]={0,0,0},tem;
clrscr（）;
printf（"\n\n\n\t\t\ttheresultsareshownbelow!
\n\n\n"）;
while
（1）
{for（i=0;i<3;i++）
{k[i][0]=f（i,xx,tt[0],tt[1],tt[2]）;
tt[i]=yy[i]+H*k[i][0];
k[i][1]=f（i,xx+H/2,tt[0],tt[1],tt[2]）;
tt[i]=yy[i]+H*k[i][1]/2;
k[i][2]=f（i,xx+H/2,tt[0],tt[1],tt[2]）;
tt[i]=yy[i]+H*k[i][2];
k[i][3]=f（i,xx+H/2,tt[0],tt[1],tt[2]）;
for（j=0;j<3;j++）
tt[j]=yy[j];
}
xx=xx+H;
for（i=0;i<3;i++）
{yy[i]=yy[i]+H/6*（k[i][0]+2*k[i][1]+2*k[i][2]+k[0][3]）;
tt[i]=yy[i];
if（fabs（xx-0.025）<=0.00005||fabs（xx-0.045）<=0.000005||fabs（xx-0.085）<=0.000005||fabs（xx-0.1）<=0.000001）
{printf（"\tyy%d（%lf）=%f",i,xx,yy[i]）;
if（i==2）
printf（"\n"）;
}
}
if（xx>0.1）
break;
}
getch（）;
}
6.3输出结果：
6.4讨论：
Runge_Kutta方的优点:
1.精度高,不必用别的方法求开始几点的函数值.
2.可根据f’（t,y）变化的情况与需要的精度自动修改步长.
3.方法稳定.
4.程序简单,存储量少.
Runge_Kutta方的缺点:
每步要计算函数值f（t,y）四次,在f（t,y），较复杂时,工作量大,且每一步缺乏可靠的检查.
另外，仅就本程序而言，在写函数的时候另外引入了一个整形变量，方便了循环体中套用，是本人的创新之处！
学生：
鲁鹏
二零零九年十二月廿七日
升级会员 