当前位置:首页 >> 工学 >>

数值分析上机实习报告


(数值分析上机实验报告)

院 专 班 姓 学

系: 业: 级: 名: 号:

矿业学院 矿业工程 2015 王 2015022 代

指导教师:

第一题
1.用
Newton 法求解方程 x7 ? 28x 4 ? 14 ? 0 ,在(0.1,1.9)的近似根(初始近似值取为区间端点, 迭代 6 次或误差小于 0.00001)。 1.1 理论依据及方法应用条件 Newton 迭代法:由一般迭代函数 ?s ( x) ? x ? ? (?1) j
j ?1 s ?1

f ( x) f ( j ) [ f ( x)] 取 s=2 时, 有 ? 2 ( x) ? x ? ' , [ f ( x)] j , j! f ( x)

可得二阶迭代序列 xk ?1 ? xk ?

f (xk ) ,此种迭代法称为 Newton 迭代法。 f ' ( xk )

定理:设函数在有限区间[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 ) ? xk ? f ( xk ) f ' ( xk )

k=0,1,2…

所产生的迭代序列{x0}平方收敛于方程 f ( x) =0 区间[a,b]上的唯一解α 。 推论:设函数 f(x)满足定理中条件Ⅰ,Ⅱ,Ⅲ,若选初值 x0 ,使 f ( x0 ) · f ' ' ( x0 ) >0,则 Newton 迭代过程
xk ?1 ? ? ( xk ) ? xk ? f ( xk ) f ' ( xk )

(k=0,1,2…) 产生的迭代序列{xk}单调收敛于 f ( x) =0 的唯一解α 。

1.2 计算程序
# include <iostream.h> # include <iomanip.h> # include <string> # include <math.h> using namespace std; double *newton (double a,double b,double eps); double newtonz (double x); void main () { double a=0.1,b=1.9,eps=0.00001,*result;

//牛顿迭代函数 //牛顿迭代子函数

//初始数据

cout<<"\n 牛顿法解方程:x^7-28x^4+14=0,在(0.1,1.9)中求近似根,初始值为区间端点,\n 误差为 0.00001。 \n"<<endl; cout<<"学号:2014021966 姓名:徐林\n"<<endl;

result = newton (a,b,eps); if (a<=result[0]&&result[0]<=b) cout<<"近似根为:"<<result[0]<<endl; if (a<=result[1]&&result[1]<=b) cout<<"近似根为:"<<result[1]<<endl; //------------------------------------------cout<<"\n"<<"结束,按任意键关闭"<<endl; getchar(); }//主函数结束 //******************************************************************* double newtonz (double x) //牛顿迭代子函数 { double x1=0.0,t; t = (7*pow(x,6)-4*28*pow(x,3)); if (t==0) exit (0); x1 = (x-((pow(x,7)-28*pow(x,4)+14)/t)); return x1; } double *newton (double a,double b,double eps) { double x0=0.0,x1=1.0,x2=0.0,re[2]; int k=0; x0 = a; while (x0>eps) { k++; x2 = x1; x1 = newtonz(x1); x0 = fabs(x1-x2); }re[0] = x1; x0 = b, k = 0; while (x0>eps) { k++; x2 = x1; x1 = newtonz(x1); x0 = fabs(x1-x2); }re[1] = x1; return re; } //牛顿迭代函数

//代入 a 迭代计算

//调用牛顿迭代子函数

//代入 b 迭代计算

//调用牛顿迭代子函数

1.3 计算结果打印

1.4 MATLAB 上机程序
function y=Newton(f,df,x0,eps,M) d=0; for k=1:M if feval(df,x0)==0 d=2;break else x1=x0-feval(f,x0)/feval(df,x0); end e=abs(x1-x0); x0=x1; if e<=eps&&abs(feval(f,x1))<=eps d=1;break end end if d==1 y=x1; elseif d==0 y='迭代 M 次失败'; else y= '奇异' end function y=df(x) y=7*x^6-28*4*x^3; End function y=f(x) y=x^7-28*x^4+14; End

>> x0=1.9; >> eps=0.00001; >> M=100; >> x=Newton('f','df',x0,eps,M); >> vpa(x,7)

1.5 问题讨论 1.需注意的是,要使用 Newton 迭代法须 f ( x) ? x7 ? 28x4 ? 14 满足定理中的条件Ⅰ ,Ⅱ,Ⅲ,以及

f ( x0 ) · f ' ' ( x0 ) >0。要用误差范围来控制循环的次数,保证循环的次数和质量,编写程序过程中要
注意标点符号的使用,正确运用适当的标点符号,Newton 迭代法是局部收敛的,在使用时应先确定 初始值,否则所得的解可能不在所要求的范围内。

(3)因为 newton 法求方程是平方收敛的,所以较为精确,但是要求出函数的导数,且必须有二阶导数。

第二题

2.已知函数值如下表:
x
f ( x) x
f ( x) f ' ( x)

1 0 6 1.7917595 f ' (1) =1

2 0.69314718 7 1.9459101

3 1.0986123 8 2.079445

4 1.3862944 9 2.1972246

5 1.6094378 10 2.3025851 f ' (10) =0.1

试用三次样条插值求 f (4.563) 及 f ' (4.563) 的近似值。 2.1 理论依据及方法应用条件 三次样条插值函数可定义为:对于[a,b]上的一个划分∏ a<x0<x1<x2<…<xn-1<xn=b.(n>=2) 如果定义在[a,b]上的函数 S(x),满足(1).在[xi,xi+1]上为 3 次多项式;(2). S(x), S'(x), S"(x)在[a,b]上连续,则称 S(x)在 [a,b]上划分 ? 的 3 次样条函数,如果对于 f ( x) ? [a, b] ,
s( x) ? f ( x) 还满足 s( xi ) ? f ( x) , i ? 0 ~ n ,则称 s( x) 为 f ( x) 的三次样条插值函数。

其基本思想是对均匀分划的插值函数的构造,三次样条函数空间中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由解的唯一性就能求得S(x)。 由 s(xi)=yi,I=1,2,…N S(xi)= s’(x0)=y0’,s’(xN)=yN’可得 cjΩ 3((xi-xj)/h)=yi
N ?1

? cjΩ 3((x-xj)/h)这样对不同插值
N ?1 j ? ?1

?

N ?1

j ? ?1

S’(x0)=1/h ? cjΩ 3’((x0-xj)/h)=y’0
j ? ?1

S’(xN)=1/h ? cjΩ 3’((xN-xj)/h)=y’N
j ? ?1

N ?1

?2 ?? ? 1 ?? ? ? ? ?

1 2 ? 1 ?

?

?

N

? 2 1

? ? ? ? ? ? M 0 ? ? d0 ? ?? ? ? ? ?? M ? ? d ? 1 1 ? ?? ? ? ? ? ? ? ?? ? ? ? ? N ? 1? ? ? ? ? ? ? 2 ? ? ? M ? ?d ? ? N? ? ? ? ? N? ?

? 6 y1 ? y0 ? d ? ? y' ) ? 0 h ( h 0 0 0 ? ? y ?y y ?y 6 j ?1 j j j ?1 ? d ? ( ? ) ? 6 f (x ,x ,x ) ? 0 h j ?1 j j ?1 ?h h h ? j ?1 j j j ?1 ? (y ?y ) 6 ? N N ?1 ] [ y' ? ?d N ? N h h ? N ?1 N ?1 ? h ? j ?1 ?? j ? h ?h ? j ?1 j ? h ? j ?? ? 1 ? ? ? j h ?h ? j j ?1 j ?

2.2 计算程序 #include<stdio.h> #include<math.h> #define N 10 /*宏定义*/ main() { float s,ds,t; float dy0=1,dy9=0.1; int j; int x[N]={1,2,3,4,5,6,7,8,9,10}; float y[N]={0,0.69314718,1.0986123,1.3862944,1.6094378, 1.7917595,1.9459101,2.079445,2.1972246,2.3025851}; int b[N]={2,2,2,2,2,2,2,2,2,2},h[N-1]; float d[N],u[N-1],v[N-1],a[N-1],c[N-1],B[N],l[N-1],p[N],X[N]; for(j=1;j<=9;j++) h[j-1]=x[j]-x[j-1]; d[0]=6/h[0]*(y[1]/h[0]-y[0]/h[0]-dy0); d[9]=6/h[8]*(dy9-y[9]/h[8]+y[8]/h[8]); for(j=1;j<=8;j++) d[j]=6/(h[j-1]+h[j])*(y[j+1]/h[j]-y[j]/h[j]-y[j]/h[j-1]+y[j-1]/h[j-1]); for(u[8]=1,j=0;j<=7;j++) u[j]=h[j-1]/(h[j-1]+h[j]); for(v[0]=1,j=1;j<=8;j++) v[j]=h[j]/(h[j-1]+h[j]); for(j=0;j<=8;j++) a[j]=u[j]; for(j=0;j<=8;j++) c[j]=v[j]; for(B[0]=b[0],j=1;j<=9;j++)
/*追赶法求解三弯矩方程*/

