数值计算方法上机实验报告
- 格式:docx
- 大小:327.52 KB
- 文档页数:33
数值代数上机实验报告试验项目名称:平方根法与改进平方根法实验内容:先用你熟悉的计算机语言将平方根法和改进平方根法编写成通用的子程序,然后用你编写的程序求解对称正定方程组Ax=b,其中,A=[101 10 1…1 10 11 10]100*100b随机生成,比较计算结果,评论方法优劣。
实验要求:平方根法与改进的平方根的解法步骤;存储单元,变量名称说明;系数矩阵与右端项的生成;结果分析。
实验报告姓名:罗胜利班级:信息与计算科学0802 学号:u200810087 实验一、平方根法与改进平方根法先用你所熟悉的计算机语言将平方根法和改进的平方根法编成通用的子程序,然后用你编写的程序求解对称正定方程组AX=b,其中系数矩阵为40阶Hilbert矩阵,即系数矩阵A的第i行第j列元素为=,向量b的第i个分量为=.平方根法函数程序如下:function [x,b]=pingfanggenfa(A,b)n=size(A);n=n(1);x=A^-1*b; %矩阵求解disp('Matlab自带解即为x');for k=1:nA(k,k)=sqrt(A(k,k));A(k+1:n,k)=A(k+1:n,k)/A(k,k);for j=k+1:n;A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endend %Cholesky分解for j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n); %前代法A=A';for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1); %回代法disp('平方根法的解即为b');endfunction [x]=ave(A,b,n) %用改进平方根法求解Ax=b L=zeros(n,n); %L为n*n矩阵D=diag(n,0); %D为n*n的主对角矩阵S=L*D;for i=1:n %L的主对角元素均为1L(i,i)=1;for i=1:nfor j=1:n %验证A是否为对称正定矩阵if (eig(A)<=0)|(A(i,j)~=A(j,i)) %A的特征值小于0或A非对称时,输出wrong disp('wrong');break;endendendD(1,1)=A(1,1); %将A分解使得A=LDL Tfor i=2:nfor j=1:i-1S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)');L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1);endD(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');endy=zeros(n,1); % x,y为n*1阶矩阵x=zeros(n,1);for i=1:ny(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)))/D(i,i); %通过LDy=b解得y的值endfor i=n:-1:1x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n)); %通过L T x=y解得x的值改进平方根法函数程序如下:function b=gaijinpinfanggenfa(A,b)n=size(A);n=n(1);v=zeros(n,1);for j=1:nfor i=1:j-1v(i)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j);end %LDL'分解B=diag(A);D=zeros(n);for i=1:nD(i,i)=B(i);A(i,i)=1;EndA=tril(A); %得到L和Dfor j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n); %前代法A=D*(A');for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1); %回代法disp('改进平方根法解得的解即为b');end调用函数解题:clear;clc;n=input('请输入矩阵维数:');b=zeros(n,1);A=zeros(n);for i=1:nfor j=1:nA(i,j)=1/(i+j-1);b(i)=b(i)+1/(i+j-1);endend %生成hilbert矩阵[x,b]=pingfanggenfa(A,b) b=gaijinpinfanggenfa(A,b)运行结果:请输入矩阵维数:40Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.570692e-020. > In pingfanggenfa at 4In qiujie at 10Matlab自带解即为x平方根法的解即为bx =1.60358.96850.85621.01950.9375-50.2500-3.0000-16.000024.0000-49.5000-30.000039.000022.0000-64.0000-12.00002.000010.2500-10.5000-1.0000-10.875083.000046.0000-98.0000-69.000068.000021.0000-50.7188-8.7500-8.0000 112.0000 6.0000 -68.7500 22.000044.0000 -28.0000 8.0000 -44.000012.0000b =1.0e+007 *0.0000-0.00000.0001-0.0004-0.00140.0424-0.29801.1419-2.73354.2539-4.30182.7733-1.19890.5406-0.36880.32850.4621-0.25130.05650.0000-0.00510.0071-0.0027-0.0031-0.00190.00090.0002-0.0002-0.00060.00040.0001-0.00020.00010.0000-0.00000.0000-0.0000-0.0000改进平方根法解得的解即为bb =1.0e+024 *0.0000-0.00000.0001-0.0012-0.0954 0.4208 -1.2101 2.0624 -1.0394 -3.3343 6.2567 -0.2463 -7.45942.80303.6990 0.7277 -1.7484 -0.4854 -3.6010 0.2532 5.1862 1.4410 0.8738 -4.5654 1.0422 4.0920 -2.7764 -2.2148 -0.8953 0.3665 4.8967 1.0416 0.1281-1.1902-2.83348.4610-3.6008实验二、利用QR分解解线性方程组:利用QR分解解线性方程组Ax=b,其中A=[16 4 8 4;4 10 8 4;8 8 12 10;4 4 10 12];b=[32 26 38 30];求解程序如下:定义house函数:function [v,B]=house(x)n=length(x);y=norm(x,inf);x=x/y;Q=x(2:n)'*x(2:n);v(1)=1;v(2:n)=x(2:n);if n==1B=0;elsea=sqrt(x(1)^2+Q);if x(1)<=0v(1)=x(1)-a;elsev(1)=-Q/(x(1)+a);endB=2*v(1)^2/(Q+v(1)^2);endend进行QR分解:clear;clc;A=[16 4 8 4;4 10 8 4;8 8 12 10;4 4 10 12]; b=[32 26 38 30];b=b';x=size(A);m=x(1);n=x(2);d=zeros(n,1);for j=1:n[v,B]=house(A(j:m,j));A(j:m,j:n)=(eye(m-j+1)-B*(v')*v)*A(j:m,j:n); d(j)=B;if j<m< p="">A(j+1:m,j)=v(2:m-j+1);endend %QR分解R=triu(A); %得到R D=A;I=eye(m,n);Q=I;for i=1:nD(i,i)=1;endH=tril(D);M=H';for i=1:nN=I-d(i)*H(1:m,i)*M(i,1:m);Q=Q*N;end %得到Qb=(Q')*b; %Q是正交阵for j=n:-1:2b(j)=b(j)/R(j,j);b(1:j-1)=b(1:j-1)-b(j)*R(1:j-1,j);endb(1)=b(1)/R(1,1); %回带法运行结果如下:R =18.7617 9.8072 15.7769 11.08640 9.9909 9.3358 7.53410 0 5.9945 9.80130 0 0 -0.5126Q =0.8528 -0.4368 -0.2297 -0.17090.2132 0.7916 -0.4594 -0.34170.4264 0.3822 0.2844 0.76890.2132 0.1911 0.8095 -0.5126b=1.000000000000001.000000000000010.9999999999999881.00000000000001实验三、Newton下山法解非线性方程组:3x-cos(yz)-=0,-81+sinz+1.06=0,exp(-xy)+20z+=0;要求满足数值解=满足或.定义所求方程组的函数:Newtonfun.mfunction F = Newtonfun(X)F(1,1)=3*X(1)-cos(X(2)*X(3))-1/2;F(2,1)=X(1)^2-81*(X(2)+0.1)^2+sin(X(3))+1.06;F(3,1)=exp(-X(1)*X(2))+20*X(3)+(10*pi-3)/3;End向量求导:Xiangliangqiudao.mfunction J=xiangliangqiudao()syms x y zX=[x,y,z];F=[3*X(1)-cos(X(2)*X(3))-1/2;X(1)^2-81*(X(2)+0.1)^2+sin(X(3))+1.06;exp(-X(1)*X(2))+20*X(3)+(10*pi-3)/3];J=jacobian(F,[x y z]);End代值函数:Jacobi.mfunction F=Jacobi(x)F=[ 3,x(3)*sin(x(2)*x(3)), x(2)*sin(x(2)*x(3));2*x(1), -162*x(2)-81/5,cos(x(3));-x(2)/exp(x(1)*x(2)),-x(1)/exp(x(1)*x(2)),20];End方程组求解:format long; %数据表示为双精度型X1=[0,0,0]';eps=10^(-8);k=1;i=1;X2=X1-Jacobi(X1)^(-1)*Newtonfun(X1);while (norm(subs(X2-X1,pi,3.1415926),2)>=eps)&&(norm(Newtonfun(X1),2)>=eps) if norm(Newtonfun(X2),2)<="" p="">X1=X2;B=inv(Jacobi(X2));C=Newtonfun(X2);X2=X2-B*C;i=i+1;elsev=1/(2^k); %引入下山因子X1=X2;B=inv(Jacobi(X2));C=Newtonfun(X2);X2=X2-v*B*C;k=k+1;endendj=i+k-1 %迭代次数X=X2 %输出结果运行结果如下:j =5X =0.500000000000000 -0.000000000000000 -0.523598775598299</m<>。
本科实验报告课程名称:计算机数值方法实验项目:方程求根、线性方程组的直接解法、线性方程组的迭代解法、代数插值实验地点:专业班级:学生姓名:指导教师:实验一方程求根}五、实验结果与分析二分法实验结果迭代法实验结果结果分析:本题目求根区间为[1,2],精度满足|x*-x n|<0.5×10-5,故二分法用公式|x*-x n|<(b-a)/ 2n,可求得二分次数并输出每次结果。
对迭代法首先要求建立迭代格式。
迭代格式经计算已输入程序之中,故直接给初值便可利用迭代法求出精度下的解。
六、讨论、心得每次的实验都是对已学过的理论知识的一种实战。
通过本次实验,我将二分法与迭代法的思路清晰化并且将其变成计算机设计语言编写出来,运用到了实际解决问题上感觉很好。
我自认为本次跟其他同学比较的优点在于我在二分法实现的时候首先利用换底公式将需要的二分次输输出,如此便很清晰明了的知道接下来每一步的意思。
迭代法给我的感觉便是高度的便捷简化,仅用几行代码便可以同样解决问题。
相比较二分法来说,我更喜欢迭代的思路。
实验二线性方程组的直接解法for(k=n-2;k>=0;k--){sum=0;for(j=k+1;j<n;j++)sum=sum+a[k][j]*x[j];x[k]=(b[k]-sum)/a[k][k];}for(i=0;i<n;i++)printf("x[%d]=%f ",i,x[i]); printf("\n"); //输出解向量x}五、实验结果与分析结果结果分析:如上图所示,输入线性方程组元数n=3,则会要求输入3*3的系数矩阵A与向量b构成的增广矩阵。
根据算法需要将系数矩阵A消元成上三角矩阵。
随后根据矩阵乘法公式变形做对应的回代。
六、讨论、心得本次实验在编写时候感觉还好,感觉将思路变成了程序设计语言,得以实现题目的要求。
但是在运行以及结果分析的时候,感觉到了本实验的一些不足之处:就是我的实验虽然可以实现不同的元数的线性方程组求解,但是缺少了分析初始条件——主元素不能为零。
1.拉格朗日插值多项式,用于离散数据的拟合#include <stdio.h>#include <conio.h>#include <alloc.h>float lagrange(float *x,float *y,float xx,int n) /*拉格朗日插值算法*/{ int i,j;float *a,yy=0.0; /*a作为临时变量,记录拉格朗日插值多项式*/a=(float *)malloc(n*sizeof(float));for(i=0;i<=n-1;i++){ a[i]=y[i];for(j=0;j<=n-1;j++)if(j!=i) a[i]*=(xx-x[j])/(x[i]-x[j]);yy+=a[i];}free(a);return yy;}main(){ int i,n;float x[20],y[20],xx,yy;printf("Input n:");scanf("%d",&n);if(n>=20) {printf("Error!The value of n must in (0,20)."); getch();return 1;} if(n<=0) {printf("Error! The value of n must in (0,20)."); getch(); return 1;} for(i=0;i<=n-1;i++){ printf("x[%d]:",i);scanf("%f",&x[i]);}printf("\n");for(i=0;i<=n-1;i++){ printf("y[%d]:",i);scanf("%f",&y[i]);}printf("\n");printf("Input xx:");scanf("%f",&xx);yy=lagrange(x,y,xx,n);printf("x=%f,y=%f\n",xx,yy);getch();}2.牛顿插值多项式,用于离散数据的拟合#include <stdio.h>#include <conio.h>#include <alloc.h>void difference(float *x,float *y,int n){ float *f;int k,i;f=(float *)malloc(n*sizeof(float));for(k=1;k<=n;k++){ f[0]=y[k];for(i=0;i<k;i++)f[i+1]=(f[i]-y[i])/(x[k]-x[i]);y[k]=f[k];}return;}main(){ int i,n;float x[20],y[20],xx,yy;printf("Input n:");scanf("%d",&n);if(n>=20) {printf("Error! The value of n must in (0,20)."); getch(); return 1;} if(n<=0) {printf("Error! The value of n must in (0,20).");getch(); return 1;} for(i=0;i<=n-1;i++){ printf("x[%d]:",i);scanf("%f",&x[i]);}printf("\n");for(i=0;i<=n-1;i++){ printf("y[%d]:",i);scanf("%f",&y[i]);}printf("\n");difference(x,(float *)y,n);printf("Input xx:");scanf("%f",&xx);yy=y[20];for(i=n-1;i>=0;i--) yy=yy*(xx-x[i])+y[i];printf("NewtonInter(%f)=%f",xx,yy);getch();}3.高斯列主元消去法,求解其次线性方程组第一种#include<stdio.h>#include <math.h>#define N 20int main(){ int n,i,j,k;int mi,tmp,mx;float a[N][N],b[N],x[N];printf("\nInput n:");scanf("%d",&n);if(n>N){ printf("The input n should in(0,N)!\n");getch();return 1;}if(n<=0){ printf("The input n should in(0,N)!\n");getch();return 1;}printf("Now input a(i,j),i,j=0...%d:\n",n-1); for(i=0;i<n;i++){ for(j=0;j<n;j++)scanf("%f",&a[i][j]);}printf("Now input b(i),i,j=0...%d:\n",n-1);for(i=0;i<n;i++)scanf("%f",&b[i]);for(i=0;i<n-2;i++){ for(j=i+1,mi=i,mx=fabs(a[i][j]);j<n-1;j++) if(fabs(a[j][i])>mx){ mi=j;mx=fabs(a[j][i]);}if(i<mi){ tmp=b[i];b[i]=b[mi];b[mi]=tmp;for(j=i;j<n;j++){ tmp=a[i][j];a[i][j]=a[mi][j];a[mi][j]=tmp;}}for(j=i+1;j<n;j++){ tmp=-a[j][i]/a[i][i];b[j]+=b[i]*tmp;for(k=i;k<n;k++)a[j][k]+=a[i][k]*tmp;}}x[n-1]=b[n-1]/a[n-1][n-1];for(i=n-2;i>=0;i--){ x[i]=b[i];for(j=i+1;j<n;j++)x[i]-=a[i][j]*x[j];x[i]/=a[i][i];}for(i=0;i<n;i++)printf("Answer:\n x[%d]=%f\n",i,x[i]); getch();return 0;}第二种#include<math.h>#include<stdio.h>#define NUMBER 20#define Esc 0x1b#define Enter 0x0dfloat A[NUMBER][NUMBER+1] ,ark;int flag,n;exchange(int r,int k);float max(int k);message();main(){float x[NUMBER];int r,k,i,j;char celect;clrscr();printf("\n\nUse Gauss.");printf("\n\n1.Jie please press Enter."); printf("\n\n2.Exit press Esc.");celect=getch();if(celect==Esc)exit(0);printf("\n\n input n=");scanf("%d",&n);printf(" \n\nInput matrix A and B:"); for(i=1;i<=n;i++){printf("\n\nInput a%d1--a%d%d and b%d:",i,i,n,i);for(j=1;j<=n+1;j++) scanf("%f",&A[i][j]);}for(k=1;k<=n-1;k++){ark=max(k);if(ark==0){printf("\n\nIt's wrong!");message();}else if(flag!=k)exchange(flag,k);for(i=k+1;i<=n;i++)for(j=k+1;j<=n+1;j++)A[i][j]=A[i][j]-A[k][j]*A[i][k]/A[k][k];}x[n]=A[n][n+1]/A[n][n];for( k=n-1;k>=1;k--){float me=0;for(j=k+1;j<=n;j++){me=me+A[k][j]*x[j];}x[k]=(A[k][n+1]-me)/A[k][k];}for(i=1;i<=n;i++){printf(" \n\nx%d=%f",i,x[i]);}message();}exchange(int r,int k){int i;for(i=1;i<=n+1;i++)A[0][i]=A[r][i];for(i=1;i<=n+1;i++)A[r][i]=A[k][i];for(i=1;i<=n+1;i++)A[k][i]=A[0][i];}float max(int k){int i;float temp=0;for(i=k;i<=n;i++)if(fabs(A[i][k])>temp){temp=fabs(A[i][k]);flag=i;}return temp;}message(){printf("\n\n Go on Enter ,Exit press Esc!");switch(getch()){case Enter: main();case Esc: exit(0);default:{printf("\n\nInput error!");message();} }}4.龙贝格求积公式,求解定积分#include<stdio.h>#include<math.h>#define f(x) (sin(x)/x)#define N 20#define MAX 20#define a 2#define b 4#define e 0.00001float LBG(float p,float q,int n){ int i;float sum=0,h=(q-p)/n;for (i=1;i<n;i++)sum+=f(p+i*h);sum+=(f(p)+f(q))/2;return(h*sum);}void main(){ int i;int n=N,m=0;float T[MAX+1][2];T[0][1]=LBG(a,b,n);n*=2;for(m=1;m<MAX;m++){ for(i=0;i<m;i++)T[i][0]=T[i][1];T[0][1]=LBG(a,b,n);n*=2;for(i=1;i<=m;i++)T[i][1]=T[i-1][1]+(T[i-1][1]-T[i-1][0])/(pow(2,2*m)-1);if((T[m-1][1]<T[m][1]+e)&&(T[m-1][1]>T[m][1]-e)){ printf("Answer=%f\n",T[m][1]); getch();return ;}}}5.牛顿迭代公式,求方程的近似解#include<stdio.h>#include<math.h>#include<conio.h>#define N 100#define PS 1e-5#define TA 1e-5float Newton(float (*f)(float),float(*f1)(float),float x0 ) { float x1,d=0;int k=0;do{ x1= x0-f(x0)/f1(x0);if((k++>N)||(fabs(f1(x1))<PS)){ printf("\nFailed!");getch();exit();}d=(fabs(x1)<1?x1-x0:(x1-x0)/x1);x0=x1;printf("x(%d)=%f\n",k,x0);}while((fabs(d))>PS&&fabs(f(x1))>TA) ;return x1;}float f(float x){ return x*x*x+x*x-3*x-3; }float f1(float x){ return 3.0*x*x+2*x-3; }void main(){ float f(float);float f1(float);float x0,y0;printf("Input x0: ");scanf("%f",&x0);printf("x(0)=%f\n",x0);y0=Newton(f,f1,x0);printf("\nThe root is x=%f\n",y0); getch();}。
数值积分上机实验报告141110038 桂贤进题一:数学上已经证明了f1 4--- dx =nJo 1+0成立,所以可以通过数值积分来求71的近似值。
1.分别使用复合梯形、复合Simpson求积公式计算11的近似值。
选择不同的h,对每种求枳公式,试将误差刻画为h的函数,并比较两方法的精度。
是否存在某个值,当低于这个值之后,再继续减小h的值,计算精度不再有所改进,为什么?2.实现Romberg求积方法,并重复上面的计算;3•实现自适应积分方法,并审;复上面的计算。
解:1.1公式分析:(a)对于复合梯形公式T"=纟[f (a) +2£f(a + 汎)]“=字⑴i=lE n(f)=-嗒f⑺⑺= ①严)(f),a v f V 方(b)对于复合Simpson公式斤m m—1SnG) = £ [/(a) + f(b) + 4》f(a + ⑵ 一1)/1) + 2》f(a + 2ih)](3)1=1 1=1. b-a b-a11 = —= --------2m n离散误差为:EQ—嘗八呢一?^炉肿vg. (4)1.2实现算法结果:分别利用梯形公式和Simpson公式计算结果如下:(下表中E丄(f) = \n-T(f儿E2(T) = |兀此处兀为MATLAB中的数,可以认为具冇足够大的精度)61/6 3.136963066471264 0.00463 3.141591780936043 8.7265e-07 8 1/8 3.138988494491089 0.00260 3.141592502458707 1.5113e-07 10 1/10 3.139925988907159 0.00167 3.141592613939215 3.9651e-08 12 1/12 3.140435246846851 0.00116 3.141592********* 1.3284e ・08 20 1/20 3.141175986954129 4.1667e-04 3.141592652969785 6.2001e-10 30 1/30 3.141407468407330 1.8519e-04 3.141592653535359 5.4434e-ll 40 1/40 3.141488486923612 1.0417e-04 3.141592653580105 9.6878e-12 50 1/50 3.141525986923254 6.6667e-05 3.141592653587253 2.5402e-12 1001/100 3.141575986923129 1.6667e-05 3.141592653589753 3.9968e-142001/2003.1415884869231304.1667e-063.141592653589793从上农中可以看出:复合Simpson 公式比复合梯形公式稱度岛,误差收敛的速度快不少。
数值计算方法实验报告一、实验目的本实验旨在探究数值计算方法在实际应用中的效果和应用范围,验证其可靠性和准确性。
二、实验原理数值计算方法是一种计算机辅助的数值分析方法,通过数学模型和计算机算法,对复杂的数学问题进行数值近似解。
其基本思想是将数学模型转化为离散的数学问题,利用计算机进行计算,得出逼近精度可控的数值解。
在本次实验中,我们主要涉及以下数值计算方法:1.插值法:通过已知数据点之间的曲线进行近似曲线的构造,从而推导出相应函数在某些点上的近似值。
2.数值积分法:将区间内的函数用多项式逼近,计算逼近多项式在该区间内的积分值,从而得出原函数在该区间内的近似积分值。
3.常微分方程数值解法:通过数值计算的方式对常微分方程进行求解,得出方程的近似解。
三、实验过程1.插值法实验(1)使用拉格朗日插值法对给定的一组数据进行插值,并画出插值曲线。
(2)使用牛顿插值法对同样的一组数据进行插值,并比较两种插值方法的效果。
2.数值积分法实验(1)使用复合梯形公式求解给定的一元函数在某个区间上的近似积分值,并与准确值进行比较。
(2)使用复合辛普森公式求解同样的一元函数在同一区间上的近似积分值,并比较两种方法的精度。
3.常微分方程数值解法实验(1)使用欧拉法求解给定的一阶常微分方程,并与准确解进行比较。
(2)使用四阶龙格-库塔法求解同样的一阶常微分方程,并比较两种方法的精度。
四、实验结果1.插值法实验结果(1)拉格朗日插值法的插值曲线如下:(2)牛顿插值法的插值曲线如下:从图中可以看出,两种插值方法的效果相似,但牛顿插值法的插值曲线更加平滑。
2.数值积分法实验结果(1)复合梯形公式的近似积分值为0.6932,准确值为0.6931。
相对误差为0.022%。
(2)复合辛普森公式的近似积分值为0.6931,准确值为0.6931。
相对误差为0.003%。
由上述结果可以看出,复合辛普森公式的精度更高。
3.常微分方程数值解法实验结果(1)欧拉法的数值解如下:(2)四阶龙格-库塔法的数值解如下:从图中可以看出,四阶龙格-库塔法的数值解更加接近于准确解。
数值分析上机实践报告一、实验目的本实验的目的是通过编写数值分析程序,掌握解决数学问题的数值计算方法,并通过实际应用来检验其有效性和准确性。
具体包括以下几个方面的内容:1.掌握二分法和牛顿迭代法的基本原理和实现方法;2.熟悉利用矩阵的LU分解和追赶法解线性方程组的过程;3.通过具体的实例应用,比较不同方法的计算效果和精度。
二、实验内容本实验分为三个部分,每个部分包括一个具体的数学问题和相应的数值计算方法。
1.问题一:求方程f(x)=x^3-5x^2+10x-80=0的近似解。
在问题一中,我们通过二分法和牛顿迭代法来求解方程的近似解,并比较两种方法的精度和收敛速度。
2.问题二:用LU分解解线性方程组。
问题二中,我们通过矩阵的LU分解方法解线性方程组Ax=b,然后和直接用追赶法解线性方程组进行对比,验证LU分解的有效性和准确性。
三、实验结果及分析1.问题一的结果分析:通过二分法和牛顿迭代法求解方程f(x)=x^3-5x^2+10x-80=0的近似解,得到的结果如下:从结果来看,两种方法得到的近似解均与真实解x≈5非常接近。
但是,通过比较可以发现,牛顿迭代法的计算速度比二分法更快,迭代的次数更少。
因此,在需要高精度近似解的情况下,牛顿迭代法是一个更好的选择。
2.问题二的结果分析:通过LU分解和追赶法解线性方程组Ax=b,得到的结果如下:-用LU分解解线性方程组得到的结果为x1≈1.0,x2≈2.0,x3≈3.0;-用追赶法解线性方程组得到的结果为x1≈1.0,x2≈2.0,x3≈3.0。
从结果来看,两种方法得到的结果完全一致,而且与真实解非常接近。
这表明LU分解方法和追赶法均可以有效地解决线性方程组问题。
但是,在实际应用中,当方程组规模较大时,LU分解方法的计算复杂度较高,因此追赶法更加适用。
四、实验总结通过本实验,我掌握了二分法和牛顿迭代法以及LU分解和追赶法的基本原理和实现方法。
通过具体的数学问题实例应用,我比较了不同方法的计算效果和精度,得出以下结论:1.在求解函数的近似解时,牛顿迭代法相对于二分法具有更快的收敛速度和更高的计算精度;2.在解决线性方程组问题时,LU分解方法在计算准确性方面与追赶法相当,但在处理较大规模的问题时,计算复杂度较高,追赶法更适合。
. 《计算方法》上机实验报告班级:XXXXXX小组成员:XXXXXXXXXXXXXXXXXXXXXXXXXXXX任课教师:XXX二〇一八年五月二十五日前言通过进行多次的上机实验,我们结合课本上的内容以及老师对我们的指导,能够较为熟练地掌握Newton 迭代法、Jacobi 迭代法、Gauss-Seidel 迭代法、Newton 插值法、Lagrange 插值法和Gauss 求积公式等六种算法的原理和使用方法,并参考课本例题进行了MATLAB 程序的编写。
以下为本次上机实验报告,按照实验内容共分为六部分。
实验一:一、实验名称及题目: Newton 迭代法例2.7(P38):应用Newton 迭代法求在附近的数值解,并使其满足.二、解题思路:设'x 是0)(=x f 的根,选取0x 作为'x 初始近似值,过点())(,00x f x 做曲线)(x f y =的切线L ,L 的方程为))((')(000x x x f x f y -+=,求出L 与x 轴交点的横坐标)(')(0001x f x f x x -=,称1x 为'x 的一次近似值,过点))(,(11x f x 做曲线)(x f y =的切线,求该切线与x 轴的横坐标)(')(1112x f x f x x -=称2x 为'x 的二次近似值,重复以上过程,得'x 的近似值序列{}n x ,把)(')(1n n n n x f x f x x -=+称为'x 的1+n 次近似值,这种求解方法就是牛顿迭代法。
三、Matlab 程序代码:function newton_iteration(x0,tol) syms z %定义自变量 format long %定义精度 f=z*z*z-z-1; f1=diff(f);%求导 y=subs(f,z,x0);y1=subs(f1,z,x0);%向函数中代值 x1=x0-y/y1; k=1;while abs(x1-x0)>=tol x0=x1;y=subs(f,z,x0); y1=subs(f1,z,x0); x1=x0-y/y1;k=k+1; endx=double(x1) K四、运行结果:实验二:一、实验名称及题目:Jacobi 迭代法例3.7(P74):试利用Jacobi 迭代公式求解方程组要求数值解为方程组的精确解. 二、解题思路:首先将方程组中的系数矩阵A 分解成三部分,即:U D L A ++=,D 为对角阵,L 为下三角矩阵,U 为上三角矩阵。
学院:自动化学院班级:自动化085姓名:学号:2011年3月一、实验的性质、目的和任务本实验是与本专业基础课《数值计算方法》相配套的,旨在巩固专业课内容和学生编程的能力。
通过实验加强对数值方法的理解和掌握,编制出适用的程序。
同时,在理论教学的基础上,注意方法处理的技巧及其与计算机的结合,;其次要通过例子,学习使用各种数值方法解决实际计算问题。
要求学生应用高级计算机语言Matlab编程完成实验。
二、实验基本要求要求熟悉高级计算机语言Matlab,以及相关上机操作说明;上机时要遵守实验室的规章制度,爱护实验设备;记录调试过程及结果。
三、实验原理应用高级计算机语言实现数值计算方法课程中的各种算法。
四、设备及器材配置主机:微机操作系统:WINDOWS 98以上软件:高级计算机语言Matlab五、考核与报告每个实验完成后交一份实验报告。
本实验作为平时成绩的一部分占学期期末总成绩的20%。
六、适用对象自动化专业七、主要参考书1.王能超编,《数值分析简明教程》,高等教育出版社,2003年,第2版2.封建湖编,《数值分析原理》,科学出版社,2001年,第1版3.冯有前编,《数值分析》,清华大学出版社,2005年,第1版4.周璐等译, John H. Mathews等编,《数值方法(MATLAB版)》,电子工业出版社,2007年,第二版实验一 采用拉格朗日方法计算插值一、 实验目的:1. 掌握多项式插值的概念、存在唯一性;2. 能够熟练地应用拉格朗日方法计算插值,并完成插值程序的设计和调试。
二、 实验内容:构造拉格朗日插值多项式()p x 逼近3()f x x =,要求:(1) 取节点01x =-,11x =求线性插值多项式1()p x ;(2) 取节点01x =-,10x =,21x =求抛物插值多项式2()p x ;(3) 取节点01x =-,10x =,21x =,32x =求三次插值多项式3()p x ;(4) 分别求1(1.3)p 、2(1.3)p 、3(1.3)p 的值,并与精确值相比较。
华 北 电 力 大 学 科 技 学 院
实 验 报 告
课程名称:数值计算方法上机实验报告 姓名:牛玺童 班级:电气 11k6 学号:111904010415 1
k k * k *
1.题目 实验一 造倒数表 造倒数表,并例求 18 的倒数。(精度为 0.0005) 2.算法原理
2.1 牛顿迭代法 牛顿迭代法是通过非线性方程线性化得到迭代序列的一种方法。
对于非线性方程 f (x) 0 ,若已知根 x* 的一个近似值 x ,将 f (x) 在 x 处展 成一阶泰勒公式后忽略高次项可得: f (x) f (xk ) f '(xk )(x xk )
右端是直线方程,用这个直线方程来近似非线性方程 f (x) 。将非线性方程 f (x) 0 的根 x* 代入 f (x* ) 0 ,即 f (xk ) f '(xk )(x xk ) 0 解出 x* x f (xk ) f '(xk )
将右端取为 xk 1 ,则 xk 1 是比 xk 更接近于 x 的近似值,即
xk 1 xk
f (xk )
f '(xk )
这就是牛顿迭代公式,相应的迭代函数是
(x) x
f (x) f '(x)
2.2 牛顿迭代法的应用 计算 1 是求 cx 1 0 的解,解出 x ,即得到 1 。取 f (x) cx 1 , f '(x) c ,
c 有牛顿迭代公式
xk 1 x
k
c cxk 1 1 c c 这样就失去了迭代的意义,达不到迭代的效果。
故重新构造方程: cx2 x 0 , 1 也是该式的解。 故取 f (x) cx2 x , c
f '(x) 2cx 1 ,则有牛顿迭代公式 2
k k k cx 2 x cx 2 xk 1 xk 2cx
,
1 2c 1 k 0,1,...
1 的值在 1 k k ~ 1 之间,取初值 x 0.1。
18
3.流程图 20 10
0
1 k
k 1 k x1 x0
x0 f x0 f x0 x1
读入 x0 , , N
f x ?0 0
x x ? 1 0
输出3
3
4.输出结果 5.结果分析 当 k 3 时,得 5 位有效数字 0.05 564。此时, x3 x4 0.00 000 0.0 005 , 故取 x* x 0.05 564 0.056 。 此种迭代格式仍存在一定的缺陷,经实验后发现当初值 x0 x* 时必收敛,但 是当 x0 x* ( 0) 时迭代结果发散, 较小尚不确定。
6.心得体会 起初对题目的理解并不是很透彻,另外对构建牛顿迭代公式理论依据不是特 别充分,比如说为什么在原有直接得到的式子两边各乘一个 x,只是试出来的。 在编程方面不够成熟。当然也加深了对牛顿迭代法的理解和应用的具体实现。 4 kk
a 11 22 33
1.题目 实验二 例 3-4 用列主元消去法求解方程组 12x1 3x2 3x3
15
18x1 3x2 x3 15
x1 x2 x3 6
并求出系数矩阵 A 的行列式的值 det A 。 2.算法原理
2.1 顺序高斯消去法 顺序高斯消去法是利用线性方程组初等变换中的一种变换,即用一个不为零 的数乘一个方程后加至另一个方程,使方程组变成同解的上三角方程组,然后再 自下而上对上三角方程组求解。这样,顺序高斯消去法可分成“消去”和“回代” 两个过程。 在用顺序高斯消去法时,在消元之前检查方程组的系数矩阵的顺序主子式, 当阶数较高时是很难做到的。若线性方程组的系数具有某种性质时,如常遇到的 对角占优方程组,自然能够用高斯消去法求解。
2.2 列选主元消去法 线性方程组只要系数矩阵非奇异,就存在惟一解,但是按顺序消元过程中可 能出现主元素 a(k ) 0 ,这时尽管系数矩阵非奇异,消元过程无法再进行,或者
即使 (k ) 0
akk ,但如果其绝对值很小,用它作除数也会导致其他元素的数量级急剧
增大和使舍入误差扩大,将严重影响计算的精度。 为避免在校园过程确定乘数时的所用除数是零或绝对值小的数,即零主元或 小主元,在每一次消元之前,要增加一个选主元的过程,将绝对值大的元素交换 到主对角线的位置上来。 列选主元是当高斯消元到第 k 步时,从 k 列的 akk 以下(包括 akk )的各元素中 选出绝对值最大的,然后通过行交换将其交换到 akk 的位置上。交换系数矩阵中的 两行(包括常数项),只相当于两个方程的位置交换了,因此,列选主元不影响 求解的结果。 列选主元消去法常用来求行列式。设有矩阵
a11 L a
1n
A M M n1 L ann
用列主元消去法将其化为上三角形矩阵,对角线上元素为 a(1) , a(2) ,L, a(3) , 于是行 5
a a La akk d k l
aik i l
alj t, akj alj , t akj
j k , k 1,..., n
bl t, bk bl ,t bk
k 1 i aik
=
=
列式 det A (1)m (1) (2) (n)
11 22 nn
其中 m 为所进行的行交换次数。这是实际中求行列式值的可靠方法。
3.流程图
d ? d
k 1 k 1 k aik akk aik
i k 1, k 2,L, n
aij aik akj aij
i, j k 1, k 2,L, n
bi aibk bi
i k 1, k 2,L, n
bn ann bn
1 n
a b i a b b
ij j i
j i 1
i n 1, n 2,L,1
读入数据 aij , bi
i, j 1, 2,L, n
输出 b1 , b2 ,L, bn 6
kk 4.输出结果 5.结果分析 采用计算机运算在计算大数据时有明显的优点,另外也需要考虑到存储。 高斯消去法的使用条件是 a(k ) 0, k 1, 2,L, n ,而列选主元法可以保证这一 条件。并且可以避免在消元过程确定乘数时所用除数是绝对值小的数,相对全选 主元的运算量小,一般也可以满足精度要求。
6.心得体会
此次上机不仅需要对原理了解透彻,而且要求的编程能力较强。在定义和思 路上没问题,只是在编程软件的使用上遇到些不稳定的问题,如头文件的使用。 在存储空间上得到了新的认识,另外发现了当代码多时流程框图的好处。编程是 一件很需要耐心的事,自己还有很大进步空间。 7
(1) (1) (2)
M
1.题目 实验三 例 3-10 用杜里特尔分解法求矩阵 A 的逆矩阵 A1 。 1 1 A 1 2 12
2.算法原理
2 1 1
2.1 杜里特尔分解法
设线性方程组 Ax b ,对系数矩阵 A 进行除不交换两行位置得初等行变换相 当于用初等矩阵 M1 左乘 A ,在对方程组第一次消元后,A 和 b 分别化为 A 和
b(2) ,即 M A(1) A
(2)
1
M b(1) b
(2)
1
1
m21 1
其中 M1 m31 1
M O mn1 1第 k 次消元时, A(k ) 和 b(k ) 分别化为 A(k 1) 和 b(k 1) ,即
M A(k ) A
(k 1)
k
M b(k )
b(k 1)
k
1
其中 1
O 1 mk 1,k O
M O
mnk 1