范文一:数值分析上机报告
第一题:
1、已知A 与b
?12.38412? 2.115237?
? -1.061074?
? 1.112336A =? -0.113584
?
? 0.718719? 1.742382?
? 3.067813? -2.031743?
1.742382 3.067813 -2.031743?-0.784165 1.112348 3.123124??
-3.12543215.567914 3.123848 2.0314541.836742-1.056781 0.336993 -1.010103?
?
-1.012345 3.12384827.108437 4.101011-3.741856 2.101023 -0.71828 -0.037585? 2.189736 2.031454 4.10101119.8979180.431637-3.111223 2.121314 1.784317?
?
1.563849 1.836742 -3.741856 0.4316379.789365-0.103458 -1.103456 0.238417? -0.784165 -1.056781 2.101023-3.111223-0.10345814.7138465 3.123789 -2.213474?
?
1.112348 0.336993-0.71828 2.121314-1.1034563.12378930.719334 4.446782? 3.123124 -1.010103-0.037585 1.7843170.238417-2.213474 4.44678240.00001?? 2.11523719.141823
-1.061074 1.112336 -0.113584 -3.125432 -1.012345 2.189736
0.7187191.563849
B =[ 2.1874369 33.992318 -25.173417 0.84671695 1.784317 -86.612343 1.1101230 4.719345 -5.6784392]T
(1)用Househloser 变换,把A 化为三对角阵(并打印B )。
(2)用超松弛法求解Bx=b(取松弛因子ω=1.4,x (0)=0,迭代9次)。 (3)用列主元素消去法求解Bx=b。
一、分析如下:
(3)用列主元素消去法求解Bx=b。
将方阵A 和向量b 写成C=(A b)。将C 的第1列中第1行的元素与其下面的此列的元素逐一进行比较,找到最大的元素c j 1,将第j 行的元素与第1行的元素进行交换,然后通过行变换,将第1列中第2到第n 个元素都消成0。将变换后的矩阵C (1)的第二列中第二行的元素与其下面的此列的元素逐一进行比较,找
(1)
到最大的元素c k 将第k 行的元素与第2行的元素进行交换,然后通过行变换,2,
将第2列中第3到第n 个元素都消成0。以此方法将矩阵的左下部分全都消成0。最终形式如下:
(n ) ?a 11 0 (A b)~ 0?
n )
g 1? a 1(n
?
* ?
?
?(n )
g n ? 0a nn
二、程序:
#include "math.h"
#include "stdio.h" /* 标准的基本库函数 头文件 */ #define ge 8 void main() {
int sign(double x); /*该函数为sign 函数*/
double t[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.713846,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} }; double p,v,s,w; int i,j,h,m,g; double u[9],x1[9],y[9],q[9],b1[9][10],x[9]; double d[9]=
{2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};
/*--------------------开始迭代----------------------*/ for(j=0;j<7;++j) household="" 变换="" */="" {="" s="0.0;" for(i="">7;++j)><9;++i) s="s+t[i][j]*t[i][j];" s="sqrt(s);" v="(t[j+1][j]">0)?(s*s+s*t[j+1][j]):(s*s-s*t[j+1][j]); for(g=0;g<9;++g) {="" if="">9;++g)><=j) u[g]="0;" else="" if="" (g="=j+1)" u[g]="t[j+1][j]+s*sign(t[j+1][j]);" else="" u[g]="t[g][j];" }="" for(m="">=j)><>
{ y[m]=0; for(h=0;h<9;++h) y[m]="y[m]+t[m][h]*u[h];" y[m]="y[m]/v;" }="" p="0;" for(i="">9;++h)><9;++i) p="p+u[i]*y[i];" p="0.5*p/v;" for(i="">9;++i)><9;++i) q[i]="y[i]-p*u[i];" for(h="">9;++i)><9;++h) for(m="">9;++h)><9;++m) t[m][h]="t[m][h]-(q[m]*u[h]+u[m]*q[h]);">9;++m)>
printf("Household:\h"); for(i=0;i<9;++i) for(j="">9;++i)><9;++j) {="" if="" (j%9="=0)" printf("\n");="" printf("%-9.5f",t[i][j]);="" }="">9;++j)>
w=1.4; /*超松弛法*/ for(i=0;i<9;i++) 开始循环="" */="" x1[i]="0;" for(i="">9;i++)><9;i++) for(j="">9;i++)><9;j++) {="" if(i="=j)" b1[i][j]="0;" else="" b1[i][j]="-t[i][j]/t[i][i];}" for(i="">9;j++)><9;i++) b1[i][9]="d[i]/t[i][i];" for(h="">9;i++)><9;h++) for(i="">9;h++)><9;i++) {="" s="0;" for(j="">9;i++)><9;j++) s="s+b1[i][j]*x1[j];" s="s+b1[i][9];" x1[i]="">9;j++)>
} for(i=0;i<9;i++) {="" if="" (i="=5)" printf("\n");="" printf("x%d="%-10.6f",i,x1[i]);" 输出结果*/="" }="" printf("\n");="" u[0]="t[0][0];" 以下是消去法*/="" y[0]="d[0];" for(i="">9;i++)><9;i++) {="" q[i]="t[i][i-1]/u[i-1];" u[i]="t[i][i]-q[i]*t[i-1][i];" y[i]="d[i]-q[i]*y[i-1];" }="" x[ge]="y[ge]/u[ge];" for(i="ge-1;i">=0;i--) x[i]=(y[i]-t[i][i+1]*x[i+1])/u[i]; for(i=0;i<9;i++) {="" if="" (i="=5)" printf("\n");="" printf("="" x%d="%-10.6f",i,x[i]);" 输出结果*/="" }="">9;i++)>
int sign(double x) { int z; /*定义变量*/ z=(x>=(1e-6)?1:-1); return(z); /* 主函数返回z 告诉系统 程序结束 */ }
三、运行结果:
四.问题讨论:
题中矩阵B 是对称正定阵,且是三对角的,所以选择合适的松驰因子ω,收敛速度是很快的。已经求出x 的第m 解,第m-1基础上,经过重新组合得新的序列,而在此新序列收敛速度加快。
第三题:
3、已知函数值如下表:
一、分析如下:
S (x ) =M j -1
(x j -x ) 36h j -1
+M j
(x -x j -1) 3
6h j -1h 2x j -x h 2x -x j -1j -1j -1
+(y j -1-M j -1)() +(y j -M j )()
6h j -16h j -1
这里h j -1=x j -x j -1 ,所以只要求出M j ,就能得出插值函数S (x )。
?21
?μ
?12λ1?μ22
求M j 的方法为:?
? ????
λ2
μN -1
??M 0??d 0?
??M ??d ???1??1?????????=??
?????
???2λN -1???????
12????M N ????d N ??
6y 1-y 0?
') d =?0h (h -y 0
00
??y j +1-y j y j -y j -16d =(-) (j =1,2, , N -1) ?j
h j -1+h j h j h j -1
?这里?
?d =6[y '-1(y -y )]
N -1
?N h N -1N h N -1N ?
h j
?μ=h j -1
λj =1-μj =?j h +h h j -1+h j j -1j ?
最终归结为求解一个三对角阵的解。
用追赶法解三对角阵的方法如下:
??b 1c 1??1??β1γ1
??a b ??l 1??c γβ222122??????
?=LU ?=?l 21?? A =?
??????
b c a βγn -1n -1n -1n -1n -1??????
?l n 1?a n b n ?βn ??????????δ1?
L δ=d ??, 则由L δ=d 得
LUx =d , 即?, 若记δ=? ??
?Ux =δ??δn ??
?1??δ1??d 1?
?l ?? ?? ?12????=?? , ? ?? ?? ???????
l n 1??δn ??d n ??
综上可得求解方程Ax=d的算法:
?β1
?????
γ1
βn -1
??x 1??δ1??? ?? ????=?? γn -1?? ?? ?
?????βn ??x n ??δn ?
αi +1?
β=b , δ=d , l =, βi +1=b i +1-l i +1c i , 11i +1?11
βi
??
i =1, 2,3, , n -1 ?δi +1=d i +1-l i +1δi
?δδ-c x
?x n =n , x i =i i i +1, i =n -1, , 2,1
βn βi ??
二、程序如下:
#include double th[10][10]; int i; double d[10]; double y[10]={0,0.69314718,1.0986123,1.3862944, 1.6094378,1.7917595,1.9459101,2.079445,2.1972246,2.3025851}, xm[10]={1,2,3,4,5,6,7,8,9,10}; double h=1;a; double x[10]; double b[10],l[10]; double c[10],c1[10],c2[10]; float f[10],f1[10]; for(i=0;i<10;i++) th[i][i]="2;" th[0][1]="">10;i++)> for(i=1;i<9;i++)>9;i++)> th[i][i+1]=0.5; th[i][i-1]=0.5; } th[9][8]=1; \*输入矩阵值*\ d[0]=6*((y[1]-y[0])/h-1)/h; for(i=1;i<> d[i]=6*((y[i+1]-y[i])/h-(y[i]-y[i-1])/h)/(2*h); d[9]=6*(0.1-(y[9]-y[8])/h)/h; \*计算d 值*\ b[0]=th[0][0]; c[0]=d[0]; for(i=1;i<10;i++)>10;i++)> l[i]=th[i][i-1]/b[i-1]; b[i]=th[i][i]-l[i]*th[i-1][i]; c[i]=d[i]-l[i]*c[i-1]; } x[9]=c[9]/b[9]; for(i=8;i>=0;i--) x[i]=(c[i]-th[i][i+1]*x[i+1])/b[i]; \*用追赶法求解三样条插值法中的M*\ printf("please input:"); \*输入x 值*\ scanf("%lf",&a); for(i=1;i<10;i++)>10;i++)> c1[i]=(x[i-1]*pow(h,2)/6+y[i]-x[i]*pow(h,2)/6-y[i-1])/h; \*求出c1值*\ c2[i]=y[i]-x[i]*pow(h,2)/6; \*求出c1值*\ f[i]=x[i-1]*pow((xm[i]-a),3)/6+x[i]*pow((a-xm[i-1]),3)/6+c1[i]*(a-xm[i])+c2[i]; \*f(x)表达式*\ f1[i]=x[i]*pow((a-xm[i-1]),2)/2+c1[i]-x[i-1]*pow((xm[i]-a),2)/2; \*f’(x) 表达式*\ } printf("\n"); \*以下为判断求输入的自变量x 的表达式*\ if(0<><=1) printf("f="%8.7f,f1=%8.7f",f[1],f1[1]);" else="">=1)><><=2) printf("f="%8.7f,f1=%8.7f",f[1],f1[1]);" else="">=2)><><=3) printf("f="%8.7f,f1=%8.7f",f[2],f1[2]);" else="">=3)><><=4) printf("f="%8.7f,f1=%8.7f",f[3],f1[3]);" else="">=4)><><=5) printf("f="%8.7f,f1=%8.7f",f[4],f1[4]);" else="">=5)><><=6) printf("f="%8.7f,f1=%8.7f",f[5],f1[5]);" else="">=6)><><=7) printf("f="%8.7f,f1=%8.7f",f[6],f1[6]);" else="">=7)><><=8) printf("f="%8.7f,f1=%8.7f",f[7],f1[7]);" else="">=8)><><=9) printf("f="%8.7f,f1=%8.7f",f[8],f1[8]);" else="" printf("f="%8.7f,f1=%8.7f",f[9],f1[9]);">=9)> 三、计算结果如下: f(4.563)=1.517932 f'(4.563)=0.249350 四、问题讨论: 此题最终能实现对任意一组数据的三次样条插值计算得出区间内任意一个点的函数值及其导数值。 第四题: 4、用Newton 法求方程 x 7-28x 4+14=0 在(0.1 , 1.9)中的近似根(初始近似 值取为区间端点,迭伐6次或误差小于0.000 01)。 一、 分析如下: 对满足下列条件的函数: i. f(a)f(b)<> ii. f”(x)在区间[a,b]上不变号; iii. f’(x)≠0; iv. |f(c)|/(b-a)≤|f’(c)|,其中c 是a,b 中使min(|f’(a)|,|f’(b)|)达到的一个。 则对任意初始近似值x 0∈[a , b ],迭代方程为: x k +1=?(x k ) =x k - f (x k ) , k =0,1,2 f '(x k ) 二、程序如下: #include double x1=1.9,x2, ε,f, f’; do{ f=pow(x1,7)-28*pow(x1,4)+14; f ’=7*pow(x1,6)-28*4*pow(x1,3); x 2=x1-f/f’; ε=x1-x 2; x 1=x2; } while (fabs(ε)>=pow(10,-6)); printf(“%f”,x 2); } 三、计算结果如下: equation is : 1.00 x^7 - 28.00 x^4 + 14.00 = 0 x = 0.8454966 四、问题讨论: 当x=1.9时运行结果在[0.1,1.9]之间,是要求的解;当x=0.1时运行结果不在[0.1,1.9] 之间,不符合要求。这并不是说当x=0.1时所得解不是原方程的解,而是所得解不在所要求的范围内。程序对不同的方程求解,或有不同的误差要求,只需对do while 循环中量做适当的修改即可,程序没有考虑f’(x)=0的情况,有可能产生死循环。 第五题: 5、用Romberg 算法求?13x x 1.4(5x +7)sin x 2dx (允许误差ε=0.000 01)。 一、本题分析如下: Romberg 算法的计算步骤如下: (1)先求出按梯形公式所得的积分值 3 T 1(0)= b -a [f (a ) +f (b )] 2 (2)把区间2等分,求出两个小梯形面积之和,记为T 1(1),即 b -a b +a [f (a ) +f (b ) +2f ()] 42 1 T 1(1)-() 2T 1(0) 4T 1(1)-T 1(0)(0)(0)这样由外推法可得T 2, T 2=。 = 124-11-() 2T 1(1)= (3)把区间再等分(即22等分),得复化梯形公式T 1(2),由T 1(1)与T 1(2)外推可得 T (1)2 4T 1(2) -T 1(1)42T 2(1)-T 2(0)(0)=,T 3=,如此,若已算出2k 等分的复化梯形公式T 1(k ) ,2 4-14-1 则由Richardson 外推法,构造新序列 T (k -1) m +1 (k ) (k -1) 4m T m -T m , m=1,2,…,l, k=1,2,…,l-m+1, =m 4-1 最后求得T l (0)+1。 (0)(4)T l (0)≈T l (0),计算T 1(l +1) ,一般可用如-T l (0)+1或|T l +1|<> 下算法: ?(0)b -a [f (a ) +f (b )]?T 1=2? 2l -1?b -a ?(l ) 1(l -1) b -a T ={T +f [a +(2i -1) ]}?11l -1∑l 222i =1? (k ) (k -1) ?4m T m -T m (k -1) ?T m +1=, m =1, 2, , l , k =1, 2, , l -m +1m 4-1?? 其具体流程如下,并全部存入第一列 (1)T 1(0) (3T ) (6)T 3(0) (2)T 1(1)(5T ) (4)T 1Romberg 算法的一些说明: 通常计算时,用固定L=N来计算,一般L=4或5即能达到要求。 应用Romberg 算法时,可以先给出所要求的误差ε,当|T l (0)-T l (0)+1|<ε就停机。 romberg="" 算法的优点是:(1)把积分化为代数运算,而实际上只需求t="" 1(l="" )="">ε就停机。> 推可得。(2)算法简单且收敛速度快。 Romberg 算法的缺点是:(1)对函数的光滑性要求较高。(2)计算新分点的值时,这些数值的个数成倍增加。 二 、程序如下: #include "stdio.h" #include"math.h" main() {int i,k,e; double ya,yb,a,b,tt,h,m,eps; double t[10][10]; a=1; b=3; h=b-a; m=0; tt=0; eps=0.001; ya=rem(a); yb=rem(b); t[0][1]=(ya+yb)*h/2; for(k=1;k<> {for(i=1;i<> {m=a+(i-0.5)*h/pow(2,k-1); tt=tt+rem(m) ; m=0; } t[k][1]=t[k-1][1]+h*tt/pow(2,k-1); t[k][1]=t[k][1]/2; tt=0; } for(k=1;k<> { for(i=0;i<> { t[i][k+1]=(pow(4,k)*t[i+1][k]-t[i][k])/(pow(4,k)-1); } if(abs(t[0][k+1]-t[0][k])<> break; } printf("Romberg算法求得的积分结果为: %e \n",t[0][k+1]); } rem(x) double x; {double y; y=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x); return(y);} 三、计算结果如下: The result: y=440.535928 四、问题讨论: (l ) (l ) 本程序可通过一个一维数组将T m 进行存储,但考虑到对程序的可读性,故将T m 存 储在一个二维数组里。 第六题: 6、用定步长四阶Runge-Kutta 法求解 ?dy 1/dt =1?dy /dt =y 3?2 ??dy 3/dt =1000-1000y 2-100y 3 ?y (0)=0?1 ?y 2(0)=0???y 3(0)=0 h=0.000 5,打印y i (0.025),y i (0.045),y i (0.085),y i (0.1),(i =1,2,3) 。 一、分析如下: 1?Y =Y +?n +1n 6(K 1+2K 2+2K 3+K 4) ??K 1=hF (x n , Y n ) ?11? K =hF (x +h , Y +K 1) ?2n n 22?11?K =hF (x +h , Y +K 2) n n ?322???K 4=hF (x n +h , Y n +K 3) 二、 程序如下: #include #include #define N 3 float y[N]={0,0,0},k[N][N+1]; float yd[N+1][N],a[N+1]={0.025,0.045,0.085,0.10}; /*以上是预处理和定义全局变量*/ main() {int i,l,j=0,n=1; float h=0.0005; do{for(l=0;l<> for(i=0;i<=n;i++) 计算每一个k="">=n;i++)> k[l][i]=0; k[0][0]=h; k[1][0]=h*y[2]; k[2][0]=h*(1000-1000*y[1]-100*y[2]); k[0][1]=h; k[1][1]=h*(y[2]+(1/3.0)*k[0][0]); k[2][1]=h*(1000-1000*(y[1]+(1/3.0)*k[1][0])-100*(y[2]+(1/3.0)*k[1][0])); k[0][2]=h; k[1][2]=h*(y[2]+(1/3.0)*k[1][0]+k[1][1]); k[2][2]=h*(1000-1000*(y[1]+(1/3.0)*k[2][0]+k[2][1])-100*(y[2]+(1/3.0)*k[2][0]+ k[2][1])); k[0][3]=h; k[1][3]=h*(y[2]+k[1][0]-k[1][1]+k[1][2]); k[2][3]=h*(1000-1000*(y[1]+k[2][0]-k[2][1]+k[2][2])-100*(y[2]+k[2][0]-k[2][1]+ k[2][2])); for(i=0;i<> y[i]=y[i]+(k[i][0]+3*k[i][1]+3*k[i][2]+k[i][3])/8; if(n==50||n==90||n==170||n==200)/*得到所需的数据*/ {for(i=0;i j++;} n++; }while(n<> printf("The result of sixth question is:\n"); for(j=0;j<> {for(i=0;i<> printf("y%d(%5.4f)=%.7f,",i+1,a[j],yd[j][i]); printf("\b\n");} } 三、 计算结果如下: y1(0.000)=0.000000 y2(0.000)=0.000000 y3(0.000)=0.000000 y1(0.025)=0.025000 y2(0.025)=0.151599 y3(0.025)=8.335426 y1(0.045)=0.045000 y2(0.045)=0.312860 y3(0.045)=7.536280 y1(0.085)=0.085000 y2(0.085)=0.560581 y3(0.085)=4.946353 y1(0.100)=0.100000 y2(0.100)=0.628881 y3(0.100)=4.180993 四、问题讨论: 在设计本程序初期,曾将再次计算时的初始值赋值给忽略了,造成只有第一个计算结果是正确的,而不能继续进行其它数据的计算,后经修改,可实现对其它数据的连续计算。在设计过程中,由于将浮点数赋值给了整数,造成了计算误差,后经修改,使计算正确进行。 本程序可以通过修改主程序所调用的函数中的表达式来实现对其它函数的任意初值条件求微分计算。 数值分析上机报告 班级:20级学隧2班 姓名:000000000 学号:00000000000 目录 1 序言 ................................................................................................................... 6 2 题目 ................................................................................................................... 7 2.1 题2 ......................................................................................................... 7 2.1.1 题目内容 ..................................................................................... 7 2.1.2 MATLAB程序 ................................................................................ 8 2.1.3 计算结果 ..................................................................................... 8 2.1.4 图形 ............................................................................................. 9 2.1.5 分析 ................................................................................................... 14 2.2 题3 ....................................................................................................... 14 2.2.1 题目内容 ................................................................................... 14 2.2.2 程序 ........................................................................................... 14 2.2.3 计算结果 ................................................................................... 14 2.2.4 图形 ........................................................................................... 15 2.2.5 分析 ........................................................................................... 16 2.3 选做题5 ............................................................................................... 16 2.3.1方法介绍 .................................................................................... 17 2.3.2计算结果及分析 ........................................................................ 17 3总结 .................................................................................................................. 18 4.附录 ............................................................................................................... 19 4.1 题1程序代码 ...................................................................................... 19 4.2 题2程序代码 ...................................................................................... 22 4.3 题3程序代码 ...................................................................................... 26 数值分析2015上机实习报告要求 1.应提交一份完整的实习报告。具体要求如下: (1)报告要排版,美观漂亮(若是纸质要有封面,封面上)要标明姓名、 学号、专业和联系电话; (2)要有序言,说明所用语言及简要优、特点,说明选用的考量; (3)要有目录,指明题目、程序、计算结果,图标和分析等内容所在位置, 作到信息简明而完全; (4)要有总结,全方位总结机编程计算的心得体会; (5)尽量使报告清晰明了,一般可将计算结果、图表及对比分析放在前面, 程序清单作为附录放在后面,程序中关键部分要有中文说明或标注,指明该部分的功能和作用。 2.程序需完好保存到期末考试后的一个星期,以便老师索取用于验证、询问或质疑部分内容。 3.认真完成实验内容,可以达到既学习计算方法又提高计算能力的目的,还可以切身体会书本内容之精妙所在,期间可以得到很多乐趣。 4.拷贝或抄袭他人结果是不良行为,将视为不合格。 5.请按任课老师要求的时间和载体(电子或纸质)提交给任课老师。 数值分析上机试题 (请在1-4题中选择两个题目,5-6中选择一个题目) 1. 分别用牛顿法,及基于牛顿算法下的Stef fensen 加速法 (1) 求ln(x +sinx )=0的根。初值x0分别取0.1, 1,1.5, 2, 4进行计算。 (2) 求sin x =0的根。初值x0分别取1,1.4,1.6, 1.8,3进行计算。 分析其中遇到的现象与问题。 2. 某过程测涉及两变量x 和y, 拟分别用插值多项式和多项式拟合给出其对应规律的近似多项式,已知xi 与yi 之间的对应数据如下,xi=1,2,…,10 yi = 34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795 103.5743 97.4847 78.2392 (1)请用次数分别为3,4,5,6的多项式拟合并给出最好近似结果f(x)。 (2)请用插值多项式给出最好近似结果 下列数据为另外的对照记录,它们可以作为近似函数的评价参考数据。 xi = Columns 1 through 7 1.5000 1.9000 2.3000 2.7000 3.1000 3.5000 3.9000 Columns 8 through 14 4.3000 4.7000 5.1000 5.5000 5.9000 6.3000 6.7000 Columns 15 through 17 7.1000 7.5000 7.9000 yi = Columns 1 through 7 42.1498 41.4620 35.1182 24.3852 11.2732 -1.7813 -12.3006 Columns 8 through 14 -18.1566 -17.9069 -11.0226 2.0284 19.8549 40.3626 61.0840 Columns 15 through 17 79.5688 93.7700 102.3677 3.用雅格比法与高斯-赛德尔迭代法解下列方程组Ax =b ,研究其收敛性,上机验证理论分 析是否正确,比较它们的收敛速度,观察右端项对迭代收敛有无影响。 (1)A 行分别为A 1=[6,2,-1],A 2=[1,4,-2],A 3=[-3,1,4];b 1=[-3,2,4]T , b 2=[100,-200,345]T , (2) A 行分别为A 1=[1,0,8,0.8],A 2=[0.8,1,0.8],A 3=[0.8,0.8,1];b 1=[3,2,1]T , b 2=[5,0,-10]T , (3)A 行分别为A 1=[1,3],A 2=[-7,1];b =[4,6]T , 4. 松弛因子对SOR 法收敛速度的影响。 用SOR 法求解方程组Ax =b , 其中 ?-41??-3? ? ? 1-41? -2? ? -2?. . . ?, B = ? ? . ?. . . ? ?1-41 ? -2? -3?1-4????? 要求程序中不存系数矩阵A , 分别对不同的阶数取w=1.1, 1.2, ...,1.9进行迭代,记录近似解x (k)达到||x (k)-x (k-1)||<10-6时所用的迭代次数k ,观察松弛因子对收敛速度的影响,并观察当w="" ≤0或w="">10-6时所用的迭代次数k> 5. 用Ru n ge-Kutt a 4阶算法对初值问题y /=-20*y ,y (0)=1按不同步长求解,用于观察稳定区 间的作用,推荐两种步长h=0.1,0.2。 注:此方程的精确解为:y =e-20x 6. 实验内容 (1) 实际验证梯形求积公式、Simpson 求积公式、Newton-Cotes 求积公式的代数精度。 (2) 针对上述三个函数和积分区间[a,b],实验观察梯形求积公式、Simpson 求积公式和 Newton-Cotes 求积公式的复化求积公式的实际计算效果。 y=exp(-x.^2).*sin(10*x)+4; a=1; b=3; y=sin(5*x)./x.^3;a=2*pi;b=4*pi; y=sin(5*x)./x.^3;a=2*pi;b=9.4248; 1 序言 MATLAB 由一系列工具组成。这些工具方便用户使用MATLAB 的函数和文件,其中许多工具采用的是图形用户界面。包括MATLAB 桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。随着MATLAB 的商业化以及软件本身的不断升级,MATLAB 的用户界面也越来越精致,更加接近Windows 的标准界面,人机交互性更强,操作更简单。而且新版本的MATLAB 提供了完整的联机查询、帮助系统,极大的方便了用户的使用。简单的编程环境提供了比较完备的调试系统,程序不必经过编译就可以直接运行,而且能够及时地报告出现的错误及进行出错原因分析。 总的来说,该软件有三大特点。一是功能强大。具有数值计算和符号计算、计算结果和编程可视化、数学和文字统一处理、离线和在线计算等功能;二是界面友善、语言自然。MATLAB 以复数处理作为计算单元,指令表达与标准教科书的数学表达式相近;三是开放性强。主要有以下优点: 1) 高效的数值计算及符号计算功能,能使用户从繁杂的数学运算分析中解脱出来; 2) 具有完备的图形处理功能,实现计算结果和编程的可视化; 3) 友好的用户界面及接近数学表达式的自然化语言,使学者易于学习和掌握; 4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ,为用户提供了大量方便实用的处理工具。 MATLAB 擅长于矩阵运算,在各种系统仿真方面应用广泛。他不同于普通的编程语言,集成有许多领域专家为各自领域开发的工具箱,直接调用即可。面向具体应用,使用MATLAB 更有针对性。这种程序在数值分析方面应用广泛,基于它的众多优点,这次我应用MATLAB 编程,把学习这个软件作为我数值分析课程的实践环节。 2 题目 2.1 题2 2.1.1 题目内容 2. 某过程测涉及两变量x 和y, 拟分别用插值多项式和多项式拟合给出其对应规律的近似多项式,已知xi 与yi 之间的对应数据如下,xi=1,2,…,10 yi = 34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795 103.5743 97.4847 78.2392 (1)请用次数分别为3,4,5,6的多项式拟合并给出最好近似结果f(x)。 (2)请用插值多项式给出最好近似结果 下列数据为另外的对照记录,它们可以作为近似函数的评价参考数据。 xi = Columns 1 through 7 1.5000 1.9000 2.3000 2.7000 3.1000 3.5000 3.9000 Columns 8 through 14 4.3000 4.7000 5.1000 5.5000 5.9000 6.3000 6.7000 Columns 15 through 17 7.1000 7.5000 7.9000 yi = Columns 1 through 7 42.1498 41.4620 35.1182 24.3852 11.2732 -1.7813 -12.3006 Columns 8 through 14 -18.1566 -17.9069 -11.0226 2.0284 19.8549 40.3626 61.0840 Columns 15 through 17 79.5688 93.7700 102.3677 2.1.2 MATLAB程序 见附录4.1 2.1.3 计算结果 2.1.3.1 插值多项式 (1)线性插值误差 Error = 10.8098 (2)三次函数插值误差 Error = 7.2441 (3)分段三次Hermit 插值误差 Error = 7.2441 (4)三次样条函数插值误差 Error = 1.0338 (5)最终结果 三次样条函数插值效果最好,误差最小。 2.1.3.2 多项式拟合 (1) 拟合次数 n = 3时 拟合误差 Error = 57.2759 3次拟合多项式 F(x) = -1.0326 x^3 + 19.3339 x^2 - 94.4787 x + 131.7944 (2) 拟合次数n = 4时 拟合误差 Error = 29.6408 4次拟合多项式 F(x) = -0.38185 x^4 + 7.368 x^3 - 42.1433 x^2 + 73.5334 x + 0.74498 (3) 拟合次数n = 5时 拟合误差 Error = 11.3278 5次拟合多项式 F(x) =0.098075 x^5 - 3.0789 x^4 + 34.502 x^3 - 163.5107 x^2 + 304.7282 x - 139.5019 (4)拟合次数 n = 6时 拟合误差 Error = 3.2509 6次拟合多项式 F(x) = 0.019359 x^6 - 0.54079 x^5 + 5.1137 x^4 - 16.8973 x^3 0.86696 x^2 + 66.375 x - 18.6991 (5)最终结果 n = 6时拟合效果最佳,误差error = 3.2509。 最佳拟合多项式为: F(x) = 0.019359 x^6 - 0.54079 x^5 + 5.1137 x^4 - 16.8973 x^3 0.86696 x^2 + 66.375 x - 18.6991 2.1.4图形 2.1.4.1 插值多项式 –– 图2-1 线性插值图像 图2-2 三次函数插值图像 图2-3 分段三次Hermit 插值图像 图2-4 三次样条函数插值图像 2.1.4.2 多项式拟合 图2-5次拟合图像 图2-6次拟合图像 图2-7次拟合图像 图2-8次拟合图像 2.1.5分析 由图2-1~2-4可以发现,插值曲线逐步变得光滑,插值误差越来越小。 由图2-5~2-8可以发现,随着拟合次数(2<><7)的不断增大,拟合曲线与检查节点间的差距越来越小,最后几乎穿过节点形成一条光滑的曲线;采用了2范数作为拟合误差控制指标,>7)的不断增大,拟合曲线与检查节点间的差距越来越小,最后几乎穿过节点形成一条光滑的曲线;采用了2范数作为拟合误差控制指标,><><7)的增大逐渐减小。说明在一定范围内,拟合曲线可以随拟合次数的增大而逼近真实曲线。 2.2="">7)的增大逐渐减小。说明在一定范围内,拟合曲线可以随拟合次数的增大而逼近真实曲线。> 2.2.1 题目内容 用雅格比法与高斯-赛德尔迭代法解下列方程组Ax =b ,研究其收敛性,上机验证理论分析是否正确,比较它们的收敛速度,观察右端项对迭代收敛有无影响。 (1)A 行分别为A 1=[6,2,-1],A 2=[1,4,-2],A 3=[-3,1,4];b 1=[-3,2,4]T , b 2=[100,-200,345]T , (2) A 行分别为A 1=[1,0,8,0.8],A 2=[0.8,1,0.8],A 3=[0.8,0.8,1];b 1=[3,2,1]T , b 2=[5,0,-10]T , (3)A 行分别为A 1=[1,3],A 2=[-7,1];b =[4,6]T , 2.2.2程序 见附录2 2.2.3计算结果 (1)雅格比法 取迭代初值矩阵=[0 0 0],经过61步迭代得出稳定结果x1 = -0.7273 x2 = 0.8081 x3 = 0.2525 (2) 高斯-赛德尔迭代法 取迭代初值矩阵=[0 0 0],经过36步迭代得出稳定结果x1 = -0.7273 x2 = 0.8081 x3 = 0.2525 2.2.4图形 图2-9 雅格比迭代法收敛图像 图2-10 高斯-赛德尔迭代法 2.2.5分析 运行程序后发现,高斯-赛德尔迭代法比起雅格比法要快得多,在取同样的初值条件下,高斯-赛德尔迭代法经过36步则迭代出稳定的结果;而雅格比法则需要61步。 另外,在右端项取b2的情况下,经过程序运行,发现经过1000次迭代仍然无法收敛,可见右端项对迭代收敛有很大影响。 2.3 选做题5 用Runge-Kutta 4阶算法对初值问题y ’=-20*y,y(0)=1按不同步长求解,用于观察稳定区间的作用,推荐两种步长h=0.1、0.2。 注:此方程的精确解为:y=e-20x 2.3.1方法介绍 通过用一些点上函数值f(x ,y)的适当线性组合,来替代Euler 法中的f(xk ,y k ) ,从而使方法的阶数更高,这就是Runge-Kutta 法的基本思想。在本次上机实习中采用的是标准4阶Runge-Kutta 法,其计算形式如下: 2.3.2计算结果及分析 利用FORTRAN 语言编写程序program R_K(附录3),计算结果如下。 表2不同步长计算结果对比 从上面结果中可以看出取步长h=0.1时候结果可以很好的收敛,与准确解之间的误差也越来越小,但是取h=0.2时候结果就发生了发散。说明该方法对步长的选择比较敏感。 3总结 通过这次实《matlab 》实验使我学习掌握了许多知识。首先是对matlab 有了一个全新的认识,其次是对matlab 的更多操作和命令的使用有了更高的掌握,最重要的事对matlab 的处理能力有了一个更高的飞跃尤其是对相关函数的使用及相关问题的处理。 通过这次上机练习,让我对数值分析所介绍的迭代求解方法及其理论有了更深层次的理解,了解了各种方法之间的优缺点,并且认识到了自己在以前的学习中所存在的问题,及时的修补了自己知识上的漏洞。同时也提高了我在编写程序上的熟练程度,所以,我认为这次上机实习是非常有收获的,给予我学习数值分 析的帮助也是非常大的。 我也希望通过这篇总结来表达自己对知道老师的感谢之情,谢谢您的不懈努力和耐心指导,才使得我再这次的实验过程中收获的这么多,也正式您的不吝教诲才使得我们在这次实验中学习和收获了许多的有用的知识和技巧,我相信在以后的学习或者工作中一定有其用武之地。过多的感谢无以言表,万分感激,百口不胜言表,至此敬礼! 4.附录 4.1 题1程序代码 (1)插值多项式源程序 clc; clear; %给定数据对 x0 = [1:1:10]; y0 = [34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 103.5743 97.4847 78.2392]; x1 = [1.5:0.4:7.9]; ycheck = [42.1498 41.4620 35.1182 24.3852 11.2732 -12.3006 -18.1566 -17.9069 ... -11.0226 2.0284 19.8549 40.3626 61.0840 79.5688 102.3677]%检查点矩阵 y1 = interp1(x0,y0,x1,'linear'); %线性插值 y2 = interp1(x0,y0,x1,'spline'); %三次样条函数插值 y3 = interp1(x0,y0,x1,'pchip'); %分段三次Hermite 插值 y4 = interp1(x0,y0,x1,'cubic'); %三次函数插值 %插值函数计算值与检算值差距计算 75.2795 -1.7813 93.7700 data1 = y1 - ycheck; data2 = y2 - ycheck; data3 = y3 - ycheck; data4 = y4 - ycheck; %求2范数作为误差指标 error1 = norm(data1,2) error2 = norm(data2,2) error3 = norm(data3,2) error4 = norm(data4,2) % 计算最小误差 error = [error1,error2,error3,error4]; min_error = error(1); for i = 1:3 if error(i+1) < min_error;="" min_error="error(i+1);" end="" end="" min_error=""> plot(x1,y1,'-b',x0,y0,'*r','markersize',5) grid on title('线性插值') figure(2) plot(x1,y2,'-b',x0,y0,'*r','markersize',5) grid on title('三次样条函数插值') figure(3) plot(x1,y3,'-b',x0,y0,'.r','markersize',15) grid on title('分段三次Hermite 插值') figure(4) plot(x1,y4,'-b',x0,y0,'.r','markersize',15) grid on title('三次函数插值') (2)多项式拟合 clc; clear; n = input('请输入拟合次数:n=') x0 = [1:1:10]; y0 = [34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795 103.5743 97.4847 78.2392]; ycheck = [42.1498 41.4620 35.1182 24.3852 11.2732 1.7813 -12.3006 -18.1566 -17.9069 ... -11.0226 2.0284 19.8549 40.3626 61.0840 79.5688 93.7700 102.3677]; %检查点矩阵 set(0,'defaultfigurecolor','w') plot(x0,y0,'b',x0,y0,'r*','markersize',10),xlabel('x'),ylabel('y') %原始点图像 p = polyfit(x0,y0,n); %求拟合多项式的系数 fx1 = poly2str(p,'x'); %以显式方式写出拟合得到的f(x) x1 = [1.5:0.4:7.9] y1 = polyval(p,x1) %由拟合多项式拟合得到的值 y2 = polyval(p,x0) %用MATLAB 库函数求2范数 data = (y2 - y0).*(y2 - y0); error = norm((y2-y0),2); %编程求2范数 mistake1 = 0; for i = 1:10 mistake1 = mistake1 + abs(data(i)); end mistake = sqrt(mistake1) ; figure(2) set(0,'defaultfigurecolor','w') plot(x1,y1,'-b',x0,y0,'.r','markersize',20) %绘出拟合图像 figure(3) set(0,'defaultfigurecolor','w') plot(x0,y2,'-b',x0,y0,'r--')%拟合图像与原图像的比较 error %拟合误差 mistake %拟合误差 fx = sprintf(fx1) %输出拟合得到的f(x)到屏幕上 4.2 题2程序代码 (1)雅可比法 clc; clear; orignal = input('请输入迭代初值矩阵:oiginal = ') n = input('请输入迭代次数:n = ') m = input('请选择右端项(1-b1;2-b2):m = ') A = [6,2,-1;1,4,-2;-3,1,4]; b1 = [-3,2,4]'; %转置为列向量 b2 = [100,-200,345]'; d = diag(A); for i = 1:3 for j = 1:3 if i == j dj(i,j) = d(i); else dj(i,j) = 0; end end end dj switch m case 1 t = b1; case 2 t = b2; end B = eye(length(A)) - inv(dj)*A; g = inv(dj)*t; x{1} = orignal'; %从矩阵A 中提取对角元素 %写成对角阵 %选择右端项 %计算迭代矩阵B %计算矩阵g for k = 1:n X{1} = B*x{1}+g; x{1} = X{1}; diedai = A*X{1} %使用迭代得到的值计算左端项 error = diedai - b1 %左端项与右端项的差距,作为迭代误差指标 if error == 0 step = k; %记录迭代步数 break; end end step x1 = X{1}(1) %方程组的根x1 x2 = X{1}(2) %方程组的根x2 x3 = X{1}(3) %方程组的根x3 (2)高斯-赛德尔迭代法 clc; clear; orignal = input('请输入迭代初值矩阵:oiginal = ') n = input('请输入迭代次数:n = ') m = input('请选择右端项(1-b1;2-b2):m = ') A = [6,2,-1;1,4,-2;-3,1,4]; b1 = [-3,2,4]'; %转置为列向量 b2 = [100,-200,345]'; d = diag(A); %从矩阵A 中提取对角元素 for i = 1:3 %写成对角阵 for j = 1:3 if i == j D(i,j) = d(i); else D(i,j) = 0; end end end D switch m case 1 t = b1; case 2 t = b2; end %计算L 矩阵 l = tril(A); L = -(l - D); %计算U 矩阵 u = triu(A); U = -(u - D); B = inv(D - L)*U g = inv(D - L)*t x{1} = orignal'; for k = 1:n X{1} = B*x{1}+g; %选择右端项 %计算迭代矩阵B %计算矩阵g x{1} = X{1}; diedai = A*X{1} %使用迭代得到的值计算左端项 error = diedai - b1 %左端项与右端项的差距,作为迭代误差指标 if error == 0 step = k; %记录迭代步数 break; end end step x1 = X{1}(1) %方程组的根x1 x2 = X{1}(2) %方程组的根x2 x3 = X{1}(3) %方程组的根x3 4.3 题3程序代码 a=1; b=3; h=0.5; y(1)=-2; x(1)=a; n=(b-a)/h+1; yy(1)=-2; for i=2:n k1=(y(i-1)^2+y(i-1))/x(i-1); k2=((y(i-1)+h*k1/2)^2+(y(i-1)+h*k1/2))/(x(i-1)+h/2); k3=((y(i-1)+h*k2/2)^2+(y(i-1)+h*k2/2))/(x(i-1)+h/2); k4=((y(i-1)+h*k3)^2+(y(i-1)+h*k3))/(x(i-1)+h); y(i)=y(i-1)+h*(k1+2*k2+2*k3+k4)/6; % 四阶Runge-Kutta 公式解 x(i)=x(i-1)+h; %有解区间的值 yy(i)=-x(i)/(x(i)-1/2); %解析解 s(i)=abs(y(i)-yy(i)); %误差项 end [x' y' yy' s'] 数值分析--上机实习报告 学 号: 12011345 姓 名: 高 剑 专 业: 建筑与土木工程 联系电话: 13550365902 任课老师: 谢玲红 二零一二年十二月 数值分析上机报告 目 录 目 录 ........................................................................................................................................... 1 第一节 序言 ................................................................................................................................. 2 第二节 插值法求函数 ............................................................................................................... 3 2.1 问题.................................................................................................................................... 3 2.2 计算原理和方法简介 ...................................................................................................... 4 2.3 计算结果............................................................................................................................ 5 2.4 问题解答............................................................................................................................ 7 第三节 松弛因子对SOR法收敛速度的影响 ........................................................................... 8 3.1 问题.................................................................................................................................... 8 3.2 松弛法简介 ................................................................................................................... 8 3.2.1 松弛法基础 ............................................................................................................. 8 3.2.3松弛法收敛条件 ...................................................................................................... 8 3.3 计算结果 .......................................................................................................................... 9 3.4 结果对比及分析 ............................................................................................................ 10 第四节 用Runge-Kutta 4阶算法 ............................................................................................... 11 4.1 Runge-Kutta 4阶算法求解的题目 .................................................................................. 11 4.2 Runge-Kutta 4阶算法求解过程及结果 .......................................................................... 11 4.3 总结.................................................................................................................................. 13 第五节 附录 ............................................................................................................................... 13 6.1 laglange插值程序 ............................................................................................................ 13 6.2 SOR迭代法程序 .............................................................................................................. 14 6.3四阶 Runge-Kutta 算法程序代码: .............................................................................. 15 1 数值分析上机报告 第一节 序言 数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科,是数学的一个分支,它以数字计算机求解数学问题的理论和方法为研究对象。为计算数学的主体部分。运用数值分析解决问题的过程:实际问题?数学模型?数值计算方法?程序设计?上机计算求出结果数值分析这门学科有如下特点:1面向计算机;2有可靠的理论分析;3要有好的计算复杂性;4要有数值实验;5要对算法进行误差分析。 本学期开设了数值分析课程,该课程讲授了数值分析绪论、非线性方程的求解、插值法、函数逼近与曲线拟合、线性方程组的直接接法、线性方程组的迭代法、数值积分和数值微分、常微分方程初值问题的数值解法等内容。为我们解决实际数学问题提供了理论基础,很多问题的求解人工计算量太大甚至无法操作,必须借助计算机运算。所以学好数值分析的关键是要加强上机操作,充分利用计算机程序语言实现数值分析的算法。 本次上机实习采用C++语言编程解决问题。C++语言是近年来在国内外得到迅速推广应用的一种计算机语言,主要有以下优点:C++是严谨、精确、数理化的。标准定义很细致,如果你对标准深入了解,非“固有不可移植特性”代码的运行时状态是可以通过静态分析确定的;而就其本身的语法思路,是层次分明的、呼应的,有着丰富的逻辑演绎色彩的;其语言结构而言,代码绝大多数是显式的、明确的。C++是活着的。它的标准在不断更新,变得更加严谨、易用,却保持着系统性。C++是支持底层操作的。虽然许多底层操作的具体效果会随机器变化有所差异(也就是固有不可移植特性),但如果必要,你依然可以使用它们编写程序、改善效率。C++的标准是庞大、层级分明的。如果你能熟练掌握全部C++标准(语法、标准库),你的代码编写工作会极其精炼、轻松;而如果你只是用它来实验,不愿学习过多细节,它的基础特性又可以认为是一个“完备集”。C++兼容C语言的几乎所有特征(部分极少的不严谨、存在歧义的特性被去除),你依旧采用C语言的思路,却可以轻易加入一个C++工程的开发。C++拥有面向对象特性,“模板”等高度抽象化特性使得大型工程可以通过它进行整体管理。它考虑了兼容性,为连接其他语言的对象文件留有接口。因此,C++更具有发展潜力。 2 数值分析上机报告 本次上机实习中使用计算机求解了线性方程组问题。每一题中都介绍了相关 理论、计算方法以及结果分析,以加强对理论知识的掌握。 第二节 插值法求函数 2.1 问题 已知:a=-5,b=5, 以下是某函数 f(x)的一些点(xk,yk), 其中 xk=a+0.1(k-1) ,k=1,..,101 xk=a+0.1k, 请用插值类方法给出函数f(x)的一个解决方案和具体结果。并通过实验考虑下 列问题 (1)Ln(x)的次数n越高,逼近f(x)的程度越好, (2)高次插值收敛性如何, (3)如何选择等距插值多项式次数 , (4)若要精度增高,你有什么想法, 比如一定用插值吗, (5)逼近某个函数不用插值方式,有何变通之举, (6)函数之间的误差如何度量,逼近的标准又是什么, (7)如何比较好的使用插值多项式呢, x = k Columns 1 through 7 -5.0000 -4.9000 -4.8000 -4.7000 -4.6000 -4.5000 -4.4000 Columns 8 through 14 -4.3000 -4.2000 -4.1000 -4.0000 -3.9000 -3.8000 -3.7000 Columns 15 through 21 -3.6000 -3.5000 -3.4000 -3.3000 -3.2000 -3.1000 -3.0000 Columns 22 through 28 -2.9000 -2.8000 -2.7000 -2.6000 -2.5000 -2.4000 -2.3000 Columns 29 through 35 -2.2000 -2.1000 -2.0000 -1.9000 -1.8000 -1.7000 -1.6000 Columns 36 through 42 -1.5000 -1.4000 -1.3000 -1.2000 -1.1000 -1.0000 -0.9000 Columns 43 through 49 -0.8000 -0.7000 -0.6000 -0.5000 -0.4000 -0.3000 -0.2000 Columns 50 through 56 -0.1000 0 0.1000 0.2000 0.3000 0.4000 0.5000 Columns 57 through 63 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 Columns 64 through 70 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 Columns 71 through 77 2.0000 2.1000 2.2000 2.3000 2.4000 2.5000 2.6000 Columns 78 through 84 3 数值分析上机报告 2.7000 2.8000 2.9000 3.0000 3.1000 3.2000 3.3000 Columns 85 through 91 3.4000 3.5000 3.6000 3.7000 3.8000 3.9000 4.0000 Columns 92 through 98 4.1000 4.2000 4.3000 4.4000 4.5000 4.6000 4.7000 Columns 99 through 101 4.8000 4.9000 5.0000 y(x)=y= kk Columns 1 through 7 25.0000 24.0100 23.0400 22.0900 21.1600 20.2500 19.3600 Columns 8 through 14 18.4900 17.6400 16.8100 16.0000 15.2100 14.4400 13.6900 Columns 15 through 21 12.9600 12.2500 11.5600 10.8899 10.2397 9.6093 8.9991 Columns 22 through 28 8.4092 7.8405 7.2941 6.7705 6.2693 5.7866 5.3144 Columns 29 through 35 4.8403 4.3522 3.8463 3.3402 2.8832 2.5554 2.4475 Columns 36 through 42 2.6154 3.0219 3.4920 3.7149 3.3232 2.0435 -0.1277 Columns 43 through 49 -2.8066 -5.2470 -6.5469 -5.9893 -3.3862 0.7365 5.2312 Columns 50 through 56 8.6985 10.0000 8.6985 5.2312 0.7365 -3.3862 -5.9893 Columns 57 through 63 -6.5469 -5.2470 -2.8066 -0.1277 2.0435 3.3232 3.7149 Columns 64 through 70 3.4920 3.0219 2.6154 2.4475 2.5554 2.8832 3.3402 Columns 71 through 77 3.8463 4.3522 4.8403 5.3144 5.7866 6.2693 6.7705 Columns 78 through 84 7.2941 7.8405 8.4092 8.9991 9.6093 10.2397 10.8899 Columns 85 through 91 11.5600 12.2500 12.9600 13.6900 14.4400 15.2100 16.0000 Columns 92 through 98 16.8100 17.6400 18.4900 19.3600 20.2500 21.1600 22.0900 Columns 99 through 101 23.0400 24.0100 25.0000 2.2 计算原理和方法简介 本题需要为了用插值法求出函数。首先绘出101个点的函数曲线图,通过 matlab运用拉格朗日插值计算出用不同次数的插值求出的函数曲线与前面的曲 4 数值分析上机报告 线图比较,从而得出问题的解答。 下面简单介绍Lagrange插值多项式,拉格朗日插值法是基于基函数的插值方法,插值多项式可表示为: n yxfxlx()()(),,jjj,0 lx()j其中称为j次基函数: ()()()()()xxxxxxxxxx,,,,,,,,,,,0111jjn,,lx(),jk()()()()()xxxxxxxxxx,,,,,,,,,,,jjjjjjjn0111,, 很容易验证它满足插值条件。 2.3 计算结果 不同次数多项式等距插值比较 2次 f = 0.6*x^2 + 0.0000000000000204281*x + 10.0 3次 f =- 0.00306668*x^3 + 1.00992*x^2 + 0.0766671*x - 0.248002 5 数值分析上机报告 4次 f = 0.0614143*x^4 + 0.00605248*x^3 - 0.935358*x^2 - 0.151312*x + 10.0 5次 f = - 0.00263324*x^5 + 0.00767868*x^4 + 0.091116*x^3 + 0.734223*x^2 - 0.632127*x + 1.84524 6次 f = - 0.0135359*x^6 - 0.00234626*x^5 + 0.525821*x^4 + 0.0733606*x^3 - 4.08561*x^2 - 0.3676*x + 10.0 6 数值分析上机报告 10次 由于10次以上多项式没有使用价值,所以函数不列出。 多项式插值次数越高,对节点数据越敏感。数据越不稳定。 2.4 问题解答 (1)Ln(x)的次数n越高,逼近f(x)的程度越好, 错误!高次插值有病态性,比如龙格现象,在本题示例中可以看出。 (2)高次插值收敛性如何, 较差!很难达到一致收,因为龙格现象,两边会有很大的突起。 (3)选择等距插值多项式次数, 根据选点数目,n个点进行插值,则多项式次数为n-1。 (4)若要精度增高,你有什么想法,比如一定用插值吗, 在本题中用分段低次插值和拟合法效果都比较好,可以提高精度,避免病态现象。 (5)逼近某个函数不用插值方式,有何变通之举, 拟合法。 (6)函数之间的误差如何度量,逼近的标准又是什么, 是具体情况使用各种范数,最近看的函数型数据分析里面有粗糙惩罚法,min MSE+拟合函数的粗糙度(一般是二阶导平方求积分) 7 数值分析上机报告 (7)如何比较好的使用插值多项式呢 分段低次插值,取特征点进行插值。 第三节 松弛因子对SOR法收敛速度的影响 3.1 问题 松弛因子对SOR法收敛速度的影响,用SOR法求解方程组Ax=b,其中: ,41-3,,,,,,,,1,41-2,,,, ,,,,...-2,,,,A= B,....,,,, ,,,,1,41-2,,,,,,,,1,4-3,,,, 要求程序中不存系数矩阵A,分别对不同的阶数取w=1.1, 1.2, ...,1.9进行迭 (K)(K)(K-1),6XX,X10代,记录近似解达到||||<时所用的迭代次数k,观察松弛因子对收敛速度的影响,并观察当w,0或w,2会有什么影响, 3.2="" 松弛法简介="">时所用的迭代次数k,观察松弛因子对收敛速度的影响,并观察当w,0或w,2会有什么影响,> 3.2.1 松弛法基础 高斯-赛德尔方法的迭代矩阵形式可以写成: ,,,,,,k,1k,1k,1,,k ,,x,x,D[b,Lx,D,Ux] 松弛法相对于高斯-赛德尔方法只是添加了一个松弛系数,所以松弛法的, ,,,,,,k,1k,1k,1,,k矩阵形式为: ,,x,x,,D[b,Lx,D,Ux] ,1,1,,k,,k,1整理得:,,,,,, x,D,,LD,,D,,Ux,D,,L,b 所以,松弛法的迭代矩阵为: ,1,1,,,, ,, M,D,,LD,,D,,UT,D,,L,b3.2.3松弛法收敛条件 ,,0xAAx,b设有方程组,为非奇异阵,且a,0,对于任何初始向量及实数,,ii 8 数值分析上机报告 ,1用松弛法收敛的充要条件是(其中)。 ,,,,M,D,,LD,,D,,U,,,M,1 松弛法收敛的必要条件是。 ,,,,0,2 3.3 计算结果 当矩阵A为二阶: W w=1.1 w=1.2 w=1.3 w=1.4 w=1.5 w=1.6 w=1.7 w=1.8 w=1.9 次数 k 8 10 13 17 23 30 43 65 140 当矩阵A为三阶: W w=1.1 w=1.2 w=1.3 w=1.4 w=1.5 w=1.6 w=1.7 w=1.8 w=1.9 次数 k 8 11 14 17 23 29 43 66 140 当矩阵A为四阶: W w=1.1 w=1.2 w=1.3 w=1.4 w=1.5 w=1.6 w=1.7 w=1.8 w=1.9 次数 k 8 11 14 18 22 31 42 66 138 当矩阵A为五阶: W w=1.1 w=1.2 w=1.3 w=1.4 w=1.5 w=1.6 w=1.7 w=1.8 w=1.9 次数 k 8 11 14 18 24 30 42 66 140 当矩阵A为六阶: W w=1.1 w=1.2 w=1.3 w=1.4 w=1.5 w=1.6 w=1.7 w=1.8 w=1.9 次数 k 9 12 15 18 24 30 43 67 142 当矩阵A为七阶: W w=1.1 w=1.2 w=1.3 w=1.4 w=1.5 w=1.6 w=1.7 w=1.8 w=1.9 次数 k 10 13 15 18 23 32 44 68 139 9 数值分析上机报告 当矩阵A为八阶: W w=1.1 w=1.2 w=1.3 w=1.4 w=1.5 w=1.6 w=1.7 w=1.8 w=1.9 次数 k 11 13 16 19 22 33 42 67 140 当矩阵A为九阶: W w=1.1 w=1.2 w=1.3 w=1.4 w=1.5 w=1.6 w=1.7 w=1.8 w=1.9 次数 k 12 14 17 19 23 31 43 67 142 当矩阵A为十阶: W w=1.1 w=1.2 w=1.3 w=1.4 w=1.5 w=1.6 w=1.7 w=1.8 w=1.9 次数 k 12 15 17 20 24 30 45 68 140 当矩阵A为十五阶: W w=1.1 w=1.2 w=1.3 w=1.4 w=1.5 w=1.6 w=1.7 w=1.8 w=1.9 次数 k 12 15 18 24 28 33 43 67 139 当矩阵A为二十五阶: W w=1.1 w=1.2 w=1.3 w=1.4 w=1.5 w=1.6 w=1.7 w=1.8 w=1.9 次数 k 12 15 18 23 30 42 50 68 152 3.4 结果对比及分析 1. 对于给定的系数矩阵A及矩阵B,对于不同阶数的系数矩阵A,随着松弛因子w在1.1—1.9之间连续递增时,由上观察得迭代次数越来越大,收敛速度越来越慢。随着矩阵A阶数的增大,迭代次数相对稳定。 ,2,2. 对于当w0或w时,分别取w=-0.5,2.5运行观察得: 迭代过程如同“死循环”,松弛法不收敛,从而证明了松弛法收敛的必要条件是: ,2,,,0w2的正确性,即当w0或w时,松弛法不收敛。 10 数值分析上机报告 第四节 用Runge-Kutta 4阶算法 4.1 Runge-Kutta 4阶算法求解的题目 /用Runge-Kutta 4阶算法对初值问题y=-20*y,y(0)=1按不同步长求解,用于观察稳定区间的作用,推荐两种步长h=0.1,0.2。 -20x注:此方程的精确解为:y=e 4.2 Runge-Kutta 4阶算法求解过程及结果 在 MATLAB 中按照标准四阶 Runge-Kutta 的算法编写程序代码,具体代码见附录。并运行程序 解:首先探讨不同步长下算法稳定的问题,取不同步长得到求解的微分方程的解同精确解的误差分别如表 1、表 2、表 3、表 4、表 5。 11 数值分析上机报告 从上面得表中可以看到,在 h=0.0001,0.001,0.01,0.1 时,绝对误差在逐步变小,因此在这种情况下认为此方法是稳定的。当 h=0.2 时,绝对误差逐渐扩大,因此在 h=0.2 时,算法不稳定。 现在来探讨不同 h 下面的精度问题,取不同的步长求 y(0.1)的值得表 6. 12 数值分析上机报告 从表 6 中可以看到步长越小,标准四阶 Runge-Kutta 算法得到的值精度越高。计算所需的步骤也越长。 4.3 总结 通过这次的编程计算,使我学到了很多东西。在完成实验内容后,使我既学习计算方法又提高计算能力的目的,还使我切身体会书本内容之精妙所在,在编程中间使我体会到了很多乐趣。数值分析作为工程计算和科学计算连接计算机的一门课程,因此他的实用性是比较强的,因此在我们掌握了一定的理论基础的情况下,就需要我们学以致用,能解决实际问题。因此这次的实验内容能使我们学到很多的东西,对提高解决实际问题的能力有了一定的提高。并对理论中的一些现象,通过上机的分析来论证了理论分析的正确性,相信这次的编程计算对我以后的学习和解决工程问题会有很大的帮助的。 第五节 附录 6.1 laglange插值程序 function laglange(xk,yk) x0=1; syms x w=1; for i=1:n w=w*(x-X(i)); end a=diff(w,x); f=0; for i=1:n b=w/(x-X(i)); l=b/(subs(a,x,X(i))); f=f+l*Y(i); end 13 数值分析上机报告 f=expand(f) f=vpa(f,6) %subs(f,x,x0) sprintf('f(%f)=%f',x0,subs(f,x,x0)) 6.2 SOR迭代法程序 # define N 3 # include # include { FILE *fp; fp=fopen("sor.txt","a"); int i,j,k,g; double p, q,s,t,eps ; double a[N][N], x[N],b[N],w[9]; for(i=0;i for(j=0;j scanf("%lf",&a[i][j]); for(i=0;i scanf("%lf",&b[i]); scanf("%lf",&eps); for(i=0;i<9;i++)>9;i++)> scanf("%lf",&w[i]); for(g=0;g<9;g++)>9;g++)> { for(i=0;i x[i]=0.0; p=eps+1.0; for(k=1;p>=eps;k++) { fprintf(fp,"k=%d\n",k); p=0.0; 14 数值分析上机报告 for(i=0;i { t=x[i]; s=0.0; for(j=0;j if(i!=j) s=s+a[i][j]* x[j]; x[i]=(1-w[g])*x[i]+(w[g]*(b[i]-s))/a[i][i]; fprintf(fp,"%15.10f",x[i]); if(i==N-1) fprintf(fp,"\n"); q=fabs(x[i]-t); if(q>p) p=q; } } } fclose(fp); } 6.3四阶 Runge-Kutta 算法程序代码: h=[0.0001,0.001,0.01,0.1] k=1 for xk=0.1:0.1:1 %取不同的x进行计算 for i=1:1:4 %取不同的步长计算每个x hj=h(i) yk=1 for xkp=hj:hj:xk %四阶标准Runge-Kutta迭代算法 k1=hj*(-20*yk); format long k1; k2=hj*(-20*(yk+k1*0.5)); k3=hj*(-20*(yk+k2*0.5)); k4=hj*(-20*(yk+k3)); 15 数值分析上机报告 ykm=yk+(k1+2*k2+2*k3+k4)*(0.16666666666667); yk=ykm; end yk=ykm; yeal=exp(-20*xkp) ; %求出精确值 ck=abs(yk-yeal); %精确值和计算值之间的绝对误差 ykeal(k)=yeal; ckk(k)=ck; yxk(k)=yk; xkk(k)=xk; k=k+1; end end D=[yxk' ykeal' ckk']; %将迭代得到的结果及误差存到矩阵D中 for i=1:1:10 %将同一个步长下不同X得到的结果集中到一个矩阵 for j=1:1:3 Ah1(i,j)=D(4*i-3,j); Ah2(i,j)=D(4*i-2,j); Ah3(i,j)=D(4*i-1,j); Ah4(i,j)=D(4*i,j); end Axk(i)=xkk(4*i-3); end Ch1=[Axk' Ah1] Ch2=[Axk' Ah2] Ch3=[Axk' Ah3] Ch4=[Axk' Ah4] Cheng=[Ch1(1,[2 3 4]);Ch2(1,[2 3 4]); Ch3(1,[2 3 4]);Ch4(1,[2 3 4])] 16 数值分析上机报告 数值分析上机报告 篇一: 数值分析实验报告 (一)(完整) 数值分析实验报告 1 2 3 4 5 篇二: 哈工大-数值分析上机实验报告 实验报告一 题目: 非线性方程求解 摘要: 非线性方程的解析解通常很难给出,因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Netn法及改进的Netn法。 前言: (目的和意义) 掌握二分法与Netn法的基本原理和应用。 数学原理: 对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法: 二分法和Netn法。 对于二分法,其数学实质就是说对于给定的待求解的方程f(x),其在[a,b]上连续,f(a)f(b) 0,且f(x)在[a,b]内仅有一个实根x*,取区间中点c,若,则c恰为其根,否则根据f(a)f(c) 0是否成立判断根在区间[a,c]和[c,b]中的哪一个,从而得出新区间,仍称为[a,b]。重复运行计算,直至满足精度为止。这就是二分法的计算思想。 Netn法通常预先要给出一个猜测初值x0,然后根据其迭代公式 xk?1?xk?f(xk) f(xk) 产生逼近解x*的迭代数列{xk},这就是Netn法的思想。当x0接近x*时收敛很快,但是当x0选择不好时,可能会发散,因此初值的选取很重要。另外,若将该迭代公式改进为 xk?1?xk?rf(xk) f (xk) 其中r为要求的方程的根的重数,这就是改进的Netn法,当求解已知重数的方程的根时,在同种条件下其收敛速度要比Netn法快的多。 程序设计: 本实验采用Matlab的M文件编写。其中待求解的方程写成functin的方式,如下 functin y=f(x); y=-x*x-sin(x); 写成如上 形式即可,下面给出主程序。 二分法源程序: clear %%%给定求解区间 b= 1.5; a=0; %%%误差 R=1; k=0;%迭代次数初值 hile (R 5e-6) ; c=(a+b)/2; if f12(a)*f12(c) a=c; else b=c; end R=b-a;%求出误差 k=k+1; end x=c%给出解 Netn法及改进的Netn法源程序: clear %%%% 输入函数 f=input( 请输入需要求解函数 , s ) %%%求解f(x)的导数 df=diff(f); %%%改进常数或重根数 miu=2; %%%初始值x0 x0=input( input initial value x0 k=0;%迭代次数 max=100;%最大迭代次数 R=eval(subs(f, x0 , x %求解f(x0),以 确定初值x0时否就是解 hile (abs(R) 1e-8) x1=x0-miu*eval(subs(f, x0 , x ))/eval(subs(df, x0 , x R=x1-x0; x0=x1; k=k+1; if (eval(subs(f, x0 , x )) 1e-10); break end if k %如果迭代次数大于给定值,认为迭代不收敛,重新输入初值 ss=input( maybe result is errr,chse a ne x0,y/n? , s if strcmp(ss, y ) x0=input( input initial value x0 k=0; else break end end end k;%给出迭代次数 x=x0;%给出解 结果分析和讨论: x2 ?0在[1,2]内的根。(??5*10?6,下同) 1. 用二分法计算方程sinx?2 计算结果为 x= 1.40441513061523; f(x)= - 3.797205105904311e-007; k=18; 由f(x)知结果满足要求,但 迭代次数比较多,方法收敛速度比较慢。 2. 用二分法计算方程x3?x?1?0在[1, 1.5]内的根。 计算结果为 x= 1.32471847534180; f(x)= 2.209494846194815e-006; k=17; 由f(x)知结果满足要求,但 迭代次数还是比较多。 3. 用Netn法求解下列方程 a) xex?1?0 x0=0.5; 计算结果为 x= 0.56714329040978; f(x)= 2.220446049250313e-016; k=4; 由f(x)知结果满足要求,而且又迭代次数只有4次看出收敛速度很快。 b) x3?x?1?0 x0=1; c) (x?1)2(2x?1)?0 x0=0.45, x0=0.65; 当x0=0.45时,计算结果为 x= 0.49999999999983; f(x)= - 8.362754932994584e-014; k=4; 由f(x)知结果满足要求,而且又迭代次数只有4次看出收敛速度很快,实际上该方程确实有真解x=0.5。 当x0=0.65时,计算结果为 x= 0.50000000000000; f(x)=0; k=9; 由f(x)知结果满足要求,实际上该方程确实有真解x=0.5,但迭代次数增多,实际上当取x0〉0.68时,x?1,就变成了方程的另一个解,这说明Netn法收敛与初值很有关系,有的时候甚至可能不收敛。 4. 用改进的Netn法求解,有2重根,取??2 (x?1)2(2x?1)?0 x0=0.55;并与 3.中的c)比较结果。 当x0=0.55时,程序死循环,无法计算,也就是说不收敛。改?? 1.5时,结果收敛为 x=0.50000087704286; f(x)= 4.385198907621127e-007; k=16; 显然这个结果不是很好,而且也不是收敛至方程的2重根上。 当x0=0.85时,结果收敛为 x= 1.00000000000489; f(x)= 2.394337647718737e-023; k=4; 这次达到了预期的结果,这说明初值的选取很重要,直接关系到方法的收敛性,实际上直接用Netn法,在给定同样的条件和精度要求下,可得其迭代次数k=15,这说明改进后的Netn法法速度确实比较快。 结论: 对于二分法,只要能够保证在给定的区间内有根,使能够收敛的,当时收敛的速度和给定的区间有关,二且总体上来说速度比较慢。Netn法,收敛速度要比二分法快,但是最终其收敛的结果与初值的选取有关,初值不同,收敛的结果也可能不一样,也就是结果可能不时预期需要得结果。改进的Netn法求解重根问题时,如果初值不当,可能会不收敛,这一点非常重要,当然初值合适,相同情况下其速度要比Netn法快得多。 实验报告二 题目: Gauss列主元消去法 摘要: 求解线性方程组的方法很多,主要分为直接法和间接法。本实验运用直接法的Guass消去法,并采用选主元的方法对方程组进行求解。 前言: (目的和意义) 1. 学习Gauss消去法的原理。 2. 了解列主元的意义。 3. 确定什么时候系数阵要选主元 数学原理: (k?1)由于一般线性方程在使用Gauss消去法求解时,从求解的过程中可以看到,若akk=0, (k?1)则必须进行行交换,才能使消去过程进行下去。有的时候即使akk?0,但是其绝对值非 常小,由于机器舍入误差的影响,消去过程也会出现不稳定得现象,导致结果不正确。因此有必要进行列主元技术,以最大可能的消除这种现象。这一技术要寻找行r,使得 (k?1)(k?1)|ark|?maxaik i?k (k?1)并将第r行和第k行的元素进行交换,以使得当前的akk的数值比0要大的多。这种列主 元的消去法的主要步骤如下: 1. 消元过程 对k=1,2,…,n-1,进行如下步骤。 1) 选主元,记 |ark|?aik i?k 若|ark|很小,这说明方程的系数矩阵严重病态,给出警告,提示结果可能不对。 2) 交换增广阵A的r,k两行的元素。 arj?akj (j=k,…,n+1) 3) 计算消元 aij?aij?aikakj/akk (i=k+1,…,n; j=k+1,……,n+1) 2. 回代过程 对k= n, n-1,…,1,进行如下计算 篇三: 数值分析上机实验报告-李宝君 《 数 值 分 析 实 验 报 告 》 指导老师: 代新敏 姓 名: 李宝君 学 号: 201X020917 专 业: 机械制造及其自动化 院 系: 机械工程学院 贵州大学201X级研究生 第一题 一 程序说明: 1 Husehlder其基本思想是: 利用初等反射阵H?E?2uuT,将矩阵的每一行向量 变换成所需方向的行向量,从而最终得到想要的三对角阵。它的乘法运算次数仅是Givenr方法的一半,且只需要作n?2次开方运算。 归纳起来,对换矩阵三对角化的算法步骤为: (1) 令A0?A,aij?aij,已知Ar?1,即Ar?1?(aij)。 (2) sr?( (1) (r) i?r?1 ?(a n (r)21/2ir )) T (r)(r)(r)(r)) ?? (3) ?r?sr2?ar(r?1,rsr., ur??0,?0,ar?1,r?sign(ar?1,r)sr,ar?2,r?an,r?. (4) yr?Ar?1ur/?r。 (5) kr? 1T uryr/?r。 2 (6) qr?yr?krur TT (7) Ar?Ar?1?(qrur?urqr), r?1,2,?n?2。 2 松弛法简称SR法,其基本思想是在GS方法已求出x(m),x(m?1)的基础上,经过 重新组合而得到新的序列,而此新序列使收敛速度加快。其算法如下: i?1 n x (m) i ?(1??)x (m?1)i ??(?bx j?1 (m)ijj??bijx(jm?1)?gi) j?i?1 当? 1时称为超松弛法。 3 用列主元素消去法求解BX=b: n阶方程组的系数矩阵为: a11 a12…a1n b1 a21 a22…a2n b2 an1 an2?ann b Gauss消去法的算法为: ij=aij/ajj (ajj!=0) j=1,2?n, i=j+1,j+2?n ? ij=aij-lik-1ak-1 i,j=k,k+1?n, k=2,3?n i=bi-lik-1bk-1 i=k,k+1?n, k=2,3?n ? xi=(bi-?aijxj)/aii i=n,n-1?1,j=i+1,i+2,?n 列主元素消去法是在Gauss消去法的基础上选主元,选取绝对值最大(或尽量大)的元素为主元,使lij绝对 值很小。 二 计算程序: #include stdafx.h #include stdi.h #include math.h #define MAX 9 vid SR(duble (*p)[MAX],int n); //定义一个函数实现SR法 vid L(duble (*p)[MAX],int n); //定义一个函数实现列主消元 vid sap(duble *a,duble*b) //定义一个函数实现交换功能 { } vid main { duble A[MAX][MAX]= { {1 2.38412, 2.115237,- 1.061074, 1.112336,-0.113584,0.718719, 1.742382, 3.067813,- 2.03duble t; t=*a;*a=*b;*b=t; return; 1743}, { 2.115237,1 9.141823,- 3.125432,- 1.015, 2.189736, 1.563849,-0.784165, 1.118, 3.123124}, {- 1.061074,- 3.125432,1 5.567914, 3.123848, 2.031454, 1.836742,- 1.056781,0.336993,- 1.010103}, { 1.112336,- 1.015, 3.123848,2 7.108437, 4.101011,- 3.741856, 2.101023,-0.71828,-0.037585}, 4317}, 8417}, {-0.113584, 2.189736, 2.031454, 4.101011,1 9.897919,0.431637,- 3.111223, 2.121314, 1.78{0.718719, 1.563849, 1.836742,- 3.741856,0.431637, 9.789365,-0.103458,- 1.103456,0.23 { 1.742382,-0.784165,- 1.056781, 2.101023,- 3.111223,-0.103458,1 4.7138465, 3.123789,- 2.213474}, 6782}, { 3.067813, 1.118,0.336993,-0.71828, 2.121314,- 1.103456, 3.123789,30.7193344, 4.44{- 2.031743, 3.123124,- 1.010103,-0.037585, 1.784317,0.238417,- 2.213474, 4.446782,40. 00001} }; duble U[MAX],Y[MAX]; duble ku,s,a; int n=MAX,i,j,k,m,sign; fr(k=0;k k++) //转化为三对角阵 { } fr(m=k+1;m m++) s+=p(A[m][k],2); //计算s sign=(A[k+1][k] 0)?-1:1; a=s+p(s,0.5)*fabs(A[k+1][k]); //计算a fr(m=0;m m++) //计算U[] { if (m =k) U[m]=0; else if(m==k+1) U[m]=A[k+1][k]+sign*p(s,0.5); else U[m]=A[m][k]; } fr(i=0;i i++) //计算Y[] { Y[i]=0.0; fr(j=0;j j++) Y[i]+=A[i][j]*U[j]; Y[i]=Y[i]/a; } ku=0.0; fr(m=0;m m++) ku+=U[m]*Y[m]; ku=ku/(2*a); fr(i=0;i i++) //转化A为B fr(j=0;j j++) A[i][j]=A[i][j]-Y[i]*U[j]-U[i]*Y[j]+2*ku*U[i]*U[j]; } printf( 转换后的B矩阵:\n fr(i=0;i i++) { fr(j=0;j j++) printf( % 9.6lf ,A[i][j]); printf( \n } SR(A,n); L(A,n); vid SR(duble (*p)[MAX],int n) { duble B[MAX][MAX]; duble x[MAX]={0},b[MAX]={ 2.1874369,3 3.992318,-2 5.173417,0.84671695, 1.784317,-8 6.613, 1.1101230, 4.719345,- 5.6784392}; int i,j,k; } fr(i=0;i i++) { b[i]=b[i]/(*(*(p+i)+i)); fr(j=0;j j++) } fr(k=0;k k++) { fr(i=0;i i++) { t=0; fr(j=0;j j++) t+=B[i][j]*x[j]; fr(j=i+1;j j++) t+=B[i][j]*x[j]; B[i][j]=(i==j)?0:-(*(*(p+i)+j))/(*(*(p+i)+i)); x[i]=(1-)*x[i]+*(t+b[i]); } printf( SR法求解X:\n fr(k=0;k k++) { printf( x[%d]=%.6lf\n ,k,x[k]); } return; vid L(duble (*p)[MAX],int n) { duble x[MAX],L[MAX],b[MAX]={ 2.1874369,3 3.992318,-2 5.173417,0.84671695, 1.784317,-8 6.613, 1.1101230, 4.719345,- 5.6784392}; duble B[MAX][MAX]; duble max; int i,k,m,j; fr(i=0;i i++) fr(j=0;j j++) B[i][j]=(*(*(p+i)+j)); fr(i=0;i i++) { max=0; fr(j=i;j j++) if(fabs(B[j][i]) =max) //找出每 列最大值 { max=fabs(B[j][i]);篇四: 数值分析实验报告(包含源程序) 课程实验报告 课程实验报告 数值分析上机实习报告 专 业:土木工程 班 级: 学 号: 姓 名: 指导老师: 联系电话: 2015.12.12 序言 随着本学期逐渐接近尾声, 我也逐渐掌握了数值分析的一些基本理论. 本次上机作业是理论与实践的结合. 本次作业使用了matlab 与c++两种语言. 其中matlab 具有编程效率高,用户使用方 便,方便的绘图功能的优点。而c++是一种基本的编程语言,在实际的工程中也有广泛的应用。 本次作业根据题目的特点,结合两种语言各自的优势,采用了不同的方法。其中牛顿法,Steffensen 加速法采用了c 语言。插值与多项式拟合使用了两种语言。Ru n ge-Kutt a 4阶算法仅使用了matlab 编程。 本次作业注重问题的计算过程,分析总结,及编程。由于所涉及原理课本均有详细陈述,在此不再赘述。 第一题 ...................................................................................................................................... 3 1.1题目 . ............................................................................................................................ 3 1.2计算过程和结果 . ........................................................................................................ 3 1.3结果分析 . .................................................................................................................... 3 第二题 ...................................................................................................................................... 4 2.1题目 . ............................................................................................................................ 4 2.2计算过程和结果 . ........................................................................................................ 4 2.3结果分析 . .................................................................................................................... 8 第三题 ...................................................................................................................................... 8 3.1题目 . ............................................................................................................................ 8 3.2问题求解及过程 . ........................................................................................................ 8 3.3结果分析 . .................................................................................................................... 9 总结 ........................................................................................................................................ 10 附 件................................................................................................................................. 11 第一题............................................................................................................................. 11 1.1.1第一问牛顿法 . ............................................................................................... 11 1.1.2第一问牛顿-Steffensen 法 . ........................................................................... 11 1.2.1第二问牛顿法 . ............................................................................................... 12 1.2.2第二问牛顿-Steffensen 法 . ........................................................................... 13 第二题............................................................................................................................. 14 2.1.1最小二乘法求解 ............................................................................................ 14 2.2.1拉格朗日差值多项式拟合 ............................................................................ 15 2.2.2牛顿插值 ........................................................................................................ 15 第三题 . ............................................................................................................................ 17 3.1.1Runge-Kutta 4阶算法 .................................................................................. 17 第一题 1.1题目 分别用牛顿法,及基于牛顿算法下的Steffensen 加速法 (1) 求ln(x +sinx )=0的根。初值x0分别取0.1, 1,1.5, 2, 4进行计算。 (2) 求sin x =0的根。初值x0分别取1,1.4,1.6, 1.8,3进行计算。 分析其中遇到的现象与问题。 1.2计算过程和结果 1. 对方程ln(x+sinx)=0,可求解x+sinx=1的解。使用牛顿法,令f (x ) =x +sin x -1,则直至 |x k +1-x k |<1?10-5时,结束迭代;然后再使用基于牛顿法的steffensen f="" '(x="" )="1+cos" x="">1?10-5时,结束迭代;然后再使用基于牛顿法的steffensen> 加速法进行计算,直至|x k +1-x k |<1?10-5时,结束迭代。其迭代结果与迭代次数如下表所示(注n1为牛顿法迭代次数,n2为基于牛顿法steffensen>1?10-5时,结束迭代。其迭代结果与迭代次数如下表所示(注n1为牛顿法迭代次数,n2为基于牛顿法steffensen> 2. 对方程sinx=0,令f (x ) =sin x ,误差仍取|x k +1-x k |<> 1.3结果分析 1从牛顿与Steffensen 加速法可以看出,牛顿—Steffensen 加速法迭代速度明显快于牛顿 加速法,从计算结果可以看出,x=0.510973处取得精确值,而在x0=0.1,1处收敛速度较快, 说明牛顿法与基于牛顿法Steffensen 加速法在单根附近有较快的收敛速度,迭代法是否收敛,与初始近似值x0的好坏很有关系。 2从sinx=0使用两种方法可以看出,Steffensen 加速法迭代速度明显快于牛顿加速法,且牛顿法和基于牛顿法Steffensen 加速法在单根附近有较快的收敛速度,迭代法是否收敛,与初始近似值x0的好坏很有关系,如在x0=1.6处,得到的收敛解很大,这是因为 f '(1.6)=cos(1.6?180/π) ≈0,很难说明f(x)在此处是否发散或者是否收敛。 第二题 2.1题目 某过程测涉及两变量x 和y, 拟分别用插值多项式和多项式拟合给出其对应规律的近似多项式,已知xi 与yi 之间的对应数据如下,xi=1,2,…,10 yi = 34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795 103.5743 97.4847 78.2392 (1)请用次数分别为3,4,5,6的多项式拟合并给出最好近似结果f(x)。 (2)请用插值多项式给出最好近似结果 下列数据为另外的对照记录,它们可以作为近似函数的评价参考数据。 xi = Columns 1 through 7 1.5000 1.9000 2.3000 2.7000 3.1000 3.5000 3.9000 Columns 8 through 14 4.3000 4.7000 5.1000 5.5000 5.9000 6.3000 6.7000 Columns 15 through 17 7.1000 7.5000 7.9000 yi = Columns 1 through 7 42.1498 41.4620 35.1182 24.3852 11.2732 -1.7813 -12.3006 Columns 8 through 14 -18.1566 -17.9069 -11.0226 2.0284 19.8549 40.3626 61.0840 Columns 15 through 17 79.5688 93.7700 102.3677 2.2计算过程和结果 第一问 编好拟合程序,计算结果如下: (1)3次多项式拟合结果: 解得拟合多项式系数矩阵 e=[131.794399999964 -94.478681429651 19.333892249411 -1.032633896659] T * 故:s 3(x ) =131.7944-94.4787x +19.3338x 2-1.0326x 3 用MA TLAB 划出画出插值多项式函数曲线如下: f -s *22 =∑??f (x i ) -s (x i ) ??=3.2805e+003 * i =1 10 2 (2)4次多项式拟合结果: 解得拟合多项式系数矩阵 e=[ 0.74497500045210 73.53340190312746 -42.14325642464905 7.36797026998528 -0.38184564393842]T * (x ) =0.7450+73.5334x -42.1433x 2+7.3680x 3-0.3818x 4: 用因此四次拟合多项式为s 4 MATLAB 划出画出插值多项式函数曲线如下: 10 2 * f -s *22 =∑??f (x i ) -s (x i ) ??=878.5741 i =1 (3)5次多项式拟合结果: 解得拟合多项式系数矩阵 e =[-139.501853316416 304.728173432415 163.510704007400 034.501972019191 3.078900034678 0.098074705118] T 因此5次拟合多项式为: *s 5(x ) =-139.5019+304.7282x -163.5107x 2+34.5020x 3-3.0789x 4+0.0981x 5 用MA TLAB 划出画出插值多项式函数曲线如下: f -s *22 =∑??f (x i ) -s (x i ) ??= 128.3196 * i =1 10 2 (4)6次多项式拟合结果: e =[-18.69909465342829 66.37503804475497 -0.86696651297073 -16.89727867937389 5.11365302218162 -0.54078603814650 0.01935941646399] T 因此6次拟合多项式为: *s 6(x ) =-18.6991+66.3750x -0.8670x 2-16.8973x 3+5.1137x 4-0.5408x 5+0.0194x 6 用MA TLAB 划出画出插值多项式函数曲线如下: f -s *22 =∑??f (x i ) -s (x i ) ??=10.5683 * i =1 22 * 22 * f -s 4 22 * 22 10 2 * 比较不同次数的多项式拟合结果,可见f -s 6 说明6次多项式的最小二乘拟合的拟合程度最好. 因此最好的拟合多项式为: *s 6(x ) =-18.6991+66.3750x -0.8670x 2-16.8973x 3+5.1137x 4-0.5408x 5+0.0194x 6 第二问 用不同的插值方法进行插值多项式的求解。 (1)拉格朗日插值 求得9次拉格朗日插值多项式为: N 9(x ) =100?(-0.25204200011955+1.13500691817544x -0.8646529666577x 2+0.50113073368872x 3-0.2240802436669636x 4+0.05932495695820x 5-0.00879208298798x 6+0.00072285644857x 7-0.00003088405259x 8+0.00000053525132x 9) 用上面的插值多项式计算出被插点数据x (参考点)对应的y : Y= 42.3840 41.4947 35.0742 24.3601 11.2792 -1.7683 -12.2977 -18.1626 -17.9118 -11.0210 2.0333 19.8565 40.3584 61.0794 79.5709 93.7788 102.3713 计算所得到的值和参考值之间的方差为: D 17 (2)牛顿插值 使用9次牛顿插值多项式可得到同样的结果: N 9(x ) =100?(-0.25204200011955+1.13500691817544x -0.8646529666577x 2+0.50113073368872x 3-0.2240802436669636x 4+0.05932495695820x 5-0.00879208298798x 6+0.00072285644857x 7-0.00003088405259x 8+0.00000053525132x 9) N 9(x ) 即为最佳插值多项式近似结果。 2.3结果分析 * 1 比较不同次数的多项式拟合结果,可见f -s 6 22 * 22 * f -s 4 22 * 22 说明次数越高最小二乘拟合的拟合程度最好。 2牛顿插值多项式跟拉格朗日插值的结果完全一致, 这说明多项式是唯一的,因为插值节点有10个,九次插值需要得到10个常系数,因此方程有唯一解,所以插值多项式是唯一的。 第三题 3.1题目 5用Runge-Kutta 4阶算法对初值问题y/=-20*y,y(0)=1按不同步长求解,用于观察稳定区间的作用,推荐两种步长h=0.1,0.2。 注:此方程的精确解为:y=e-20x 3.2问题求解及过程 (1)给定初值x 为0,选取h=0.1进行计算,并对迭代过程和迭代图像进行分析。计算结 果如下: 请输入a 的值:0 请输入b 的值:1 请输入步长h :0.1 x[ 1]= 0.000000 y[ 1]= 1.000000 yi[ 1]= 1.000000 error[ 1]= 0.000000 x[ 2]= 0.100000 y[ 2]= 0.333333 yi[ 2]= 0.135335 error[ 2]= 0.197998 x[ 3]= 0.200000 y[ 3]= 0.111111 yi[ 3]= 0.018316 error[ 3]= 0.092795 x[ 4]= 0.300000 y[ 4]= 0.037037 yi[ 4]= 0.002479 error[ 4]= 0.034558 x[ 5]= 0.400000 y[ 5]= 0.012346 yi[ 5]= 0.000335 error[ 5]= 0.012010 x[ 6]= 0.500000 y[ 6]= 0.004115 yi[ 6]= 0.000045 error[ 6]= 0.004070 x[ 7]= 0.600000 y[ 7]= 0.001372 yi[ 7]= 0.000006 error[ 7]= 0.001366 x[ 8]= 0.700000 y[ 8]= 0.000457 yi[ 8]= 0.000001 error[ 8]= 0.000456 x[ 9]= 0.800000 y[ 9]= 0.000152 yi[ 9]= 0.000000 error[ 9]= 0.000152 x[10]= 0.900000 y[10]= 0.000051 yi[10]= 0.000000 error[10]= 0.000051 h=0.1 h=0.2 (2)给定初值x 为0,, 分别选取h=0.2进行计算,并对迭代过程和迭代图像进行分析。计算结果如下 请输入a 的值:0 请输入b 的值:2 请输入步长h :0.2 x[ 1]= 0.000000 y[ 1]= 1.000000 yi[ 1]= 1.000000 error[ 1]= 0.000000 x[ 2]= 0.200000 y[ 2]= 5.000000 yi[ 2]= 0.018316 error[ 2]= 4.981684 x[ 3]= 0.400000 y[ 3]= 25.000000 yi[ 3]= 0.000335 error[ 3]= 24.999665 x[ 4]= 0.600000 y[ 4]=125.000000 yi[ 4]= 0.000006 error[ 4]=124.999994 x[ 5]= 0.800000 y[ 5]=625.000000 yi[ 5]= 0.000000 error[ 5]=625.000000 x[ 6]= 1.000000 y[ 6]=3125.000000 yi[ 6]= 0.000000 error[ 6]=3125.000000 x[ 7]= 1.200000 y[ 7]=15625.000000 yi[ 7]= 0.000000 error[ 7]=15625.000000 x[ 8]= 1.400000 y[ 8]=78125.000000 yi[ 8]= 0.000000 error[ 8]=78125.000000 x[ 9]= 1.600000 y[ 9]=390625.000000 yi[ 9]= 0.000000 error[ 9]=390625.000000 x[10]= 1.800000 y[10]=1953125.000000 yi[10]= 0.000000 error[10]=1953125.000000 由图象和数据中可以看出,当步长为0.2时,迭代是不收敛的。 3.3结果分析 (1)用标准的四阶 Runge-Kutta 算法计算常微分方程时,不同的步长算法的稳定性有影响。不够稳定的步长下面的计算,误差会越来越大,结果失真严重。 (2)一般情况下,步长越小,标准的四阶 Runge-Kutta 算法的稳定性越高,精度也越高。 总结 通过这次上机练习,让我对数值分析一些方法及其理论有了更深层次的理解,通过对比,使我认识到各种方法的优缺点,并使我认识到学习数值分析中的不足,时的修补了自己知识上的漏洞,以便在复习与以后的工作中得以加强。同时,通过本次上机实习,复习了本科阶段所学习的C++编程,提高了编程的熟练度,并且学习了一门新的语言-matlab ,这些对我以后解决实际遇到的工程问题以及数学问题开辟了一个新的道路。 附 件 第一题 1.1.1第一问牛顿法 #include #include #define EPS 1.0e-5 int main() { double x=0.0; double y=0.0; int n=0; cout<><> cin>>x; do { y=x; x=x-(x+sin(x)-1)/(1+cos(x)); n++; } while(fabs(x-y)>=EPS); cout<><><><><> cout<><><> } 1.1.2第一问牛顿-Steffensen 法 #include #include #define EPS 1.0e-5 int main() { double x=0.0; double y=0.0; double z=0.0; double a=0.0; int n=0; cout<><> cin>>x; do { a=x; y=x-(x+sin(x)-1)/(1+cos(x)); z=y-(y+sin(y)-1)/(1+cos(y)); x=x-(y-x)*(y-x)/(z-2*y+x); n++; } while(fabs(x-a)>=EPS); cout<><><><><> cout<><><> } 1.2.1第二问牛顿法 #include #include #define EPS 1.0e-5 int main() { double x=0.0; double y=0.0; int n=0; cout<><> cin>>x; do { y=x; x=x-tan(x); n++; } while(fabs(x-y)>=EPS); cout<><><><><> cout<><><> } 1.2.2第二问牛顿-Steffensen 法 #include #include #define EPS 1.0e-6 int main() { double x=0.0; double y=0.0; double z=0.0; double a=0.0; int n=0; cout<><> cin>>x; do { a=x; y=x-tan(x); z=y-tan(y); x=x-(y-x)*(y-x)/(z-2*y+x); n++; } while(fabs(x-a)>=EPS); cout<><><><><> cout<><><> } 第二题 2.1.1最小二乘法求解 function [A norm]=min2mul(x,y,n,x_count,y_count) % x,y:原始数据点、A:拟合多项式各项常数项系数,n :拟合次数 x_count、y_count 拟合参考点 norm:残差的二范数 for i=1:n+1 for j=1:n+1 a(i,j)=0; for k=1:length(x) a(i,j)=a(i,j)+x(k)^(i-1)*x(k)^(j-1); end end end for i=1:n+1 b(i)=0; for k=1:length(y) b(i)=b(i)+x(k)^(i-1)*y(k); end end A=a\b'; for i=1:length(xi) y_count_cal(i)=0; for j=1:n+1 y_count_cal(i)=y_count_cal(i)+A(j)*x_count(i)^(j-1); end end for i=1:length(x) y_cal(i)=0; for j=1:n+1 y_cal(i)=y_cal(i)+A(j)*x(i)^(j-1); end end Y=[y y_count]; Y_cal=[y_cal y_count_cal]; norm=sum((Y-Y_cal).^2) 2.2.1拉格朗日差值多项式拟合 C++ #include #include #include using namespace std; int main(){ int n,i,j; cout<><> cin>>n; float X[n+1],Y[n+1],x,y,h[n],l; cout<><> for(i=0;i<> cin>>X[i]>>Y[i]; } for(i=0;i<=n;i++){ 求拉格朗日l="" 的分子项乘以相应的y="" 值(方便写出插值多项式)="" h[i]="">=n;i++){> for(j=0;j<> if(j!=i) h[i]*=(X[i]-X[j]); } cout<><<'><'> } while(cin>>x){ //用插值近似函数计算 y=0; for(i=0;i<> l=1; for(j=0;j<> if(j!=i) l=l*(x-X[j])/(X[i]-X[j]); } y=y+l*Y[i]; } cout<><><> } return 0; } 2.2.2牛顿插值 #include #include using namespace std; #include int main(){ int n,i,j; float h,x,y,x0,t,temp1,temp2,T; cout<><> cin>>n; float X[n],Y[n]; cout<><> for(i=0;i<> cin>>X[i]>>Y[i]; } h=X[1]-X[0]; //步长 while(cin>>x){ cout<><> cin>>x0; t=(x-x0)/h; y=0; for(i=0;i<> temp1=Y[i]; for(j=i+1;j<> temp2=Y[j]; Y[j]=temp2-temp1; temp1=temp2; } } T=1; for(i=1;i cout<><><><<'><'> cout<> for(i=0;i<> y+=T*Y[i]; T*=t; t--; } cout<><><> } return 0; } 第三题 3.1.1Runge-Kutta 4阶算法 a=input(' 请输入a 的值:'); b=input(' 请输入b 的值:'); h=input(' 请输入步长h:'); x=[a:h:b]; y(1)=1; yi(1)=exp(-20*x(1)); error(1)=0; n=(b-a)/h+1; for i=2:n k1=-20*y(i-1)*h; k2=-20*(y(i-1)+k1/2)*h; k3=-20*(y(i-1)+k2/2)*h; k4=-20*(y(i-1)+k3)*h; y(i)=y(i-1)+(k1+2*k2+2*k3+k4)/6; yi(i)=exp(-20*x(i)); error(i)=y(i)-yi(i); fprintf('x[%2.0f]=%10.6f y[%2.0f]=%10.6f yi[%2.0f]=%10.6f error[%2.0f]=%10.6f\n',i-1,x(i-1),i-1,y(i-1),i-1,yi(i-1),i-1,error(i-1)) end plot(x,y,'g' ,x,yi, 'r' ) 范文二:数值分析上机报告
范文三:数值分析上机报告
范文四:数值分析上机报告
范文五:数值分析上机报告