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

《数值分析》上机实习报告


数值分析程序设计作业

一.housholder 变换.求方程的解 二.求特征值 一,程序要求: (1)用 Household 变换,把 A 化为三对角阵 B(并打印 B). (2)用超松弛法求解 Bx=b(取松弛因子ω=1.4,X(0)=0,迭代 9 次) (3)用列主元消去法求解 Bx=b。 (4) 用对分法求 B(即 A)的小于 23 且最接近 23 的特征值的近似值(误差小于 0.1) (5) 用反幂法求 B 的该特征值的更精确近似值及相应的特征向量. 二,解题算法 1:Household 法 (1) 令 A0=A, aij(1)=aij,已知 Ar-1 即 Ar-1=(aij(r)) (2) Sr=( (air(r ))2)1/2 (3) αr=Sr2+|a(r)r+1,r|Sr ur=[0,0,a(r)r+1,r+Sign(a(r)r+1,r)Sr,a(r)r+2,r,…,anr(r)]T (4) yr=Ar-1ur/αr (5) kr= urTyr/2αr (6) qr=yr-krur (7) Ar=Ar-1-(qrurT+urqrT) r=(1,2,……n-2) 2:超松弛法 其基本思想是迭代,即在高斯方法已求出 x(m),x(m-1)的基础上,组合新的序列,从而加快 收敛速度。其算式是: xi(m)=(1-ω)xi(m-1)+ω( bijxi(m)+ xj(m-1)+gi) 其中ω是超松弛因子,当ω>1 时,可以加快收敛速度。 3:高斯消去法: 对矩阵作恰当的调整,选取绝对值尽量大的元素作为主元素。然后把矩阵化 为上三角阵,再进行回代,求出方程的解。 对分法和反幂法求解特征值及特征向量 4. 对分法:对于对称的三对角方阵,其特征值必落在[-‖B‖, ‖B‖]中,其中‖B‖是其任意 范数.通过同号数的计算,可知在[0, ‖B‖]中特征根的个数,在不断对分,就可得特征根所在 的区间,如在[lk,uk]中有一根λk,且[lk,uk]很小,就可取 lk+uk/2 作为λk 的近似值. 5. 反幂法:基于乘幂法中求 A-1 不太容易,因此不直接用 A-1xk=xk+1,而使用 Axk+1=xk.用求 解线性方程组求出 xk+1,再标准化,就可得到绝对值最小的特征值与相应的特征向量,这就是 反代法(也称反幂法).

三,程序清单 #include "math.h" #include "io.h"

#include "stdio.h" #include "string.h" #define N 9 float sign(float x); void guss(float a[N][N]); float fanmi(float a[N][N]); float sor(float a[N][N],float b[N]); float duifen(); void housholder(); int f(float x); float e[N],o[N]; float approx; float b[N]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4. 719345,-5.6784392}; float a[N][N]= { 12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031 743, 2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.12 3124, -1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.0 10103, 1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.03 7585, -0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784 317, 0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238 417, 1.742382,-0.784156,-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.4467 82, -2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.0 0001

}; void housholder() { float Ur[N],Yr[N],qr[N]; int i,j,l,k; float sr1=0,sr=0.0,ar=0.0,*p,kr=0.0; for(i=0;i<N-2;i++) { for(j=i+1;j<N;j++) sr1+=a[j][i]*a[j][i]; sr=sqrt(sr1); ar=sr1+fabs(a[i+1][i])*sr; for(l=0;l<N;l++) { if (l<=i) Ur[l]=0.0; else if(l==i+1) Ur[l]=a[l][i]+sign(a[l][i])*sr; else Ur[l]=a[l][i]; } for(k=0;k<N;k++) { Yr[k]=0; } for(k=0;k<N;k++) for(l=0;l<N;l++) Yr[k]+=a[k][l]/ar*Ur[l]; for(k=0;k<N;k++) kr+=Ur[k]*0.5/ar*Yr[k]; for(k=0;k<N;k++) qr[k]=Yr[k]-kr*Ur[k]; for(k=0;k<N;k++) for(l=0;l<N;l++) a[k][l]-=(qr[k]*Ur[l]+Ur[k]*qr[l]); sr1=0;kr=0;ar=0.0; } system("cls"); for(k=0;k<N;k++) { for(l=0;l<N;l++) printf("%.4f ,",a[k][l]); printf("\n"); } for(i=0;i<N;i++) o[i]=a[i][i]; for(i=0;i<N-1;i++) e[i]=a[i][i+1]; } float sor(float a[N][N],float b[N]) /*超松驰法求解的子函数*/ { float a1[N][N]; int i,j,k,m; float x[N],temp[N][N]; float w=1.4,h=0,g=0;

for(m=0;m<N;m++) x[m]=0; for(i=0;i<N;i++) for(j=0;j<N;j++) a1[i][j]=a[i][j]; for(i=0;i<N;i++) for(j=0;j<N;j++) temp[i][j]=0.0; for(i=0;i<N;i++) for(j=0;j<N;j++) temp[i][j]=-a1[i][j]/a1[i][i]; for(i=0;i<N;i++) b[i]/=a1[i][i]; for(k=0;k<N;k++) { for(i=1;i<=N;i++) { for(j=1;j<=N;j++) { if(j<=i-1) h+=temp[i-1][j-1]*x[j-1]; else if(j==i) h+=0 ; else h+=temp[i-1][j-1]*x[j-1]; } g=(h+b[i-1])*w; x[i-1]=(1-w)*x[i-1]+g; h=0.0;g=0.0; } } printf("超松驰法结果::\n"); printf("\n"); for(i=0;i<N;i++) printf("%f,",x[i]); /*输出超松驰法求解的特征向量*/ return 0; } float sign(float x) /*子函数符号函数*/ { if (x>0.0) return 1.0; else if (x<0.0) return -1.0; else return 0.0; } void guss(float a[N][N]) /*高斯法求解的子函数*/ { float a2[N][N]; int i,j,k; float l[N]={0,},x[N],h=0; float b[N]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4. 719345,-5.6784392}; for(i=0;i<N;i++)

for(j=0;j<N;j++) a2[i][j]=a[i][j]; for(k=0;k<N;k++) { for(j=k+1;j<N;j++) l[j]=a2[j][k]/a2[k][k]; for(i=k+1;i<N;i++) for(j=k+1;j<N;j++) a2[i][j]-=l[i]*a2[k][j]; for(i=k+1;i<N;i++) a2[i][k]=0; for(i=k+1;i<N;i++) b[i]-=l[i]*b[k]; } x[N-1]=b[N-1]/a2[N-1][N-1]; for(i=N-2;i>-1;i--) { for(j=i+1;j<N;j++) h+=a2[i][j]*x[j]; x[i]=(b[i]-h)/a2[i][i]; h=0; } printf("\n 高斯法结果:\n"); printf("\n"); for(i=0;i<N;i++) printf("%f,",x[i]); /*输出高斯法求解的特征向量*/ printf("\n"); return; } float duifen() /*对分法求取方程的解的子函数*/ { int i,j; float ss1=0,ss2=23,ss=0,h=0,g=0; float l,ll; h=f(ss2); for(i=0;1;i++) { ss1=ss; g=f(ss); ss=0.5*(ss+23); l=0.5*(ss1+23); if(g-h==1) break; } do{ if(f(ss)-h==0) { h=f(ss);ss2=ss; } else ss1=ss; ss=0.5*(ss1+ss2); ll=0.5*(ss1+ss2); if(fabs(ll-l)<0.01) break;

else l=ll; } while(1); return ll; } int f(float x) /*求取每次新的对分点的函数值*/ { int i,j,hh=0,n=9; float sa[9]; sa[0]=o[0]-x; if(sa[0]==0) sa[1]=e[0]*e[0]; else sa[1]=o[1]-x-e[0]*e[0]/sa[0]; for(i=2;i<n;i++) { if (sa[i-1]*sa[i-2]!=0) sa[i]=o[i]-x-e[i-1]*e[i-1]/sa[i-1]; else if(sa[i-2]==0) sa[i]=o[i]-x; if(sa[i-1]==0) sa[i]=-1; } for(i=0;i<n;i++) { if(sa[i]>=0) hh+=1; else hh+=0; } return hh; } float fanmi(float a[N][N]) /*反幂法求解方程的解的子函数*/ { float mk=0,t,t1; float z[N],yy[N],a3[N][N]; int i,j,k,m,tt; float temp[N][N]; float h=0; float l[N]={0,}; t1=duifen(); /*调用对分法求取方程的解*/ for(i=0;i<N;i++) z[i]=1; for(i=0;i<N;i++) for(j=0;j<N;j++) a3[i][j]=a[i][j]; for(i=0;i<N;i++) a3[i][i]-=t1; for(tt=0;tt<N;tt++) { for(m=0;m<N;m++) yy[m]=0; for(k=0;k<N;k++) { for(j=k+1;j<N;j++) l[j]=a3[j][k]/a3[k][k];

for(i=k+1;i<N;i++) for(j=k+1;j<N;j++) a3[i][j]-=l[i]*a3[k][j]; for(i=k+1;i<N;i++) a3[i][k]=0; for(i=k+1;i<N;i++) z[i]-=l[i]*z[k]; } yy[N-1]=z[N-1]/a3[N-1][N-1]; for(i=N-1;i>-1;i--) { for(j=i+1;j<N;j++) h+=a3[i][j]*yy[j]; yy[i]=(z[i]-h)/a3[i][i]; h=0; } for(m=0;m<N;m++) if(fabs(mk)<fabs(yy[m])) mk=yy[m]; for(m=0;m<N;m++) z[m]=yy[m]/mk; if (tt<(N-1)) mk=0; } t=t1+1/mk; printf("\n 反幂法的结果是:"); printf("\nd=%f\n",t); /*输出反幂法的解*/ for(m=0;m<N;m++) printf("%f,",z[m]); /*输出反幂法的解的特征向量*/ return 0; } main() { housholder(); sor(a,b); /*调用超松驰法求解*/ guss(a); /*调用高斯法求解*/ approx=duifen(a); /*调用对分法求解*/ printf("\n 对分法的结果是:"); printf("\nd=%f",approx); /*输出对分法近似解*/ fanmi(a); /*调用反幂法求解*/ printf("\n 第一题完"); } 四,运行结果 housholder 变换为: 12.3841 ,-4.8931 ,0.0000 ,0.0000 ,0.0000 ,0.0000 ,0.0000 ,0.0000 ,0.0000 , -4.8931 ,25.3984 ,6.4941 ,0.0000 ,-0.0000 ,0.0000 ,-0.0000 ,-0.0000 ,0.0000 , 0.0000 ,6.4941 ,20.6115 ,8.2439 ,0.0000 ,0.0000 ,0.0000 ,0.0000 ,0.0000 , 0.0000 ,0.0000 ,8.2439 ,23.4228 ,-13.8801 ,0.0000 ,0.0000 ,0.0000 ,0.0000 , 0.0000 ,-0.0000 ,0.0000 ,-13.8801 ,29.6983 ,4.5345 ,0.0000 ,-0.0000 ,0.0000 , 0.0000 ,0.0000 ,0.0000 ,0.0000 ,4.5345 ,16.0061 ,4.8814 ,0.0000 ,0.0000 , 0.0000 ,0.0000 ,0.0000 ,0.0000 ,0.0000 ,4.8814 ,26.0134 ,-4.5036 ,0.0000 ,

0.0000 ,-0.0000 ,0.0000 ,0.0000 ,-0.0000 ,0.0000 ,-4.5036 ,21.2540 ,4.5045 , 0.0000 ,0.0000 ,0.0000 ,0.0000 ,0.0000 ,0.0000 ,0.0000 ,4.5045 ,14.5341 , 超松驰法的结果:

1.073410,2.272582,-2.856601,2.292513,2.112164,-6.422586,1.357799,0.634258,-0.587 043, 高斯法的结果: 1.075800,2.275746,-2.855515,2.293098,2.112633,-6.423838,1.357920,0.634243,-0.587 267, 对分法的结果: d=21.916260 反幂法的结果: d=21.928205 0.157064,-0.306357,0.282195,0.285922,0.198636,0.533750,0.462840,1.000000,0.61185 1, 五.问题讨论 问题讨论 1.SOR 方法的矩阵形式为: X(m)=(E-ωL)-1((1-ω)E+ωR)x(m-1)+(E-ωL)-1ωg 若记 Lω=(E-ωL)-1((1-ω)E+ωR), 收敛的充要条件是 S(Lω)<1.且若 A 为对称正定阵, SOR 则当松弛因子ω满足 0<ω<2 时, 方法收敛。 SOR 此题中矩阵 B 是对称正定阵, 且是三对角的, 所以选择合适的松驰因子ω,收敛速度是很快的。 2.对分法简单可靠,数值稳定性较高,对于求少量几个特征值特别适宜. 3.反幂法是结合对分法使用的,所以只要近似值较恰当,其收敛速度是惊人的.只要迭代两次 就可得到较满意的结果.但在运用中需把一般矩阵化为 Hessebberg 阵,然后进行反迭代。

三. 用三次样条插值求函数值

一,程序要求 已知十点函数值及起点和终点的导数值,要求用三次样条插值求 f(4.563)及 f'(3)的近似值. 二,程序算法 其基本思想是对均匀分划的插值函数的构造,三次样条函数空间中不取 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).本题是属于问题一,由 s(xi)=yi,I=1,2,…N s'(x0)=y0',s'(xN)=yN'可得 S(xi)= cjΩ3(xi-xj/h)=yi S'(x0)=1/hcjΩ3'(x0-xj/h)=y'0 S'(xN)=1/hcjΩ3'(xN-xj/h)=y'N 通过列主元消去法解出 cj,在反代回去求出所需的函数值. 三,程序清单 #include<stdio.h> #include<math.h> #define N 12 float a[N][N]={{-2,-4,0,0,0,0,0,0,0,0,0,0}, {1,4,1,0,0,0,0,0,0,0,0,0}, {0,1,4,1,0,0,0,0,0,0,0,0}, {0,0,1,4,1,0,0,0,0,0,0,0}, {0,0,0,1,4,1,0,0,0,0,0,0}, {0,0,0,0,1,4,1,0,0,0,0,0}, {0,0,0,0,0,1,4,1,0,0,0,0}, {0,0,0,0,0,0,1,4,1,0,0,0}, {0,0,0,0,0,0,0,1,4,1,0,0}, {0,0,0,0,0,0,0,0,1,4,1,0}, {0,0,0,0,0,0,0,0,0,1,4,1}, {0,0,0,0,0,0,0,0,0,0,4,2}}; float x[N],c[N],sum=0; float y[N]={1,0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445, 2.1972246,2.3025851,0.1};/*以上是预处理和定义全局变量*/ main() {int i; float h=1.0,b[N],e[N],l[N],f[N],d[N],k[N],d1[N],d2[N],k1[N],k2[N],t,r; x[1]=1;x[0]=x[1]-h; for(i=1;i<N-1;i++) x[i+1]=1+i*h; b[0]=2*h*y[0]; b[N-1]=2*h*y[N-1]; for(i=1;i<N-1;i++) b[i]=6*y[i]; b[0]=b[0]-b[1];b[N-1]=b[N-1]+b[N-2]; e[0]=a[0][0];f[0]=b[0]; for(i=0;i<N-1;i++) {l[i+1]=a[i+1][i]/e[i]; e[i+1]=a[i+1][i+1]-l[i+1]*a[i][i+1]; f[i+1]=b[i+1]-l[i+1]*f[i];}

c[N-1]=f[N-1]/e[N-1]; for(i=N-2;i>=0;i--) c[i]=(f[i]-a[i][i+1]*c[i+1])/e[i]; for(i=0;i<N;i++) k[i]=4.563-x[i]; for(i=0;i<N;i++) {if(fabs(k[i])>=2) d[i]=0; else if(fabs(k[i])<=1) d[i]=(1/2.0)*pow(fabs(k[i]),3)-k[i]*k[i]+(2/3.0); else if((fabs(k[i])<2)&&(fabs(k[i])>1)) d[i]=(-1/6.0)*pow(fabs(k[i]),3)+k[i]*k[i]-2*fabs(k[i])+(4/3.0);} for(i=0;i<N;i++) sum=sum+c[i]*d[i]; printf(" f(4.563)的值为:\n"); printf("f(4.536)=%10.8f\n",sum); for(i=0;i<N;i++) k2[i]=3.5-x[i]; for(i=0;i<N;i++) {if(fabs(k2[i])>=1.5) d2[i]=0; else if(fabs(k2[i])<=(1/2.0)) d2[i]=-(k2[i]*k2[i])+(3/4.0); else if((fabs(k2[i])<=(3/2.0))&&(fabs(k2[i])>=(1/2.0))) d2[i]=(1/2.0)*k2[i]*k2[i]-(3/2.0)*fabs(k2[i])+(9/8.0);} sum=0; for(i=0;i<N;i++) sum=sum+c[i]*d2[i]; for(i=0;i<N;i++) k1[i]=2.5-x[i]; for(i=0;i<N;i++) {if(fabs(k1[i])>=(3/2)) d1[i]=0; else if(fabs(k1[i])<=1/2.0) d1[i]=-(k1[i]*k1[i])+(3/4.0); else if((fabs(k1[i])<=(3/2.0))&&(fabs(k1[i])>=(1/2.0))) d1[i]=(1/2.0)*k1[i]*k1[i]-(3/2.0)*fabs(k1[i])+(9/8.0);} t=0; for(i=0;i<N;i++) t=t+c[i]*d1[i]; r=sum-t; printf("f'(3)的值为: "); printf("f'(3)=%10.8f\n",r); } 四,运行结果 f(4.563)的值为: f(4.563)=1.51793242 f'(3)的值为:

f'(3)=0.33496869 五.问题讨论 样条插值效果比 Lagrange 插值好,由于样条插值不必经过所有点, 所以没有 Runge 现象.而且 插值函数比较光滑。 四. Newton 法解方程

一,程序要求 用牛顿法求方程 x7-28x4+14=0 在(0.1 1.9)中的近似根(初始近似值取为区间端点),迭代 6 次或误差小于 0.00001. 二, 程序算法 设函数在有限区间[a,b]上二阶导数存在,且满足条件 ⅰ f(a)f(b)<0 ⅱ f"(x)在区间[a,b]上不变号. ⅲ f'(x)≠ ⅳ |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… 所产生的迭代序列{xk}平方收敛于方程在区间[a,b]上的唯一解α 三,程序清单 #include "math.h" #include "stdio.h" main() { int i; float t=1e-6,b,x=0.9; for(i=0;1;i++) { b=x; x-=(pow(x,7)-28*pow(x,4)+14)/(7*pow(x,6)-112*pow(x,3)); if(fabs(x-b)<t) break; } printf("方程的解为:\n"); printf("\nx=%f",x); } 四,运行结果 方程的解为: x=0.845497 五.问题讨论: Newton 法收敛速度比较快,是平方收敛,但它是局部收敛。

五. 用 Romberg 法求积分

一.程序要求: 用 Romberg 算法求∫313xx1.4(5x+7)sinx2dx(允许误差ε=0.00001). 二.程序算法: 引入记号 Tm(k)称为 Romberge 序列,当|Ti(0)-Ti+1(0)|< ε就停止计算: T1(0)=(b-a)(f(a)+f(b))/2 T1(l)=1/2[T1(l-1)+(b-a)/2l-1 f(a+(2i-1)b-a/2)] Tm+1(k-1)=(4mTm(k)-Tm(k-1))/4m-1 三.程序清单: #include "math.h" int p(int x,int y) /*子函数求取 x 的值*/ { int i,A=1; for(i=0;i<y;i++) A*=x; return A; } float f(float x) /*子函数求取 3xx1.4(5x+7)sinx 2 的值*/ { float A; A=pow(3.0,x)*pow(x,1.4)*(5*x+7)*sin(x*x); return A; } main() { int i=1,j,k,n=10; float T1[10],a=1.0,b=3.0,s=0.0,h=0.00001; T1[0]=0.5*(b-a)*(f(a)+f(b)); /*调用子函数 F(x)*/ for(j=1;j<n-1;j++) { for(k=1;k<=p(2,j-1);k++) /*调用子函数 P(x,y)*/ s+=f(a+(2*k-1)*(b-a)/pow(2,j)); /*调用子函数 F(x)*/ T1[j]=0.5*(T1[j-1]+(b-a)*s/pow(2,j-1)); s=0.0; } T1[9]=(4*T1[1]-T1[0])/(float)3; while(fabs(T1[9]-T1[0])>h) {

T1[0]=T1[9]; for(j=1;j<n-1-i;j++) T1[j]=(pow(4,i)*T1[j+1]-T1[j])/(pow(4,i)-1); T1[9]=(pow(4,i+1)*T1[1]-T1[0])/(pow(4,i+1)-1); i++; } system("cls"); printf("积分结果为:\n"); printf("%f",T1[9]); } 四.计算结果: ∫313xx1.4(5x+7)sinx2dx=440.536011 五.问题讨论: Romberge 算法的优点是: 1.把积分化为代数运算,而实际上只需求 T1(i),以后用递推可得. 2.算法简单且收敛速度快,一般 4 或 5 次即能达到要求. Romberge 算法的缺点是: 1.对函数的光滑性要求较高, 2.计算新分点的值时,这些数值的个数成倍增加。

六. Runge_Kutta 求解微分方程组 一.程序要求 用定步长四阶法求解 y1'=1,y2'=y3 y3'=1000-1000y2-100y3 (y1(0)=y2(0)=y3(0)=0) h=0.0005,打印 yi(0.025),yi(0.045),yi(0.085),yi(0.1),(i=1,2,3) 二.程序算法 本题采用古典形式:令 N=4 yn+1=yn+h/6(k1+2k2+2k3+k4) k1=fn k2=f(tn+0.5h,yn+0.5hk1) k3=f(tn+0.5h,yn+0.5hk2) k4=f(tn+h,yn+hk3) 三.程序清单: #include"math.h" #define N 3 float y[N]={0,0,0},d[N][N+1]; float b[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<N;l++) for(i=0;i<=N;i++)/*计算每一个 d 值*/ d[l][i]=0; d[0][0]=h; d[1][0]=h*y[2]; d[2][0]=h*(1000-1000*y[1]-100*y[2]); d[0][1]=h; d[1][1]=h*(y[2]+(1/3.0)*d[0][0]); d[2][1]=h*(1000-1000*(y[1]+(1/3.0)*d[1][0])-100*(y[2]+(1/3.0)*d[1][0])); d[0][2]=h; d[1][2]=h*(y[2]+(1/3.0)*d[1][0]+d[1][1]); d[2][2]=h*(1000-1000*(y[1]+(1/3.0)*d[2][0]+d[2][1])-100*(y[2]+(1/3.0)*d[2][0]+ d[2][1])); d[0][3]=h; d[1][3]=h*(y[2]+d[1][0]-d[1][1]+d[1][2]); d[2][3]=h*(1000-1000*(y[1]+d[2][0]-d[2][1]+d[2][2])-100*(y[2]+d[2][0]-d[2][1]+ d[2][2])); for(i=0;i<N;i++) y[i]=y[i]+(d[i][0]+3*d[i][1]+3*d[i][2]+d[i][3])/8; if(n==50||n==90||n==170||n==200)/*得到所需的数据*/ {for(i=0;i<N;i++) b[j][i]=y[i]; j++;} n++; }while(n<=200); printf("计算结果为:\n"); for(j=0;j<N+1;j++) {for(i=0;i<N;i++) printf("y%d(%5.4f)=%.7f,",i+1,a[j],b[j][i]); printf("\b\n");} } 四.计算结果如下: y1(0.025)=0.0250000,y2(0.025)=0.1280306,y3(0.025)=7.7846808 y1(0.045)=0.0450000,y2(0.045)=0.2871802,y3(0.045)=7.7598104 y1(0.085)=0.0850000,y2(0.085)=0.5511296,y3(0.085)=5.3303890 y1(0.100)=0.1000001,y2(0.100)=0.6248569,y3(0.100)=4.4874177 五.问题讨论: Runge_Kutta 方的优点: 1.精度高,不必用别的方法求开始几点的函数值. 2.可根据 f'(t,y)变化的情况与需要的精度自动修改步长. 3.程序简单,存储量少. 4.方法稳定. Runge_Kutta 方的缺点: 每步要计算函数值 f(t,y)四次,在 f(t,y)较复杂时,工作量大, 且每一步缺乏可靠的 检查.


相关文章:
西南交大数值分析上机实习报告.doc
西南交大数值分析上机实习报告_理学_高等教育_教育专区 暂无评价|0人阅读|0次下载|举报文档西南交大数值分析上机实习报告_理学_高等教育_教育专区。研究生课程上机...
《数值分析》上机实习报告.doc
《数值分析》上机实习报告 - 数值分析程序设计作业 一.housholder 变
《数值分析》上机实习报告.doc
《数值分析》上机实习报告_学习总结_总结/汇报_实用文档。贵州大学《数值分析》上机实习报告 数值分析上机实验报告(C 语言) 姓班学院 名: 级:2006 级号:...
数值分析上机实习报告.doc
数值分析上机实习报告 - 数值分析上机实习报告 姓名:刘童超 学号:12011174 班级:结构 1 班 联系电话:15928948862 序言 数值分析,顾名思义,就是对于一些比较复杂的...
《数值分析》上机实习报告.doc
《数值分析》上机实习报告 - 数值分析上机实验报告 (C 语言) 姓班学院 名:
西南交通大学研究 数值分析上机实习报告2012.doc
西南交通大学研究 数值分析上机实习报告2012 - 数值分析上机实习报告要求 1
数值分析上机实习报告.doc
数值分析上机实习报告 - 研究生数值分析课程,雅可比迭代法和高斯-塞得尔迭代法的
《数值分析》上机实验报告.doc
《数值分析》上机实验报告 - 数值分析上机实验报告 《数值分析》上机实验报告 1
2014数值分析上机实习题.doc
2014数值分析上机实习题 - 1 信件主题请采用学号在前,姓名在后的格式,即形如: 201422222 张三数值实验报告 2. 报告文档采用.doc 或.pdf 格式, 也请按上述...
数值分析上机实习报告题目.doc
数值分析上机实习报告题目 - 数值分析上机实习报告要求 1.应提交一份完整的实习
谷根代数值分析--上机实习报告.doc
谷根代数值分析--上机实习报告_理学_高等教育_教育专区。《数值分析》谷根代、杨
数值分析上机报告.doc
数值分析上机报告 - 数值分析上机报告 姓学专 名: 号: 业: 联系电话: Matlab 是一种用于算法开发、 本次数值分析上机实习采用 Matlab 数学软件。 数据可视化、...
数值分析上机实习报告(西南交通大学).doc
数值分析上机实习报告(西南交通大学) - 数值分析上机实习报告 姓名: 学号: 专业: 电话: 大地测量学与测量工程 序 言 1. 所用程序语言:本次数值分析上机实习...
(最新版)哈工大-数值分析上机实验报告.doc
(最新版)哈工大-数值分析上机实验报告 - 创业计划,研究报告,项目建议书,项目建设,项目可行性研究报告,可行性研究报告,项目研究报告,项目设计
数值分析上机实习报告.doc
数值分析上机实习报告 上机实习一一、理论依据及算法 Householde
数值分析上机实习报告.doc
数值分析上机实习报告 - 指导教师: 指导教师: 姓学专名: 号: 业: 联系电话: 联系电话: 上海交通大学 数值分析实习报告 目序 录 言...
数值分析上机题答案.doc
数值分析上机题答案_院校资料_高等教育_教育专区。数值分析上机题的答案很详细 数值分析上机题姓名:武均 学号:142648 习题 1 17. (上机题)舍入误差与有效数 设...
数值分析上机实习报告10.doc
数值分析上机实习报告10 - 数值分析上机实习报告要求 1.应提交一份完整的实习
2014级硕士研究生数值分析上机实习报告.pdf
2014级硕士研究生数值分析上机实习报告 - 哈尔滨工业大学(威海)实验报告纸 2014 级硕士研究生数值分析上机实习 (第一次) 姓名: 学号: 学院: 实习题目: 分别用...
东南大学 数值分析上机题作业 MATLAB版.doc
东南大学 数值分析上机题作业 MATLAB版 - 2015.1.9 上机作业题报告 1.Chapter1 1.1 题目 设S = 1 =2 2 ?1,其精确值为 1 3 1 1 (...
更多相关标签: