数值分析上机作业
姓名：黄金金 学号：220112244
指导老师：江风
习题1
舍入误差与有效数 设∑=-=N
j N j S 2211，其精确值为)11123(21+--N N 。
（1）编制按从小到大的顺序
11131121222-+???+-+-=N S N ，计算N S 的通用程序。 （2）编制按从大到小的顺序
1211)1(111222-+???+--+-=
N N S N ，计算N S 的通用程序。
（3）按两种顺序分别计算642101010,,S S S ，并指出有效位数。（编制程序时用单精度）
（4）通过本上机题，你明白了什么？
c++代码：
#include
#include
void main(){
float S1,S2,S=0.0，N,j;
fstream myfile; //定义文件
myfile.open("chap0.txt",ios::in|ios::out);
for (int i=0;i<3;i++)
{
do //读入三组数据
{
cout<<"请输入N 值:";
cin>>N;
cout< myfile<<"当N = "< }while (N<2); S=(float) ((1.5-(1.0/N)-(1.0/(N+1)))/2); myfile<<"精确值S=:"< S1=0.0;S2=0.0; for (j=2;j<=N;j++) S1=S1+(float) (1.0/(j*j-1.0)); for (j=N;j>=2;j--) S2=S2+(float) (1.0/(j*j-1.0)); myfile<<"按从大到小的顺序求和结果S1=:"< myfile<<"按从小到大的顺序求和结果S2=:"< } myfile.close(); } 输出结果： 当N = 100 时 精确值S=:0.740049 按从大到小的顺序求和结果S1=:0.740049 有效位数6 按从小到大的顺序求和结果S2=:0.74005 有效位数5 当N = 10000时 精确值S=:0.7499 按从大到小的顺序求和结果S1=:0.749852 有效位数 4 按从小到大的顺序求和结果S2=:0.7499 有效位数 4 当N =1000000时 精确值S=:0.749999 按从大到小的顺序求和结果S1=:0.749852 有效位数 4 按从小到大的顺序求和结果S2=:0.749999 有效位数 6 体会：通过本题，我认识到了舍入误差对计算结果的影响，所以必须选择合适的计算方法， 尽量减少这种误差的影响。 习题2 NEWTON 迭代法 （1）给定初值0x 及容许误差ε，编制Newton 法解方程0)(=x f 根的通用程序。 （2）给定方程03/)(3=-=x x x f ，易知其有三个根3,0,3* 3*2*1==-=x x x 。 1．由Newton 方法的局部收敛性可知存在0>δ，当),(0δδ-∈x 时，Newton 迭代序列 收敛于根*2x 。试确定尽可能大的δ。 2.试取若干初值，观察当),1(),1,(),,(),,1(),1,(0∞-----∞∈δδδδx 时Newton 序列是 否收敛以及收敛于哪一个根。 （3）通过本上机题，你明白了什么？ c++代码： #include #include #define h 0.000001 float xo=0.0,x1,x2,e,a,b; float f(float x){ //原函数03 3 =-x x float f1=(float) (a*x*x*x+b*x); return f1; } float g(float x){ //导函数 float g1=(float) (3*a*x*x+b); return g1; } int test(float xk){ //求收敛于0的最大的δ int r,flag=0; float xkplus1; while (flag==0) { xkplus1 = xk - f(xk)/g(xk); if (xkplus1 == xk) flag=1; else xk = xkplus1; }; if (xk==0) r=1; else r=0; return r; } void main() { cout<<"请输入初值X0:"; cin>>x1; cout< cout<<"请输入容许误差E:"; cin>>e; cout< cout<<"请输入三次系数A,一次系数B:"; cin>>a>>b; cout< x2=(float) (x1-f(x1)/g(x1));// 牛顿迭代法求根 while (fabs(x2-x1)>=e) { x1=x2; x2=(x1-f(x1)/g(x1)); } cout<<"方程的根X 为:"< while (test(xo)==1){xo+=h;} // 求根x=0最大的迭代区间 cout<<"X=0的最大迭代区间为(-"< } 输出结果： X=0的最大迭代区间为 (-0.774597,0.774597); =ε0.0000000005 ∈0x (-∞,-1) ; 收敛于* 1x . ∈0x （-1,-δ）; 收敛于* 3x . ∈0x (-δ,δ); 收敛于* 2x . ∈0x (δ,1); 收敛于* 1x . ∈0x (1, ∞); 收敛于* 3x . 体会：牛顿迭代的优势在于收敛速度较快。 习题3 列主元三角分解法 对于某电路的分析，归结为求解线性方程组RI ＝V 。其中R ＝ ?? ? ? ? ? ? ? ? ? ? ? ?? ??????????????------------------2920000000227005000000413000000003047700000507573000090003079100000000103190000011093513000100001331 V T ＝ （-15，27，-23，0，-20，12，-7，7，10）T （1）编制n 阶线性方程组Ax ＝b 的列主元三角分解的通用程序； （2）用所编的程序解线性方程组RI ＝V ，并打印出向量，保留5位有效数； （3）本题编程中，你提高了哪些编程能力？ c++代码： #include float a[99][100],s[99],temp[100],x[100],y[100]; int n,i,j,max; int min(int x,int y){ //取x ，y 中的小者； int z; if (x else z=y; return z; } float sigma(int i,int j){ //对已经算出的行列作积和； float v=0.0; for (int q=0;q return v; } void main(){ cout<<"请输入阶数n:"; cin>>n; cout< for (i=0;i { cout<<"请输入增广矩阵的第"< for (j=0;j<=n;j++) cin>>a[i][j]; cout< } for (i=0;i max=0; //找第一个主元 for (i=1;i for (j=0;j<=n;j++) //交换行 { temp[j]=a[0][j]; a[0][j]=a[max][j]; a[max][j]=temp[j]; } for (i=1;i for (i=1;i { for (j=i;j max=i; //找出列主元 for (j=i+1;j if (max != i) { for (j=0;j<=n;j++) { temp[j]=a[i][j]; a[i][j]=a[max][j]; a[max][j]=temp[j]; } } for (j=i;j<=n;j++) a[i][j]=a[i][j]-sigma(i,j); for (j=i+1;j //求解 x[n-1]=a[n-1][n]/a[n-1][n-1]; for (i=n-2;i>=0;i--) { for (j=n-1,y[i]=0.0;j>i;j--) y[i]=x[j]*a[i][j]+y[i]; x[i]=(a[i][n]-y[i])/a[i][i]; } //输出 cout<<"X=["; for (i=0;i { cout< i++; if(i } cout<<"]T"< } 计算结果： x=[-0.28923,0.34544.-0.71281,-0.22061,-0.4304,0.15431,-0.057823,0.20105]T 体会：本题对矩阵的处理大量的使用了C语言中的For循环和数组，对数组和For循环的联合使用有个深刻体会。 习题4 3次样条插值函数 （1）编制求第一型3次样条插值函数的通用程序； 端点条件为y 0'＝0.8，y '＝0.2。用所编的程序求出车门的3次样条插值函数S（x）， 并打印出S（i＋0.5）（i＝0，1，…，9）。 c++代码： #include #include float x1[100],f1[100],f2[99],f3[98],m[100],a[100][101],x,d[100]; float c[99],e[99],h[99],u[99],w[99],y_0,y_n,arr,s; int i,j,k,n,q; void selectprint(float y) { if ((y>0)&&(y!=1)) cout<<"+"< else if (y==1) cout<<"+"; else if (y<0) cout< } void printY(float y){ if (y!=0) cout< } float calculation(float x){ for (j=1;j<=n;j++) if (x<=x1[j]) { s=(float)(f1[j-1]+c[j-1]*(x-x1[j-1])+m[j-1]/2.0*(x-x1[j-1])*(x-x1[j-1])+e[j-1]* (x-x1[j-1])*(x-x1[j-1])*(x-x1[j-1])); break; } return s; } void main() { do{ cout<<"请输入n值:"; cin>>n; if ((n>100)||(n<1)) cout<<"请重新输入整数(1..100):"< } while ((n>100)||(n<1)); cout<<"请输入Xi(i=0,1,...,"< for (i=0;i<=n;i++) cin>>x1[i]; cout< cout<<"请输入Yi(i=0,1,...,"< for (i=0;i<=n;i++) cin>>f1[i]; cout< cout<<"请输入第一型边界条件Y'0，Y'n:"; cin>>y_0>>y_n; cout< for (i=0;i for (i=1;i for (i=1;i for (i=0;i d[0]=6*(f2[0]-y_0)/h[0]; d[n]=6*(y_n-f2[n-1])/h[n-1]; for (i=0;i<=n;i++) for (j=0;j<=n;j++) if (i==j) a[i][j]=2; else a[i][j]=0; a[0][1]=1; a[n][n-1]=1; for (i=1;i { a[i][i-1]=u[i]; a[i][i+1]=w[i]; } for (i=0;i<=n;i++) a[i][n+1]=d[i]; for (i=1;i<=n;i++) //用追赶法解方程,得M[i] { arr=a[i][i-1]; for (j=0;j<=n+1;j++) a[i][j]=a[i][j]-a[i-1][j]*arr/a[i-1][i-1]; } m[n]=a[n][n+1]/a[n][n]; for (i=n-1;i>=0;i--) m[i]=(a[i][n+1]-a[i][i+1]*m[i+1])/a[i][i]; for (i=0;i c[i]=(float) (f2[i]-(1.0/3.0*m[i]+1.0/6.0*m[i+1])*h[i]); for (i=0;i e[i]=(m[i+1]-m[i])/(6*h[i]); for (i=0;i { cout<<"X属于区间["< cout<<"S(x)="; printY(f1[i]); if (c[i]!=0) { selectprint(c[i]); cout<<"(x"; printY(-x1[i]); cout<<")"; } if (m[i]!=0) { selectprint((float)(m[i]/2.0)); for (q=0;q<2;q++) { cout<<"(x"; printY(-x1[i]); cout<<")"; } } if (e[i]!=0) { selectprint(e[i]); for (q=0;q<3;q++) { cout<<"(x"; printY(-x1[i]); cout<<")"; } } cout< } cout<<"针对本题，计算S(i+0.5)(i=0,1,...,9):"< for (i=0;i<10;i++) { if ((i+0.5<=x1[n])&&(i+0.5>=x1[0])) { calculation((float)(i+0.5)); cout<<"S("< } else cout< }; cout< } 输出结果： X属于区间[0,1]时S(x)=2.51+0.8(x)-0.0014861(x)(x)-0.00851395(x)(x)(x) X属于区间[1,2]时 S(x)=3.3+0.771486(x-1)-0.027028(x-1)(x-1)-0.00445799(x-1)(x-1)(x-1) X属于区间[2,3]时 S(x)=4.04+0.704056(x-2)-0.0404019(x-2)(x-2)-0.0036543(x-2)(x-2)(x-2) X属于区间[3,4]时 S(x)=4.7+0.612289(x-3)-0.0513648(x-3)(x-3)-0.0409245(x-3)(x-3)(x-3) X属于区间[4,5]时 S(x)=5.22+0.386786(x-4)-0.174138(x-4)(x-4)+0.107352(x-4)(x-4)(x-4) X属于区间[5,6]时 S(x)=5.54+0.360567(x-5)+0.147919(x-5)(x-5)-0.268485(x-5)(x-5)(x-5) X属于区间[6,7]时 S(x)=5.78-0.149051(x-6)-0.657537(x-6)(x-6)+0.426588(x-6)(x-6)(x-6) X属于区间[7,8]时 S(x)=5.4-0.184361(x-7)+0.622227(x-7)(x-7)-0.267865(x-7)(x-7)(x-7) X属于区间[8,9]时 S(x)=5.57+0.256496(x-8)-0.181369(x-8)(x-8)+0.0548728(x-8)(x-8)(x-8) X属于区间[9,10]时 S(x)=5.7+0.058376(x-9)-0.0167508(x-9)(x-9)+0.0583752(x-9)(x-9)(x-9) 计算S(i+0.5)(i=0,1,...,9): S(0.5)=2.90856 S(1.5)=3.67843 S(2.5)=4.38147 S(3.5)=4.98819 S(4.5)=5.38328 S(5.5)=5.7237 S(6.5)=5.59441 S(7.5)=5.42989 S(8.5)=5.65976 S(9.5)=5.7323 体会：精确度更高 习题5 常微分方程初值问题数值解 (1)编制RK4方法的通用程序； (2)编制AB4方法的通用程序（由RK 4提供初值）； (3)编制AB4-AM4预测校正方法的通用程序（由RK 4提供初值）； (4)编制带改进的AB4-AM4预测校正方法通用程序（由RK 4提供初值）； (5)对于初值问题 ???=≤≤-=3)0()5.10('22y x y x y 取步长1.0=h ，应用（1）～（4）中的四种方法进行计算，并将计算结果和精确解 )1/(3)(3x x y +=作比较； (6)通过本上机题，你能得到哪些结论？ 源程序： #include #include const double h=0.1; /*步长，可改为要求的值*/ const double a=0; /*区间[a,b]左边界值a ，可改改为要求的值*/ const double b=1.5; /*区间[a,b]右边界值b ，可改改为要求的值*/ const double y_0=3.0; /*初值，可改改为要求的值*/ double f(double,double); /*返回f(x,y)函数值*/ double RK4(int); /*RK4方法的通用程序*/ double AB4(int); /*AB4方法的通用程序*/ double AB4_AM4(int); /*AB4-AM4方法的通用程序*/ double improved_AB4_AM4(int); /*带改进的AB4-AM4方法的通用程序*/ void main(){ int i; cout.precision(7); cout<<"\n 常微分方程初值问题的通用程序\n"< double n=(b-a)/h; int m=int(n); /*区间[a,b]上的等分数 */ cout<<"xi"< cout< cout.precision(7); for(i=0;i<=3;i++){ double y=3/(1+(a+i*h)*(a+i*h)*(a+i*h)); if(i==0){ cout< cout< } else{ cout< cout< } } for(i=4;i<=m;i++){ double y=3/(1+(a+i*h)*(a+i*h)*(a+i*h)); cout< cout< } cout< } double f(double x,double y){ /*返回f(x,y)函数值*/ return -x*x*y*y; } double RK4(int m){ /*RK4方法的通用程序*/ int i; double y,y1=y_0; double k1,k2,k3,k4; for(i=0;i k1=f(a+i*h,y1); /* a+i*h即xi */ k2=f(a+i*h+0.5*h,y1+0.5*h*k1); k3=f(a+i*h+0.5*h,y1+0.5*h*k2); k4=f(a+i*h+h,y1+h*k3); y=y1+h/6*(k1+2*k2+2*k3+k4); y1=y; } return y; } double AB4(int m){ /*AB4方法的通用程序*/ int i; double y; double k1,k2,k3,k4; double y0=y_0,y1=RK4(1),y2=RK4(2),y3=RK4(3); for(i=3;i k1=f(a+i*h,y3); k2=f(a+(i-1)*h,y2); k3=f(a+(i-2)*h,y1); k4=f(a+(i-3)*h,y0); y=y3+h/24*(55*k1-59*k2+37*k3-9*k4); y0=y1,y1=y2,y2=y3,y3=y; } return y; } double AB4_AM4(int m){ /*AB4-AM4方法的通用程序*/ int i; double y,yp; double k5,k1,k2,k3,k4; double y0=y_0,y1=RK4(1),y2=RK4(2),y3=RK4(3); for(i=3;i k2=f(a+i*h,y3); k3=f(a+(i-1)*h,y2); k4=f(a+(i-2)*h,y1); k5=f(a+(i-3)*h,y0); yp=y3+h/24*(55*k2-59*k3+37*k4-9*k5); k1=f(a+(i+1)*h,yp); y=y3+h/24*(9*k1+19*k2-5*k3+k4); y0=y1,y1=y2,y2=y3,y3=y; } return y; } double improved_AB4_AM4(int m){ /*带改进的AB4-AM4方法的通用程序*/ int i; double yc,yp,y; double k5,k1,k2,k3,k4; double y0=y_0,y1=RK4(1),y2=RK4(2),y3=RK4(3); for(i=3;i k2=f(a+i*h,y3); k3=f(a+(i-1)*h,y2); k4=f(a+(i-2)*h,y1); k5=f(a+(i-3)*h,y0); yp=y3+h/24*(55*k2-59*k3+37*k4-9*k5); k1=f(a+(i+1)*h,yp); yc=y3+h/24*(9*k1+19*k2-5*k3+k4); y=251.0/270*yc+19.0/270*yp; y0=y1,y1=y2,y2=y3,y3=y; } return y; } 运行结果： 常微分方程初值问题的通用程序 xi 精确解y(xi) RK4 AB4 AB4_AM4 im_AB4_AM4 0 3 3 3 3 3 0.1 2.997003 2.997003 2.997003 2.997003 2.997003 0.2 2.97619 2.97619 2.97619 2.97619 2.97619 0.3 2.92113 2.921129 2.921129 2.921129 2.921129 0.4 2.819549 2.819547 2.818389 2.819678 2.819588 0.5 2.666667 2.666663 2.664672 2.666876 2.666713 0.6 2.467105 2.4671 2.465203 2.467252 2.467097 0.7 2.233805 2.233799 2.233079 2.233731 2.233682 0.8 1.984127 1.984123 1.984951 1.983787 1.983885 0.9 1.735107 1.735107 1.737043 1.734607 1.734808 1 1.5 1.500006 1.502195 1.499516 1.499732 1.1 1.287001 1.287013 1.288763 1.286657 1.286821 1.2 1.099707 1.099722 1.100724 1.099533 1.099622 1.3 0.9383797 0.9383975 0.9387105 0.9383425 0.9383673 1.4 0.8012821 0.8013004 0.8011349 0.8013274 0.8013114 1.5 0.6857143 0.6857321 0.6853346 0.6857961 0.6857604 Press any key to continue 习题6 抛物线方程Crank-Nicolson 格式 1）编制用Crank-Nicolson 格式求抛物线方程 ??? a -t u/2u/?2x = f(x,t) (0 u(x,0) = )x (? (0)1x ≤≤ u(0,t) = ()x α ，u （1,t ）=()t β (0 数值解的通用程序。 2）就a=1, f(x,t)=0,()x ?=exp(x),()t α=exp(t),()t β=exp(1+t),M=40,N=40,输入点（0.2， 1.0），（0.4，1.0），（0.6，1.0），（.8，1.0）4点处u （x,t ）的近似值。 3） 已知所给方程的精确解为u(x,t)=exp(x+t),将步长反复二分，从（0.2，1.0），（0.4， 1.0），（0.6，1.0），（0.8，1.0）4点处精确解与数值解的误差观察当步 长缩小一半时， 误差以什么规律缩小。 算法： 本题选择空间步长h=0.025,时间步长ζ=0.025 1.置初值； 2.通过点(i-1,k),(i-1,k+1),(i,k),(i+1,k),(i+1,k+1)求解点(I,k+1)的值； 3.输出结果 原程序： #include #include float h=0.025,k=0.025; int m=40;int n=40; float y[40][40],r=a*k/(h*h); void Input() { int i,j; cout<<"Loading Input Data..."< for(i=0;i { for(j=0;j if (i==j) a[i][j]=1+r; for(j=0;j if ((j=i+1)||(i=j+1)) a[i][j]=-r/2; } } int main() { Input(); //read data int r,i; for(k=0;k<(m-1);k++) { int select(); //select main element r=select(); void exchange(int g); exchange(r); //exchange void analyze(); analyze(); //analyze } void ret(); ret(); // replace back cout<<"The solution vector is below:"< cout<<"x["< } int select() { int f,t; float max; f=k; float si(int u,int v); max=float(fabs(si(k,k))); for(t=(k+1);t<(m-1);t++) { if(max { max=float(fabs(si(t,k))); f=t; } } return f; } float si(int u,int v) { float sum=0; int q; for(q=0;q sum+=a[u][q]*a[q][v]; sum=a[u][v]-sum; return sum; } void exchange(int g) { int t; float temp; for(t=0;t { temp=a[k][t]; a[k][t]=a[g][t]; a[g][t]=temp; } } void analyze() { int t; float si(int u,int v); for(t=k;t a[k][t]=si(k,t); for(t=(k+1);t a[t][k]=(float)(si(t,k)/a[k][k]); } void ret() { int t,z;float sum; x[m-1]=(float)a[m-1][m]/a[m-1][m-1]; for(t=(m-2);t>-1;t--) { sum=0; for(z=(t+1);z sum+=a[t][z]*x[z]; x[t]=(float)(a[t][m]-sum)/a[t][t]; } } 运行结果: ===Numerical Solution of Partial Differential Equations=== NumericalResult AccurateValue Errors （0.2,1.0） 3.31947 3.32012 0.00045 （0.4,1.0） 4.05432 4.05520 0.00088 （0.6,1.0） 4.95238 4.95303 0.00065 （0.8,1.0） 4.04897 6.04965 0.00067 分析与结论: Crank-Nicolson格式，步长越小，精度高。 昆明理工大学工科研究生《数值分析》上机实验 学院：材料科学与工程学院 专业：材料物理与化学 学号：2011230024 姓名: 郑录 任课教师：胡杰 P277-E1 1.已知矩阵A= 10787 7565 86109 75910 ?? ?? ?? ?? ?? ??,B= 23456 44567 03678 00289 00010 ?? ?? ?? ?? ?? ?? ?? ?? ,错误！未找到引用源。 = 11/21/31/41/51/6 1/21/31/41/51/61/7 1/31/41/51/61/71/8 1/41/51/61/71/81/9 1/51/61/71/81/91/10 1/61/71/81/91/101/11?????????????????? （1）用MA TLAB函数“eig”求矩阵全部特征值。 （2）用基本QR算法求全部特征值（可用MA TLAB函数“qr”实现矩阵的QR分解）。解：MA TLAB程序如下： 求矩阵A的特征值： clear; A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; E=eig(A) 输出结果： 求矩阵B的特征值： clear; B=[2 3 4 5 6;4 4 5 6 7;0 3 6 7 8;0 0 2 8 9;0 0 0 1 0]; E=eig(B) 输出结果： 求矩阵错误！未找到引用源。的特征值： clear; 错误！未找到引用源。=[1 1/2 1/3 1/4 1/5 1/6; 1/2 1/3 1/4 1/5 1/6 1/7; 1/3 1/4 1/5 1/6 1/7 1/8; 1/4 1/5 1/6 1/7 1/8 1/9;1/5 1/6 1/7 1/8 1/9 1/10; 1/6 1/7 1/8 1/9 1/10 1/11]; E=eig(错误！未找到引用源。) 输出结果： （2）A= 10 7877565861097 5 9 10 第一步：A0=hess(A);[Q0,R0]=qr(A0);A1=R0*Q0 返回得到： 第二部：[Q1,R1]=qr(A1);A2=R1*Q1 数值分析上机实验报告 选题：曲线拟合的最小二乘法 指导老师： 专业： 学号： 姓名： 课题八曲线拟合的最小二乘法 一、问题提出 从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量y 与时间t 的拟合曲线。 二、要求 1、用最小二乘法进行曲线拟合； 2、近似解析表达式为()33221t a t a t a t ++=?； 3、打印出拟合函数()t ?，并打印出()j t ?与()j t y 的误差，12,,2,1 =j ； 4、另外选取一个近似表达式，尝试拟合效果的比较； 5、*绘制出曲线拟合图*。 三、目的和意义 1、掌握曲线拟合的最小二乘法； 2、最小二乘法亦可用于解超定线代数方程组； 3、探索拟合函数的选择与拟合精度间的关系。 四、计算公式 对于给定的测量数据(x i ,f i )(i=1,2,…，n )，设函数分布为 ∑==m j j j x a x y 0)()(? 特别的，取)(x j ?为多项式 j j x x =)(? (j=0, 1,…，m ) 则根据最小二乘法原理，可以构造泛函 ∑∑==-=n i m j i j j i m x a f a a a H 1 10))((),,,(? 令 0=??k a H (k=0, 1,…，m ) 则可以得到法方程 ???? ??????? ?=????????????????????????),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ????????????????????? 求该解方程组，则可以得到解m a a a ,,,10 ，因此可得到数据的最小二乘解 ∑=≈m j j j x a x f 0)()(? 曲线拟合：实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。 五、结构程序设计 在程序结构方面主要是按照顺序结构进行设计，在进行曲线的拟合时，为了进行比较，在程序设计中，直接调用了最小二乘法的拟合函数polyfit ，并且依次调用了plot 、figure 、hold on 函数进行图象的绘制，最后调用了一个绝对值函数abs 用于计算拟合函数与原有数据的误差，进行拟合效果的比较。 大作业 三 1. 给定初值 0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用 程序. 解：Matlab 程序如下： 函数m 文件：fu.m function Fu=fu(x) Fu=x^3/3-x; end 函数m 文件：dfu.m function Fu=dfu(x) Fu=x^2-1; end 用Newton 法求根的通用程序Newton.m clear; x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0) while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0) 东南大学数值分析上机作业 汇总 -标准化文件发布号：（9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII 数值分析上机报告 院系： 学号： 姓名： 目录 作业1、舍入误差与有效数 (1) 1、函数文件cxdd.m (1) 2、函数文件cddx.m (1) 3、两种方法有效位数对比 (1) 4、心得 (2) 作业2、Newton迭代法 (2) 1、通用程序函数文件 (3) 2、局部收敛性 (4) （1）最大δ值文件 (4) （2）验证局部收敛性 (4) 3、心得 (6) 作业3、列主元素Gauss消去法 (7) 1、列主元Gauss消去法的通用程序 (7) 2、解题中线性方程组 (7) 3、心得 (9) 作业4、三次样条插值函数 (10) 1、第一型三次样条插值函数通用程序： (10) 2、数据输入及计算结果 (12) 作业1、舍入误差与有效数 设∑ =-=N j N j S 2 2 11 ，其精确值为?? ? ??---1112321N N . （1）编制按从小到大的顺序1 1 131121222-? ??+-+-=N S N ，计算N S 的通用程序； （2）编制按从大到小的顺序()1 21 11111222-???+--+-=N N S N ，计算N S 的通用程序； （3）按两种顺序分别计算642101010,,S S S ,并指出有效位数； （4）通过本上机你明白了什么？ 程序： 1、函数文件cxdd.m function S=cxdd(N) S=0; i=2.0; while (i<=N) S=S+1.0/(i*i-1); i=i+1; end script 运行结果（省略>>）： S=cxdd(80) S= 0.737577 2、函数文件cddx.m function S=cddx (N) S=0; for i=N:-1:2 S=S+1/(i*i-1); end script 运行结果（省略>>）： S=cddx(80) S= 0.737577 3、两种方法有效位数对比 2015.1.9 上机作业题报告 USER 1.Chapter1 1.1题目 设S N = 1 j 2?1 N j =2 ，其精确值为 )1 1123(21+--N N 。 （1）编制按从大到小的顺序1 1 1311212 22-+??+-+-=N S N ，计算S N 的通用程序。 （2）编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 （3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） （4）通过本次上机题，你明白了什么？ 1.2程序 1.3运行结果 1.4结果分析 按从大到小的顺序，有效位数分别为：6，4，3。 按从小到大的顺序，有效位数分别为：5，6，6。 可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。因此，采取从小到大的顺序累加得到的结果更加精确。 2.Chapter2 2.1题目 （1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。 （2）给定方程03 )(3 =-=x x x f ,易知其有三个根3,0,3321= *=*-=*x x x ○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收敛于根x2*。试确定尽可能大的δ。 ○2试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。 （3）通过本上机题，你明白了什么？ 2.2程序 大作业 三 1. 给定初值 0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用程序. 解：Matlab 程序如下： 函数m 文件：fu.m function Fu=fu(x) Fu=x^3/3-x; end 函数m 文件：dfu.m function Fu=dfu(x) Fu=x^2-1; end 用Newton 法求根的通用程序Newton.m clear; x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0) while flag==1 sigma=k*eps; x0=sigma; k=k+1; m=0; flag1=1; while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0) end end fprintf('最大的sigma 值为：%f\n',sigma); 2．求下列方程的非零根 5130.6651()ln 05130.665114000.0918 x x f x x +?? =-= ?-???解： Matlab 程序为： （1）主程序 clear clc format long x0=765; N=100; errorlim=10^(-5); x=x0-f(x0)/subs(df(),x0); n=1; Matlab作业3（数值分析） 机电工程学院（院、系）专业班组 学号姓名实验日期教师评定 1.计算多项式乘法(x2+2x+2)(x2+5x+4)。 答： 2. （1）将(x-6)(x-3)(x-8)展开为系数多项式的形式。（2）求解在x=8时多项 式(x-1)(x-2) (x-3)(x-4)的值。 答：（1） （2） 3. y=sin(x)，x从0到2π，?x=0.02π，求y的最大值、最小值、均值和标准差。 4.设ｘ＝[0.00.30.8 1.1 1.6 2.3]'，y=[0.500.82 1.14 1.25 1.35 1.40]'，试求二次多项式拟合系数，并据此计算x1=[0.9 1.2]时对应的y1。解：x=[0.0 0.3 0.8 1.1 1.6 2.3]'; %输入变量数据x y=[0.50 0.82 1.14 1.25 1.35 1.40]'; %输入变量数据y p=polyfit(x,y,2) %对x,y用二次多项式拟合,得到系数p x1=[0.9 1.2]; %输入点x1 y1=polyval(p,x1) %估计x1处对应的y1 p = -0.2387 0.9191 0.5318 y1 = a) 1.2909 5.实验数据处理：已知某压力传感器的测试数据如下表 p为压力值，u为电压值，试用多项式 d cp bp ap p u+ + + =2 3 ) ( 来拟 合其特性函数，求出a,b,c,d，并把拟合曲线和各个测试数据点画在同一幅图上。解： >> p=[0.0,1.1,2.1,2.8,4.2,5.0,6.1,6.9,8.1,9.0,9.9]; u=[10,11,13,14,17,18,22,24,29,34,39]; x=polyfit(p,u,3) %得多项式系数 t=linspace(0,10,100); y=polyval(x,t); %求多项式得值 plot(p,u,'*',t,y,'r') %画拟和曲线 x = 0.0195 -0.0412 1.4469 9.8267 ￥ 数值分析思考题1 1、讨论绝对误差（限）、相对误差（限）与有效数字之间的关系。 2、相对误差在什么情况下可以用下式代替 3、查阅何谓问题的“病态性”，并区分与“数值稳定性”的不同点。 4、取 ，计算 ，下列方法中哪种最好为什么（1）(3 3-，（2）(2 7-，（3） ()3 1 3+ ，（4） ()6 1 1 ，（5）99- , 数值实验 数值实验综述：线性代数方程组的解法是一切科学计算的基础与核心问题。求解方法大致可分为直接法和迭代法两大类。直接法——指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法，因此也称为精确法。当系数矩阵是方的、稠密的、无任何特殊结构的中小规模线性方程组时，Gauss消去法是目前最基本和常用的方法。如若系数矩阵具有某种特殊形式，则为了尽可能地减少计算量与存储量，需采用其他专门的方法来求解。 Gauss消去等同于矩阵的三角分解，但它存在潜在的不稳定性，故需要选主元素。对正定对称矩阵，采用平方根方法无需选主元。方程组的性态与方程组的条件数有关，对于病态的方程组必须采用特殊的方法进行求解。 数值计算方法上机题目1 1、实验1. 病态问题 实验目的： 算法有“优”与“劣”之分，问题也有“好”和“坏”之别。所谓坏问题就是问题本身的解对数据变化的比较敏感，反之属于好问题。希望读者通过本实验对此有一个初步的体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 $ r e x x e x x ** * ** - == 141 . ≈)61 目录 第一章非线性方程求根 (3) 1.1迭代法 (3) 1.2牛顿法 (4) 1.3弦截法 (5) 1.4二分法 (6) 第二章插值 (7) 2.1线性插值 (7) 2.2二次插值 (8) 2.3拉格朗日插值 (9) 2.4分段线性插值 (10) 2.5分段二次插值 (11) 第三章数值积分 (13) 3.1复化矩形积分法 (13) 3.2复化梯形积分法 (14) 3.3辛普森积分法 (15) 3.4变步长梯形积分法 (16) 第四章线性方程组数值法 (17) 4.1约当消去法 (17) 4.2高斯消去法 (18) 4.3三角分解法 (20) 4.4雅可比迭代法 (21) 4.5高斯—赛德尔迭代法 (23) 第五章常积分方程数值法 (25) 5.1显示欧拉公式法 (25) 5.2欧拉公式预测校正法 (26) 5.3改进欧拉公式法 (27) 5.4四阶龙格—库塔法 (28) 数值计算方法 第一章非线性方程求根 1.1迭代法 程序代码： Private Sub Command1_Click() x0 = Val(InputBox("请输入初始值x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = (Exp(2 * x0) - x0) / 5 If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求f(x)=e2x-6x=0在x=0.5附近的根（ep=10-10） 1.2牛顿法 程序代码： Private Sub Command1_Click() b = Val(InputBox("请输入被开方数x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = x0 - (x0 ^ 2 - b) / (2 * b) If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求56的值。（ep=10-10） 1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+-=-，从0I 的几个近似值出发，计算20I ； 解：易得：0I =ln6-ln5=0.1823, 程序为： I=0.182; for n=1:20 I=(-5)*I+1/n; end I 输出结果为：20I = -3.0666e+010 （2） 粗糙估计20I ，用n I I n n 51 5111+- =--，计算0I ； 因为 0095.05 6 0079.01020 201 020 ≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2 1 20=+= I 程序为：I=0.0087; for n=1:20 I=(-1/5)*I+1/(5*n); end I 0I = 0.0083 （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 首先分析两种递推式的误差；设第一递推式中开始时的误差为000I I E '-=，递推过程的舍入误差不计。并记n n n I I E '-=，则有01)5(5E E E n n n -==-=-Λ。因为=20E 20020)5(I E >>-，所此递推式不可靠。而在第二种递推式中n n E E E )5 1(5110-==-=Λ，误差在缩小， 所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制， 即算法是否数值稳定。 2． 求方程0210=-+x e x 的近似根，要求4 1105-+?<-k k x x ，并比较计算量。 （1） 在[0，1]上用二分法； 程序：a=0;b=1.0; while abs(b-a)>5*1e-4 c=(b+a)/2; 数值分析报大作业 班级：铁道2班 专业：道路与铁道工程 姓名：蔡敦锦 学号：13011260 一、序言 该数值分析大作业是通过C语言程序编程在Microsoft Visual C++ 6.0编程软件上运行实现的。本来是打算用Matlab软间来计算非线性方程的根的。学习Matlab也差不多有一个多月了，感觉自己编程做题应该没什么问题了；但是当自己真心的去编程、运行时才发现有很多错误，花了一天时间修改、调试程序都没能得到自己满意的结果。所以，我选择了自己比较熟悉的C程序语言来编程解决非线性的求值问题，由于本作业是为了比较几种方法求值问题的收敛速度和精度的差异，选择了一个相对常见的非线性函数来反映其差异，程序运行所得结果我个人比较满意。编写C语言，感觉比较上手，程序出现问题也能比较熟练的解决。最终就决定上交一份C程序语言编程的求值程序了！ 二、选题 本作业的目的是为了加深对非线性方程求根方法的二分法、简单迭代法、、牛顿迭代法弦截法等的构造过程的理解；能将各种方法的算法描述正确并且能够改编为程序并在计算机上实现程序的正确合理的运行，能得到自己满意的结果，并且能调试修改程序中可能出现的问题和程序功能的增减修改。本次程序是为了比较各种方法在求解同一非线性方程根时，在收敛情况上的差异。 为了达到上面的条件我选择自己比较熟悉的语言—C语言来编程，所选题目为计算方程f(x)=x3-2x-5=0在区间[2,3]内其最后两近似值的差的绝对值小于等于5 ?的根的几种方法的比较。 110- 本文将二分法、牛顿法、简单迭代法、弦截法及加速收敛法这五种方法在同一个程序中以函数调用的方式来实现，比较简洁明了，所得结果能很好的比较，便于分析；发现问题和得出结论。 高等数值分析上机作业 目录 上机作业1 (1) 上机作业2 (5) 上机作业3 (10) 上机作业4 (13) 上机作业5 (16) 上机作业6 (19) 上机作业7 (20) 上机作业8 (29) 第8章 函数逼近与曲线拟合 上机作业1： 最佳平方逼近 8-11.设()[]1,1,arcsin -∈=x x x f ， （1） 在{}32,,,1x x x span =φ中求()x f 的最佳平方逼近多项式； （2） 在{})(),(),(),(3210x T x T x T x T span =φ中求()x f 的最佳平方逼近多项式。 解：（1） 基于幂函数的最佳平方逼近 简单原理： 对于],[)(b a C x f ∈及一个线性无关函数组的集合 {}],,[)(,),(),(10b a C x x x span n ?=???φ 若存在,φ∈*S 使得 ()dx x S x f x x S x f x S x f b a x S x S ?-=-=-∈∈* 2)(2 2)(2 2 )]()([min )()(min )()(ρφ φ ，则称()x S *是 ()x f 在子集[]b a ,?φ中的最佳平方逼近函数。 取(),,,1,0,n j x x j j ==?就有{}n x x x span ,,,,12 =φ。对于任意的()φ∈x S ，有()∑==n j j j x a x S 0，()x S 为次数n ≤的多项式。 令)(x f 在},,,1{32x x x span =φ中的最佳平方逼近函数为 φ∈+++=3 *2**1*0*3 2)(x a x a x a a x S 通过求解法方程 ???? ? ? ? ??=??????? ????????? ??),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),() ,(321032103323130 3322212023121110130201000????????????????????????????????????f f f f a a a a 其中.arcsin )(,)(,)(,)(,1)(332210x x f x x x x x x x =====???? 数值计算方法上机题目1 1、实验1. 病态问题 实验目的： 算法有“优”与“劣”之分，问题也有“好”和“坏”之别。所谓坏问题就是问题本身的解对数据变化的比较敏感，反之属于好问题。希望读者通过本实验对此有一个初步的体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出： 考虑一个高次的代数多项式 ∏=-= ---=20 1)()20)...(2)(1()(k k x x x x x p （E1-1） 显然该多项式的全部根为l ，2，…，20，共计20个，且每个根都是单重的（也称为简 单的）。现考虑该多项式方程的一个扰动 0)(19 =+x x p ε (E1-2) 其中ε是一个非常小的数。这相当于是对（E1-1）中19 x 的系数作一个小的扰动。我们希望比较（E1-1）和（E1-2）根的差别，从而分析方程（E1-1）的解对扰动的敏感性。 实验内容： 为了实现方便，我们先介绍两个 Matlab 函数：“roots ”和“poly ”，输入函数 u ＝roots （a ） 其中若变量a 存储1+n 维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,...,,+n a a a ，则输出u 的各分量是多项式方程 0...1121=++++-n n n n a x a x a x a 的全部根，而函数 b=poly(v) 的输出b 是一个n ＋1维变量，它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“Poly ”是两个互逆的运算函数. ve=zeros(1,21); ve(2)=ess; roots(poly(1:20))+ve) 上述简单的Matlab 程序便得到（E1-2）的全部根，程序中的“ess ”即是（E1-2）中的ε。 实验要求： （1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。如果扰动项的系数ε很小，我们自然感觉（E1-1）和(E1-2)的解应当相差很小。计算中你有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？ （2）将方程（E1-2）中的扰动项改成18 x ε或其他形式，实验中又有怎样的现象出现？ 《数值分析》计算实习题目 第一题： 1. 算法设计方案 （1）1λ，501λ和s λ的值。 1)首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，数值小的为所求最小特征值λ1，数值大的为是所求最大特征值λ501。 2)使用反幂法求λs ，其中需要解线性方程组。因为A 为带状线性方程组，此处采用LU 分解法解带状方程组。 （2）与140k λλμλ-5011=+k 最接近的特征值λik 。 通过带有原点平移的反幂法求出与数k μ最接近的特征值 λik 。 （3）2cond(A)和det A 。 1）1=n λλ2cond(A)，其中1λ和n λ分别是按模最大和最小特征值。 2）利用步骤（1）中分解矩阵A 得出的LU 矩阵，L 为单位下三角阵,U 为上三角阵，其中U 矩阵的主对角线元素之积即为det A 。 由于A 的元素零元素较多，为节省储存量，将A 的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A 中的元素。 2.全部源程序 #include 第二章 1. 题目：运用MATLAB编程实现牛顿迭代 2. 实验操作 1、打开MATLAB程序软件。 2、在MATLAB中编辑如下的M程序。 function [p1,err,k,y]=newton(f,df,p0,delta,max) %f 是要求根的方程(f(x)=0); %df 是f(x)的导数； %p0是所给初值,位于x*附近； %delta是给定允许误差； %max是迭代的最大次数； %p1是newton法求得的方程的近似解； %err是p0的误差估计； %k是迭代次数； p0 for k=1:max p1=p0-feval('f',p0)/feval('df',p0); err=abs(p1-p0); p0=p1; k p1 err y=feval('f',p1) if (err p0 = 1.2000 k =1 p1=1.1030 err=0.0970 y=0.0329 k= 2 p1=1.0524 err=0.0507 y=0.0084 k =3 p1=1.0264 err=0.0260 y=0.0021 k =4 p1=1.0133 err=0.0131 y=5.2963e-004 k =5 p1=1.0066 err=0.0066 y=1.3270e-004 k =6 p1=1.0033 err=0.0033 y=3.3211e-005 k =7 p1=1.0017 err=0.0017 y=8.3074e-006 k =8 p1=1.0008 err=8.3157e-004 y = 2.0774e-006 k =9 p1=1.0004 err=4.1596e-004 y =5.1943e-007 k=10 p1=1.0002 err=2.0802e-004 y= 1.2987e-007 k=11 p1=1.0001 err=1.0402e-004 y =3.2468e-008 k=12 p1=1.0001 err=5.2014e-005 y=8.1170e-009 k=13 p1=1.0000 err=2.6008e-005 y= 2.0293e-009 k=14 p1=1.0000 err=1.3004e-005 y=5.0732e-010 k=15 p1 =1.0000 err=6.5020e-006 y=1.2683e-010 k=16 p1 =1.0000 err=3.2510e-006 y=3.1708e-011 k=17 p1 =1.0000 err=1.6255e-006 y =7.9272e-012 k=18 p1 =1.0000 err =8.1279e-007 y= 1.9820e-012 ans = 1.0000 结果说明：经过18次迭代得到精确解为1，误差为8.1279e-007。 研究生《数值分析》课程作业（二） 姓名： 学号： 专业： 1、据如下函数值表，建立二次的Lagrange 插值多项式及Newton 插值多项式。 20012222()()()()()()() (1)(2)(0)(2)(-0)(1)59 3143 (01)(02)(10)(12(20)(21)22 L x f x l x f x l x f x l x x x x x x x x x =++-----=? +?+?=-+------解： 二次 l agr ange插值 ） Newton 插值多项式： 200100120122()()[,](-)[,,](-)(-) 5559 32(0)(0)(1)32()3 2222 N x f x f x x x x f x x x x x x x x x x x x x x x =++=-?-+--=-+-=-+ ()y f x =2、已知单调连续函数在如下采样点处的函数值 *()0[2,4],f x x =求方程在内根的近似值使误差尽可能小。 解：1 ()()y f x x f y -==解： 对的反函数进行二次插值 1110201122012010210122021(0)(0)(0)(0)(0)(0) (0)() ()() ()()()()()() (0 2.25)(05)(03)(05)(03)(0 2.25) 2 3.54( 3 2.25)(35)(2.253)(2.255)(53)(5 2.25) y y y y y y L f y f y f y y y y y y y y y y y y y ---------=++--------+-+-=? +?+? ----+-+- 2.945 ≈()(1)01(1)1()[,]()(,),()[,],() ()()()() (1)! ,n n n n n n n n f x a b f x a b a x x x b L x x a b f R x f x L x x n a b x ξωξ+++≤<<<≤∈=-=+∈ 3、证明：设在上连续，在内存在，节点是满足拉格朗日插值条件的多项式，则对任何插值余项 这里（）且依赖于。 0110101(0,1,,)()()0()()()()()()()()[,]()()()()()()() (),,,(k n n k n n n n n n x k n R x R x R x K x x x x x x x K x x K x x x a b t f t L t K x t x t x t x t x x x x t ωφφφ+===---==----- 证由条件知节点是的零点，即。于是其中是与有关的待定函数。 现把看成上的固定点，作函数 根据插值条件和余项定义，知在点及处均为零。故明：1111)[,]2()[,]1()()[,]()(,)(,),()()(1)!()0 ()()(,),(1)! n n n n a b n t a b n t t a b n t a b a b f n K x f K x a b x n φφφφξφξξξξ++++'+'''+∈=-+==∈+（） （） （）（）在上有个零点，根据罗尔定理，在内至少有个零点。对再应用罗尔定理，可知在内至少 有个零点。依次类推，在上至少有一个零点，记为 使 于是 ， 且依赖于于是得到插值余项。 证毕。 44、试用数据表建立不超过次的埃尔米特插值多项式。 解：（用重节点的均差表建立埃尔米特多项式） 数值分析上机第四次作业 实验项目：共轭梯度法求解对称正定的线性方程组 实验内容：用共轭梯度法求解下面方程组 (1) 123421003131020141100155x x x x -?????? ? ? ?--- ? ? ?= ? ? ?-- ? ? ?-???? ?? 迭代20次或满足()(1) 1110k k x x --∞-<时停止计算。 (2) Ax b =,A 是1000阶的Hilbert 矩阵或如下的三对角矩阵， A[i,i]=4，A[i,i-1]=A[i-1,i]=-1，i=2,3,..,n b[1]=3, b[n]=3, b[i]=2，i=2,3,…,n-1 迭代10000次或满足()()710k k r b Ax -=-≤时停止计算。 （3）*考虑模型问题，方程为 222222(),(,)(0,1)(0,1)(0,)1,(1,),01(,0)1,(,1),01 xy y x u u x y e x y D x y u y u y e y u x u x e x ??+=+∈=???==≤≤==≤≤ 用正方形网格离散化，若取1/,10h N N ==，得到100n =的线性方程组，并用共轭梯度法（CG 法）求解，并对解作图。 实验要求：迭代初值可以取01(,1,...,)ij u i j N ==，计算到32||||10k r -≤停止．本 题有精确解(,)xy u x y e =，这里k u 表示以k ij u 为分量的向量， u 表示在相应点(,)i j 上取值作为分量的向量． 实验一： （1） 编制函数子程序CGmethod 。 function [x,k]=CGmethod(A,b) n=length(A);x=zeros(n,1);r=b-A*x;rho=r'*r; k=0; while rho>10^(-12) & k<20 k=k+1; if k==1 p=r; else beta=rho/rho1; p=r+beta*p; end 大作业 三 1. 给定初值0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用 程序. 解：Matlab 程序如下： 函数m 文件：fu.m function Fu=fu(x) Fu=x^3/3-x; end 函数m 文件：dfu.m function Fu=dfu(x) Fu=x^2-1; end 用Newton 法求根的通用程序Newton.m clear; x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0) while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0) 北京航空航天大学 数值分析大作业八 学院名称自动化 专业方向控制工程 学号 学生姓名许阳 教师孙玉泉 日期2014 年11月26 日 一．题目 关于x , y , t , u , v , w 的方程组(A.3) ???? ?? ?=-+++=-+++=-+++=-+++79 .0sin 5.074.3cos 5.007.1cos sin 5.067.2cos 5.0y w v u t x w v u t y w v u t x w v u t (A.3) 以及关于z , t , u 的二维数表（见表A-1）确定了一个二元函数z =f (x , y )。 表A-1 二维数表 t z u 0 0.4 0.8 1.2 1.6 2 0 -0.5 -0.34 0.14 0.94 2.06 3.5 0.2 -0.42 -0.5 -0.26 0.3 1.18 2.38 0.4 -0.18 -0.5 -0.5 -0.18 0.46 1.42 0.6 0.22 -0.34 -0.58 -0.5 -0.1 0.62 0.8 0.78 -0.02 -0.5 -0.66 -0.5 -0.02 1.0 1.5 0.46 -0.26 -0.66 -0.74 -0.5 1. 试用数值方法求出f (x , y ) 在区域}5.15.0,8.00|), {≤≤≤≤=y x y x D （上的近似表达式 ∑∑===k i k j s r rs y x c y x p 00 ),( 要求p (x , y )以最小的k 值达到以下的精度 ∑∑==-≤-=10020 7210)],(),([i j i i i i y x p y x f σ 其中j y i x i i 05.05.0,08.0+==。 2. 计算),(),,(* ***j i j i y x p y x f (i =1,2,…,8 ; j =1,2,…,5) 的值，以观察p (x , y ) 逼 近f (x , y )的效果，其中j y i x j i 2.05.0,1.0**+==。数值分析上机作业
数值分析上机作业
数值分析大作业-三、四、五、六、七
东南大学数值分析上机作业汇总
数值分析上机题课后作业全部-东南大学
数值分析大作业三 四 五 六 七
Matlab作业3(数值分析)答案
数值分析作业思考题汇总
数值计算方法大作业
(完整版)数值计算方法上机实习题答案
数值分析大作业
高等数值分析上机作业
数值分析上机作业1-1
北航数值分析大作业第一题幂法与反幂法
数值分析作业
研究生《数值分析》课程作业(二) (含答案)
数值分析上机第四次作业
数值分析大作业 三、四、五、六、七
北航数值分析报告大作业第八题