B[j]=b[j]-a[j]/B[j-1]*c[j-1]; for(j=1;j<=9;j++) l[j]=a[j]/B[j-1]; for(j=1;j<=9;j++) p[j]=d[j]-l[j]*p[j-1]; X[9]=p[9]/B[9]; for(j=8;j>=0;j--) X[j]=p[j]/B[j]-c[j]*X[j+1]/B[j]; t=4.563; s=X[3]*pow((x[4]-t),3)/6/h[3]+X[4]*pow((t-x[3]),3)/6/h[3]+ (y[3]-X[3]*h[3]*h[3]/6)*(x[4]/h[3]-t/h[3])+ (y[4]-X[4]*h[3]*h[3]/6)*(t/h[3]-x[3]/h[3]);
*/

/*解 f(x)的值

ds=-X[3]*pow((x[4]-t),2)/2/h[3]+X[4]*pow((t-x[3]),2)/2/h[3](y[3]-X[3]*h[3]*h[3]/6)/h[3]+(y[4]-X[4]*h[3]*h[3]/6)/h[3];
的值*/

/* 解 f’ (x) /*打印结果*/

printf("s=%f\nds=%f\n ",s,ds); }

2.3 计算结果打印

2.4 MATLAB 上机程序
function Q=san(ssss,p) Q=zeros(2,1); x=[1;2;3;4;5;6;7;8;9;10]; y=[0;0.69314718;1.0986123;1.3862944;1.6094378;1.7917595;1.9459101;2.079445;2.1972246;2.3025851]; h=zeros(10,1); d=zeros(10,1); u=zeros(10,1); v=zeros(10,1); r=zeros(10,1); l=zeros(10,1); z=zeros(10,1); m=zeros(10,1); for t=1:1:9; h(t)=x(t+1)-x(t); end d(1)=6/h(1)*((y(2)-y(1))/h(1)-1); d(10)=6/h(9)*(0.1-(y(10)-y(9))/h(9));

for t=1:1:8 u(t+1)=h(t)/(h(t)+h(t+1)); v(t+1)=1-u(t+1); d(t+1)=6/(h(t)+h(t+1))*((y(t+2)-y(t+1))/(x(t+2)-x(t+1))-(y(t+1)-y(t))/(x(t+1)-x(t))); end u(10)=1;v(1)=1;r(1)=d(1); for t=2:1:10 l(t)=u(t)/r(t-1); r(t)=d(t)-l(t)*v(t-1); end z(1)=d(1); for t=2:1:10 z(t)=d(t)-l(t)*z(t-1); end m(10)=z(10)/r(10); for t=9:-1:1 m(t)=(z(t)-v(t)*m(t+1))/r(t); end for t=1:1:10 if p>=t&&p<(t+1) Q(:,1)=feval(ssss,p,t,x,m,h,y);break end end function Q=ssss(p,t,x,m,h,y) Q=zeros(2,1); Q(1,1)=((power((x(t+1)-p),3)*m(t)+power((p-x(t)),3)*m(t+1))/6+(y(t)-m(t)*h(t)*h(t)/6)*(x(t+1)-p)+(y(t+1)-m(t+1)*h(t)*h(t)/6)*( p-x(t)))/h(t); Q(2,1)=(-(power((x(t+1)-p),2)*m(t)+power((p-x(t)),2)*m(t+1))/2+(y(t)-m(t)*h(t)*h(t)/6)+(y(t+1)-m(t+1)*h(t)*h(t)/6))/h(t); end

2.5 问题讨论 1.若要用追赶法求解三对角方程组,三对角阵需要满足: Ai (i=1,2,…,n)均非奇异,保证 A 有唯一的 Doolittle 分解; ai ci ≠0; 2.样条插值效果比 Lagrange 插值好, 三次样条插值的解存在且唯一,近似误差较小.并且没有 Runge 现象。

第三题
3.用 Romberg 算法求 ?1 3 x x1.4 (5x ? 7) sin x 2 dx (允许误差ε
3

=0.00001)。

? (0) b ? a 3.1 理论依据及方法应用条件 ? f (a) ? f (b)? ?T1 ? 2 数值积分的 Romberg 算法计算步骤如下 : ? l ?1 ? ? 1 b ?a 2 b ? a ?? ? (l ) ? ( l ?1) T ? T ? f a ? (2i ? 1) l ? ? ?1 ?1 l ?1 ? ? 2? 2 i ?1 ? 2 ?? ? ( k ? 1 ) ? m (k ) ?T ( k ?1) ? 4 Tm ? Tm , m ? 1,2,? ? ?, l k ? 1,2,? ? ?, l ? m ? 1 m ?1 ? 4m ? 1 ?
当 Tl
(0)

? Tl ?1

(0)

? ? 时,就停机。

3.2 计算程序 #include<stdio.h> #include<math.h> #define N 9 float f(float x) { float y; y=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x); return(y); } main() { float T[N+1][N+1],h[N+1],a=1,b=3,m[N+1]; int i,l; T[1][0]=(b-a)*(f(a)+f(b))/2; l=1; while(l<=N) { m[l]=0; for(i=1;i<=(pow(2,l-1));i++) m[l]+=f(a+(2*i-1)*(b-a)/pow(2,l)); T[1][l]=(T[1][l-1]+(b-a)*m[l]/pow(2,l-1))/2; l++; } i=1; while(i<=N) { for(l=1;l<=N-i+1;l++) T[i+1][l-1]=(pow(4,i)*T[i][l]-T[i][l-1])/(pow(4,i)-1); h[i]=T[i][0]-T[i+1][0]; if(fabs(h[i])<=1e-5) break; i++;
/*定义函数 f(x)*/

数值分析上机实习报告

2012 年 12 月

} printf("The answer is:%f",T[i+1][0]); } 3.3 计算结果打印

3.4 MATLAB 上机程序
function [T,n]=mromb(f,a,b,eps) if nargin<4,eps=1e-6;end h=b-a; R(1,1)=(h/2)*(feval(f,a)+feval(f,b)); n=1;J=0;err=1; while (err>eps) J=J+1;h=h/2;S=0; for i=1:n x=a+h*(2*i-1); S=S+feval(f,x); end R(J+1,1)=R(J,1)/2+h*S; for k=1:J R(J+1,k+1)=(4^k*R(J+1,k)-R(J,k))/(4^k-1); end err=abs(R(J+1,J+1)-R(J+1,J)); n=2*n; end R; T=R(J+1,J+1);

format long f=@(x)(3.^x)*(x.^1.4)*(5*x+7)*sin(x*x); [T,n]=mromb(f,1,3,1.e-5)

3.5 问题讨论 1、Romberge 算法的优点是: a、把积分化为代数运算,而实际上只需求 T1(i),以后用递推可得。 b、算法简单且收敛速度快,一般 4 或 5 次即能达到要求。 c、节省存储量,算出的可存入。 2、Romberge 算法的缺点是: a、对函数的光滑性要求较高。 b、计算新分点的值时,这些数值的个数成倍增加。

第 11 页

数值分析上机实习报告

2012 年 12 月

第四题
4.用定步长四阶 Runge-Kutta 法求解

dy1 / dt ? 1

dy2 / dt ? y3 dy3 / dt ? 1000? 1000y2 ? 100y3
y1 (0) ? 0 y2 (0) ? 0

y3 (0) ? 0
h ? 0.0005 ,打印 yi (0.025) , yi (0.045) , yi (0.085) , yi (0.1) , (i ? 1,2,3)

4.1 理论依据及方法应用条件 Runge-Kutta 法的基本思想: y j ?1 不是按 Taylor 公式展开, 而是先写成 t n 处附近的值的线性组合 (有 待定系数),再按 Taylor 公式展开,然后确定待定常数。 四阶古典 Runge-Kutta 公式:
? 1 ? ?Yn ?1 ? Y n ? 6 ( K1 ? 2 K 2 ? 2 K 3 ? K 4 ) ? ? K1 ? hF ( xn , Yn ) ? 1 1 ? ? K 2 ? hF ( xn ? h, Yn ? K1 ) 2 2 ? 1 1 ? ? K 3 ? hF ( xn ? 2 h, Yn ? 2 K 2 ) ? ? K ? hF ( x ? 1 h, Y ? 1 K ) 4 n n 3 ? 2 2 ?

4.2 计算程序 #include <stdio.h> int main() { int i; double h=0.0005; double k1,k2,k3,k4; double y1=0.0,y2=0.0,y3=0.0; for(i=1;i<=200;i++) { k1=k2=k3=k4=h*1.0; y1+=(k1+2*k2+2*k3+k4)/6; k1=k2=k3=k4=h*y3; y2+=(k1+2*k2+2*k3+k4)/6;
第 12 页

数值分析上机实习报告

2012 年 12 月

k1=h*(1000-1000*y2-100*y3); k2=h*(1000-1000*y2-100*(y3+0.5*k1)); k3=h*(1000-1000*y2-100*(y3+0.5*k2)); k4=h*(1000-1000*y2-100*(y3+k3)); y3+=(k1+2*k2+2*k3+k4)/6; if(i==50) { printf("\ny1(0.025)=%f continue; } if(i==90) { printf("\ny1(0.045)=%f continue; } if(i==170) { printf("\ny1(0.085)=%f continue; } if(i==200) printf("\ny1(0.100)=%f } }

y2(0.025)=%f

y3(0.025)=%f",y1,y2,y3);

y2(0.045)=%f

y3(0.045)=%f",y1,y2,y3);

y2(0.085)=%f

y3(0.085)=%f",y1,y2,y3);

y2(0.100)=%f

y3(0.100)=%f\n\n",y1,y2,y3);

4.3 计算结果打印

4.4 MATLAB 上机程序
function Y=R_K(df1,a,b,h) m=(b-a)/h; Y=zeros(3,1); S=zeros(3,1); K=zeros(3,4);
第 13 页

数值分析上机实习报告

2012 年 12 月

x=a;y1=a;y2=a;y3=a; for n=1:m K(:,1)=feval(df1,x,y1,y2,y3); x=x+0.5*h;S(:,1)=Y+0.5*h.*K(:,1); y1=S(1,1);y2=S(2,1);y3=S(3,1); K(:,2)=feval(df1,x,y1,y2,y3); S(:,1)=Y+0.5*h.*K(:,2); y1=S(1,1);y2=S(2,1);y3=S(3,1); K(:,3)=feval(df1,x,y1,y2,y3); x=x+0.5*h;S(:,1)=Y+h.*K(:,3); y1=S(1,1);y2=S(2,1);y3=S(3,1); K(:,4)=feval(df1,x,y1,y2,y3); Y=Y+h.*(K(:,1)+2.*K(:,2)+2.*K(:,3)+K(:,4))/6; end function Z=df1(x,y1,y2,y3) Z=zeros(3,1); Z(1)=0*x+0*y1+0*y2+0*y3+1; Z(2)=0*x+0*y1+0*y2+1*y3; Z(3)=0*x+0*y1-1000*y2-100*y3+1000; end

4.5 问题讨论 1.定步长四阶 runge-kutta 法稳定,精度高,可根据有 y ' ? f (t , y ) 变化的情况与需要的精度自动修 改步长,误差小且程序简单,存储量少。 2.但是 Runge-Kutta 法需要每步都计算函数值 f (t , y ) 四次,在函数较复杂时,工作量就会变得较大, 可靠性有待核查。

第 14 页

数值分析上机实习报告

2012 年 12 月

第五题
5.已知 A 与 b
12.38412 2.115237 -1.061074 1.112336 2.115237 19.141823 -3.125432 -1.012345 2.189736 1.563849 -0.784156 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 -86.612343 0.718719 1.563849 1.836742 -3.741856 0.431637 9.789365 -0.103458 -1.103456 0.238417 1.1101230 1.742382 -0.784165 -1.056781 2.101023 -3.111223 -0.103458 14.713846 3.123789 -2.213474 4.719345 3.067813 1.112348 0.336993 -0.71828 2.121314 -1.103456 3.123789 30.719334 4.446782
T

-2.031743 3.123124 -1.010103 -0.037585 1.784317 0.238417 -2.213474 4.446782 40.00001

A=

-0.113584 0.718719 1.742382 3.067813 -2.031743

b=(2.1874369

33.992318

-25.173417 0.84671695 1.784317

-5.6784392)

(1).用列主元素消去法求解 Bx=b. 5.1
(1)理论依据 用 Househloder 变换,把 A 化为三对角阵。 设 A=(aij), aij=aji , A=(a1 , a2 , … , an),ai=(a1i , a2i , … , ani) ,取第一列 a1 中 a11 不动,把(a21 , a31 , … ,an1) 化为(±s1 , 0 , 0 , … , 0) ,这里 s1
T T T

?

?a
i ?2

n

2 i1

,前面的正负表示镜像可取为两个方向相反的向量,这增加计算的灵活性。

记 d1=(a11 ,±s1 , 0, … ,0) ,则 d1=H1a1 , 这里 H1
T

? E ? 2?1?1T , ?1 ?

a1 ? d1 u ? 1 , a1 ? d1 u1

T u1 ? (a1 ? d1 )T ? (0, a21 ? s1 , a31 ,..., an1 ),

u1

2

2 2 ? a1 ? d1 ? (a21 ? s1 ) 2 ? a31 ? ... ? an 1,

2

2 2 2 2 ? a21 ? s12 ? 2a21s1 ? a31 ? ... ? an 1 ? 2s1 ? 2a21s1 .

为了增加计算的灵活性,避免同号数相减,选取符号使 2a21s1>0,这样只需要把上式右端第二项取为 2sign(a21)a21 s1 即可,从而
T u1 ? (a1 ? d1 )T ? (0, a21 ? sign(a21 )s1, a31,..., an1 ).

归纳起来算法步骤为:
(1) (r) (1)令A 0 ? A,aij ? aij ,已知A r-1,即A r-1 ? (a ij ).

? n ? (r ) (2) sr ? ? ? (air )2) ? . ? i ? r ?1 ?
r) (3)ar ? sr2 ? ai(? 1,r s r . r) (r ) (r ) (r ) T ur ? [0,..., 0, ai(? 1,r ? sin g ( ai ?1,r ) sr , ai ? 2,r ,..., anr ] .

1/ 2

(4) yr ? A r ?1ur / ar . 1 (5)kr ? urT yr / ar . 2 (6)qr ? yr ? kr ur . (7)A r ? A r-1 ? (qr urT +ur qrT ) , r=1 , 2 , ... , n-2.
第 15 页

数值分析上机实习报告

2012 年 12 月

(2)程序主体 #include <stdio.h> #include <math.h> void main( ) { int i,j,m,r,sign; double A0[9][9],s,z,p,n,v,h,l,y[9],u[9],k,q[9],X[9],x[9]={0,0,0,0,0,0,0,0,0},B[9][9],g[9]; double A[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.12312 4}, {-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.0101 03}, {1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.03758 5}, {-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.741865,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.213 474}, {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.0000 1} }; double b[9]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345, -5.6784392}; for(i=0;i<9;i++) { for(j=0;j<9;j++) { A0[i][j]=A[i][j]; } } for(r=0;r<7;r++) { s=0; for(i=r+1;i<9;i++) { s=s+A[i][r]*A[i][r];
第 16 页

数值分析上机实习报告

2012 年 12 月

} s=sqrt(s); l=s*s+fabs(A[r+1][r])*s; if(A[r+1][r]>0) sign=1; else if(A[r+1][r]<0) sign=-1; for(i=0;i<9;i++) { if(i<=r) u[i]=0; else if(i==r+1) u[i]=A[r+1][r]+sign*s; else u[i]=A[i][r]; } for(i=0;i<9;i++) { y[i]=0; for(j=0;j<9;j++) { y[i]=y[i]+A[i][j]*u[j]; } y[i]=y[i]/l; } k=0; for(i=0;i<9;i++) { k=k+u[i]*y[i]; } k=k/(2*l); for(i=0;i<9;i++) { q[i]=y[i]-k*u[i]; } for(i=0;i<9;i++) { for(j=0;j<9;j++) { A[i][j]=A[i][j]-q[i]*u[j]-q[j]*u[i]; } } } printf("\nHousehold 变换矩阵 B 为:\n"); for(i=0;i<9;i++) { printf("\n"); for(j=0;j<9;j++) { printf("%f,",A[i][j]); } } }
第 17 页

数值分析上机实习报告

2012 年 12 月

(3)结果截屏

(4) MATLAB 上机程序
unction [x]=mgauss2(A,b,flag) if nargin<3,flag=0;end n=length(b); for k=1:(n-1) [ap,p]=max(abs(A(k:n,k))); p=p+k-1; if p>k A([k p],:)=A([p k],:); b([k p],:)=b([p k],:); end m=A(k+1:n,k)/A(k,k); A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n); b(k+1:n)=b(k+1:n)-m*b(k); A(k+1:n,k)=zeros(n-k,1); if flag~=0,Ab=[A,b],end end x=zeros(n,1); x(n)=b(n)/A(n,n); for k=n-1:-1:1 x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k); end

format long A=[12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743;
第 18 页

数值分析上机实习报告

2012 年 12 月

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]; b=[2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392]'; x=mgauss2(A,b);x=x'

(5)问题讨论 1、矩阵 A 是严格对角占优阵,而且还是实数对称阵,满足 Householder 变换的条件,并且得到 的 B 就是与 A 相似的三对角阵。 2、经过 Householder 变换矩阵有良好的性质,便于计算分析。

第 19 页


赞助商链接
相关文章:
数值分析上机实习报告
数值分析上机实习报告 学姓专 号: 名: 业: 联系电话: 联系电话: 2011 年 11 月 (数值分析)上机实习报告 第 I 页 序 言 1. 程序语言:我的数值分析上机实...
西南交通大学研究 数值分析上机实习报告2012
数值分析上机实习报告数值分析上机实习报告隐藏>> 数值分析上机实习报告要求 1.应提交一份完整的实习报告。具体要求如下: (1)要有封面,封面上要标明姓名、学号、专...
西南交通大学数值分析上机实习报告要求
西南交通大学数值分析上机实习报告要求 - 数值分析 2016 上机实习报告要求 1.应提交一份完整的实习报告。具体要求如下: (1) 报告要排版,美观漂亮;要标明姓名、...
数值分析上机实习报告
数值分析上机实习报告 姓名: 学号: 专业: 联系电话: 目 录 序言... 3 (一) 雅可比迭代法和高斯-塞得尔迭代法的收敛性和收敛速度...
数值分析上机实习报告
disp('最小二乘多项式拟合');disp(i+2); disp('次多项式的方差为:');disp(S(i)); end end 数值分析上机实习报告要求 13 1.应提交一份完整的实习报告。...
2014级硕士研究生数值分析上机实习报告(答案)
2014级硕士研究生数值分析上机实习报告(答案)_理学_高等教育_教育专区。哈尔滨工业大学(威海)实验报告纸 2014 级硕士研究生数值分析上机实习 (第一次) 姓名: 学号...
数值分析上机实习报告
数值分析上机实习报告_理学_高等教育_教育专区。数值分析上机实习报告 上机习题解答 指导教师: 指导教师: 姓学专名: 号: 业: 联系电话: 联系电话: 上海交通大学...
数值分析上机实习报告
n 阶方程组的系数矩阵为: 二 计算程序: 《 数值分析上机实习题 1 househloder 变换将 A 化为三对角阵 B #include <stdio.h> #include <math.h> main...
2014级硕士研究生数值分析上机实习报告
2014级硕士研究生数值分析上机实习报告_理学_高等教育_教育专区。哈尔滨工业大学(威海)实验报告纸 2014 级硕士研究生数值分析上机实习 (第一次) 姓名: 学号: 学院...
数值分析上机实习报告
数值分析上机实习题 if(fabs(k)>=1.5) h=0; else if(fabs(k)<=0.5) h=-k*k+3.0/4; else h=(k*k)/2-(3*fabs(k))/2+9.0/8; ...
更多相关标签: