数值分析上机题(matlab版)(东南大学)
数值分析上机题
![数值分析上机题](https://img.taocdn.com/s3/m/9cb09ddf76eeaeaad1f3305c.png)
数值分析上机题数值分析上机题写在前面:五道上机题分别为:习题2、习题3、习题4、习题5和习题6。
所有程序均在matlab 语言编写,并且在matlab2011a 的环境下编译通过,每个编程题均写了一个函数保存为函数名相同的m 文件。
可在matlab 命令窗口中输入函数名调用每个函数,查看程序运行结果。
习题220.(上机题)Newton 迭代法(1)给定初值0x 及容许误差ε,编制Newton 法解方程()0f x =根的通用程序。
(2)给定方程3()/30f x x x =-=,易知其有三个根13x *=-,20x *=,33x *=。
1.由Newton 方法的局部收敛性可知存在0δ>,当0(,)x δδ∈-时,Newton 迭代序列收敛于根2x *。
试确定尽可能大的δ。
2.试取若干初始值,观察当0(,1)x ∈-∞-,(1,)δ--,(,)δδ-,(,1)δ,(1,)∞时Newton 序列是否收敛以及收敛于哪一个根。
(3)通过本上机题,你明白了什么?******************************************************************************* 编程思路Newton 迭代法的基本原理如下:)()(1k k k k x f x f x x '-=+ 取初值0x ,通过不断迭代求取满足误差ε的根。
M 文件源代码:function Nfun(x0,error); k=1; syms xdy=subs(diff(x*x*x/3-x),x,x0); x1=x0-(x0*x0*x0/3-x0)/dy; while abs(x1-x0)>error x0=x1;dy=subs(diff(x*x*x/3-x),x,x0); x1=x0-(x0*x0*x0/3-x0)/dy; k=k+1;endfprintf('xk=%f\t k=%d\n',x1,k);%error为允许误差,这里取保留4位有效数字,0.000005在matlab命令窗口中输入Nfun函数名,且包含两个参数x0和error,x0是迭代初值,error 是允许误差。
东南大学数值分析上机题C语言编程源程序(前三章)
![东南大学数值分析上机题C语言编程源程序(前三章)](https://img.taocdn.com/s3/m/5b7a5a1ac5da50e2524d7f35.png)
第一章作业#include<stdio.h>const int N=100;int main(){int n;float S1,S2;double S3;S1=0;S2=0;for(n=2;n<N+1;n++){S1=S1+1.0/(n*n-1);}printf("the sum S1 is %d\n",S1);for(n=N;n>1;n--){S2=S2+1.0/(n*n-1);}printf("the sum S2 is %d\n",S2);S3=0.5*(1.5-1.0/N-1.0/(N+1));printf("S3 is %d\n",S3);}修改N的值即可第二章作业(1)float df(float x){float df;df=x*x-1;return df;}main(){float x0,x1,a;int k=0;printf("请输入初值x0=");scanf_s("%f",&x0);do{a=-f(x0)/df(x0);x1=x0+a;k++;x0=x1;}while(fabs(a)>eps);printf("迭代次数k=%d,根值x0=%f",k,x0);}(2)#include<stdio.h>#include<math.h>#define eps 0.01float f(float x){float f;f=x*x*x/3-x;return f;}float df(float x){float df;df=x*x-1;return df;}float ddf(float x){float ddf;ddf=2*x;return ddf;}main(){float x0,x1,a,dt=0;int k=0;do{dt=dt+eps;k++;}while((f(-dt)*f(dt)<0)&&(df(dt)!=0)&&(2*dt>=-f(-dt)/df(-dt))&&(2*dt>=f(dt)/df(dt)));printf("请输入x0:");scanf_s("%f",&x0);if(x0<-1)do{a=-f(x0)/df(x0);x1=x0+a;k++;x0=x1;}while(fabs(a)>eps);else if(x0>-1&&x0<-dt)do{a=-f(x0)/df(x0);x1=x0+a;k++;x0=x1;}while(fabs(a)>eps);else if(x0>-dt&&x0<dt)do{a=-f(x0)/df(x0);x1=x0+a;k++;x0=x1;}while(fabs(a)>eps);else if(x0>dt&&x0<1)do{a=-f(x0)/df(x0);x1=x0+a;k++;x0=x1;}while(fabs(a)>eps);else if (x0>1)do{a=-f(x0)/df(x0);x1=x0+a;k++;x0=x1;}while(fabs(a)>eps);printf("收敛域=%f,迭代次数=%d\n",dt,k);printf("根=%f\n",x0);}第三章作业#include<stdio.h>#include<math.h>main(){int n,i,j,k,t;float a[10][11],s[10],x[10],s2,max,sum,b;printf("请输入n= ");scanf("%d",&n);printf("请输入矩阵A= \n");for(i=0;i<=n;i++)for(j=0;j<=(n+1);j++)scanf("%f",&a[i][j]);//变换成高斯消去法矩阵//for(k=0;k<n;k++){//找出第k列中绝对值最大的一项//j=0;for(i=k;i<=n;i++){s[j]=a[i][k];j++;}max=fabs(s[0]);t=k;for(i=1;i<j;i++){if (fabs(s[i])>max){max=fabs(s[i]);t=k+i;}}// printf("t=%d\n",t);//将在列中有最大的元素的一行和第k行交换//for(j=0;j<=n+1;j++){b=a[k][j];a[k][j]=a[t][j];a[t][j]=b;}/*printf("***\n");for(i=0; i<=n; i++){for(j=0; j<=(n+1); j++)printf("%f ",a[i][j]);printf("\n");}*///将第k行k列的元素变为零//for(i=(k+1);i<=n;i++){s2=a[i][k]/a[k][k];for(j=k;j<=(n+1);j++)a[i][j]=a[i][j]-s2*a[k][j];}/*printf("###\n");for(i=0; i<=n; i++){for(j=0; j<=(n+1); j++)printf("%f ",a[i][j]);printf("\n");}*/printf("\n");}printf("\n");//解方程--//for(i=n;i>=0;i--){sum=0;x[n+1]=a[n][n+1]/a[n][n];for(j=i;j<=n;j++)sum+=a[i-1][j]*x[j+1];x[i]=(a[i-1][n+1]-sum)/a[i-1][i-1];}printf("方程的解为:");for(i=0;i<=n;i++)printf("%f,",x[i+1]);}。
上机数值计算练习题及答案.docx
![上机数值计算练习题及答案.docx](https://img.taocdn.com/s3/m/3922e93525c52cc58bd6beb9.png)
习题31、在MATLAB 中,先运行指令A=magic(3), B=[l,2,l;3,4,3;5,6,7], C=reshape(l:6,3,2)生成阵列A 珂,B3X2,C3X2 ,然后根据运行结果回答以下问题:(1)计算A*B,B*A,这两个乘积相同吗?(2)计算A\B, B/A,左除、右除结果相同吗?(3)计算B(:,[1,2]).*C和C.*B(:, [1,2]),这两个乘积相同吗?(4)计算A\A和AAA,这两个计算结果相同吗?(5)计算A\eye(3)和inv(A),这两个计算结果相同吗?(提示:根据对计算结果的目测回答问题)解答如下:运行指令:A=magic(3), B=[l,2,l;3,4,3;5,6,7], C=reshape(l:6,3,2)得到结果:8 1 63 5 74 9 2B =1 2 13 4 35 6 7C =1 42 53 6(1)计算A*B,并得到结果如下:» A*Bans =41 56 5353 68 6741 56 45计算B*A, 并得到结果如下:»B*Aans =18 20 2248 50 5286 98 86从以上计算结果可以得出结论:A*BJ (2)计算A\B ,并得到结果如下:» A\Bans =0.0333 0.1000 0.16110.5333 0.6000 0.74440.0333 0.1000 -0.1722计算B/A, 并得到结果如下:B/Aans =0.0056 0.0889 0.17220.1389 0.2222 0.30560.2333 0.7333 0.2333 与B*A这两个乘积不同。
从以上计算结杲可以得出结论:左除、右除结杲不同。
(3)计算B(:,[1,2]).弋,并得到结果如下:A =» B(:,[1,2]).*C ans =1 8 6 20 15 36计算C.*B(:, [1,2]),并得到结果如下: » CFB(:, [1,2]) ans =1 6 20 15 36从以上计算结果可以得出结论:B(: J1,2]).*C 和C ・*B(:, [1,2])的两个乘积相同。
东南大学数值分析上机实验题(下)
![东南大学数值分析上机实验题(下)](https://img.taocdn.com/s3/m/2b22853a5ef7ba0d4b733b83.png)
数值分析上机报告XX:学号:2013年12月22日第四章38.(上机题)3次样条插值函数(1)编制求第一型3次样条插值函数的通用程序;端点条件为'0y =0.8,'10y =0.2。
用所编制程序求车门的3次样条插值函数S(x),并打印出S(i+0.5)(i=0,1,…9)。
解:(1)#include <iostream.h> #include <math.h>floatx1[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<<"+"<<y; else if (y==1) cout<<"+"; else if (y<0) cout<<y; }void printY(float y){ if (y!=0) cout<<y; }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):"<<endl;}while ((n>100)||(n<1));cout<<"请输入Xi(i=0,1,...,"<<n<<"):";for (i=0;i<=n;i++) cin>>x1[i];cout<<endl;cout<<"请输入Yi(i=0,1,...,"<<n<<"n):";for (i=0;i<=n;i++) cin>>f1[i];cout<<endl;cout<<"请输入第一型边界条件Y'0,Y'n:";cin>>y_0>>y_n;cout<<endl;for (i=0;i<n;i++) h[i]=x1[i+1]-x1[i];for (i=1;i<n;i++) u[i]=(float) (h[i-1]/(h[i-1]+h[i]));for (i=1;i<n;i++) w[i]=(float) (1.0-u[i]);for (i=0;i<n;i++) f2[i]=(f1[i+1]-f1[i])/h[i]; //一阶差商for (i=0;i<n-1;i++) f3[i]=(f2[i+1]-f2[i])/(x1[i+2]-x1[i]); //二阶差商for (i=1;i<n;i++) d[i]=6*f3[i-1]; //求出所有的d[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<n;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<n;i++) //计算S(x)的表达式c[i]=(float) (f2[i]-(1.0/3.0*m[i]+1.0/6.0*m[i+1])*h[i]);for (i=0;i<n;i++)e[i]=(m[i+1]-m[i])/(6*h[i]);for (i=0;i<n;i++){cout<<"X属于区间["<<x1[i]<<","<<x1[i+1]<<"]时"<<endl<<endl;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<<endl<<endl;}cout<<"针对本题,计算S(i+0.5)(i=0,1,...,9):"<<endl;for (i=0;i<10;i++){if ((i+0.5<=x1[n])&&(i+0.5>=x1[0])){calculation((float)(i+0.5));cout<<"S("<<i+0.5<<")="<<s<<endl;}else cout<<i+0.5<<"超出定义域"<<endl;};cout<<endl;}(2)编制的程序求车门的3次样条插值函数S(x):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(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.7237S(6.5)=5.59441 S(7.5)=5.42989 S(8.5)=5.65976 S(9.5)=5.7323第六章21.(上机题)常微分方程初值问题数值解(1)编制RK4方法的通用程序;(2)编制AB4方法的通用程序(由RK4提供初值);(3)编制AB4-AM4预测校正方法的通用程序(由RK4提供初值);(4)编制带改进的AB4-AM4预测校正方法的通用程序(由RK4提供初值);(5)对于初值问题h=,应用(1)~(4)中的四种方法进行计算,并将计算结果和精确解取步长0.13=+作比较;y x x()3/(1)(6)通过本上机题,你能得到哪些结论?解:#include<iostream.h>#include<fstream.h>#include<stdlib.h>#include<math.h>ofstream outfile("data.txt");//此处定义函数f(x,y)的表达式//用户可以自己设定所需要求得函数表达式double f1(double x,double y){double f1;f1=(-1)*x*x*y*y;return f1;}//此处定义求函数精确解的函数表达式double f2(double x){double f2;f2=3/(1+x*x*x);return f2;}//此处为精确求函数解的通用程序void accurate(double a,double b,double h){double x[100],accurate[100];x[0]=a;int i=0;outfile<<"输出函数准确值的程序结果:\n";do{x[i]=x[0]+i*h;accurate[i]=f2(x[i]);outfile<<"accurate["<<i<<"]="<<accurate[i]<<'\n';i++;}while(i<(b-a)/h+1);}//此处为经典Runge-Kutta公式的通用程序void RK4(double a,double b,double h,double c) {int i=0;double k1,k2,k3,k4;double x[100],y[100];y[0]=c;x[0]=a;outfile<<"输出经典Runge-Kutta公式的程序结果:\n"; do{x[i]=x[0]+i*h;k1=f1(x[i],y[i]);k2=f1((x[i]+h/2),(y[i]+h*k1/2));k3=f1((x[i]+h/2),(y[i]+h*k2/2));k4=f1((x[i]+h),(y[i]+h*k3));y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;outfile<<"y"<<"["<<i<<"]="<<y[i]<<'\n';i++;}while(i<(b-a)/h+1);}//此处为4阶Adams显式方法的通用程序void AB4(double a,double b,double h,double c) {double x[100],y[100],y1[100];double k1,k2,k3,k4;y[0]=c;x[0]=a;outfile<<"输出4阶Adams显式方法的程序结果:\n";for(int i=0;i<=2;i++){x[i]=x[0]+i*h;k1=f1(x[i],y[i]);k2=f1((x[i]+h/2),(y[i]+h*k1/2));k3=f1((x[i]+h/2),(y[i]+h*k2/2));k4=f1((x[i]+h),(y[i]+h*k3));y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;}int j=3;y1[0]=y[0];y1[1]=y[1];y1[2]=y[2];y1[3]=y[3];do{x[j]=x[0]+j*h;y1[j+1]=y1[j]+(55*f1(x[j],y1[j])-59*f1(x[j-1],y1[j-1])+37*f1(x[j-2],y1[j-2])-9*f1(x[j-3],y1[ j-3]))*h/24;outfile<<"y1"<<"["<<j<<"]="<<y1[j]<<'\n';j++;}while(j<(b-a)/h+1);}//主函数void main(void){double a,b,h,c;cout<<"输入上下区间、步长和初始值:\n";cin>>a>>b>>h>>c;accurate(a,b,h);RK4(a,b,h,c);AB4(a,b,h,c);}float si(int u,int v){float sum=0; int q;for(q=0;q<k;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<n;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<n;t++)a[k][t]=si(k,t);for(t=(k+1);t<m;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<m;z++)sum+=a[t][z]*x[z];x[t]=(float)(a[t][m]-sum)/a[t][t];}}(5)由经典Runge-Kutta公式得出的结果列在下面的表格中,以及精确值y(x i)和精确值和数值解的误差:由AB4方法得出的结果为:Y1[0]=3 y1[1]=2.997 y1[2]=2.97619 y1[3]=2.92113 y1[4]=2.81839 y1[5]=2.66467 y1[6]=2.4652 y1[7]=2.23308 y1[8]=1.98495 y1[9]=1.73704 y1[10]=1.50219 y1[11]=1.28876 y1[12]=1.10072 y1[13]=0.93871 y1[14]=0.801135y1[15]=0.685335(6)本次上机作业让我知道了在遇到复杂问题中,无法给出解析式的情况下,要学会灵活使用各种微分数值解法,并且能计算出不同方法的精度大小。
数值分析上机题Matlab(东南大学)3
![数值分析上机题Matlab(东南大学)3](https://img.taocdn.com/s3/m/dff07fbec77da26925c5b03d.png)
0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 0.38 0.40 0.42 0.44 0.46 0.48 0.50 0.52 0.54 0.56 0.58 0.60 0.62 0.64 0.66 0.68 0.70 0.72
152 139 128 119 110 103 96 90 85 80 76 72 68 65 62 59 56 53 51 49 47 45 43 41 39 38
========================================================================================================================
======================================================================================================================================================================== 习题 3_36 ======================================================================================================================================================================== Omega n x1 x2 x3 x4 x5 x6 x7 x8 x9
-0.71279 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71280 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281 -0.71281
数值计算上机实习题目(matlab编程)
![数值计算上机实习题目(matlab编程)](https://img.taocdn.com/s3/m/92abafa0d5d8d15abe23482fb4daa58da0111c98.png)
数值计算上机实习题目(matlab编程)非线性方程求根一、实验目的本次实验通过上机实习,了解迭代法求解非线性方程数值解的过程和步骤。
二、实验要求1、用迭代法求方程230x x e -=的根。
要求:确定迭代函数?(x),使得x=?(x),并求一根。
提示:构造迭代函数2ln(3)x ?=。
2、对上面的方程用牛顿迭代计算。
3、用割线法求方程3()310f x x x =--=在02x =附近的根。
误差限为410-,取012, 1.9x x ==。
三、实验内容1、(1)首先编写迭代函数,记为iterate.mfunction y=iterate(x)x1=g(x); % x 为初始值。
n=1;while(abs(x1-x)>=1.0e-6)&(n<=1000) % 迭代终止的原则。
x=x1;x1=g(x);n=n+1;endx1 %近似根n %迭代步数(2)后编制函数文件?(x),记为g.mfunction y=g(x)y=log(3*x.^2);(3)设初始值为0、3、-3、1000,观察初始值对求解的影响。
将结果记录在文档中。
>>iterate(0)>>iterate(3) 等等2、(1)首先编制牛顿迭代函数如下,记为newton.mfunction y=newton(x0)x1=x0-fc(x0)/df(x0); % 牛顿迭代格式n=1;while(abs(x1-x0)>=1.0e-6)&(n<=1000000) % 迭代终止的原则。
x0=x1;x1=x0-fc(x0)/df(x0);n=n+1;endx1 %近似根n %迭代步数(2)对题目中的方程编制函数文件,记为fc.mfunction y=fc(x)y=3*x.^2-exp(x)编制函数的导数文件,记为df.mfunction y=df(x)y=6*x-exp(x)(3)在MATLAB 命令窗计算,当设初始值为0时,newton(0);给定不同的初始值,观察用牛顿法求解时所需要的迭代步数,并与上面第一题的迭代步数比较。
matlab上机练习题答案(可编辑修改word版)
![matlab上机练习题答案(可编辑修改word版)](https://img.taocdn.com/s3/m/29963dc8581b6bd97e19ea3b.png)
a ⎣ ⎦ 1. 计算 a = 5⎦ >> a=[6 9 3;2 7 5]; >> b=[2 4 1;4 6 8]; >> a.*b ans =⎣4 8⎦>> d=deconv([3 13 6 8],[1 4]) d =312⎡2 4 7 4⎤⎡8⎤1236 36 求欠定方程组⎢9 3 5 6⎥ x = ⎢5⎥ 的最小范数解84240⎡4 9 2⎤⎡37⎤⎣ >> a=[2 4 7 4;9 3 5 6]; >> b=[8 5]'; ⎦ ⎣ ⎦2. 对于 AX = B ,如果 A = ⎢7 6 4⎥ , B = ⎢26⎥ ,求解 X 。
>> x=pinv(a)*b⎢ ⎥ ⎢ ⎥x =>> A=[4 9 2;7 6 4;3 5 7]; >> B=[37 26 28]’; >> X=A\B X =⎢⎣3 5 7⎥⎦ ⎢⎣28⎥⎦-0.2151 0.4459 0.7949 0.27077 用符号函数法求解方程 a t 2+b*t +c=0-0.5118 >> r=solve('a*t^2+b*t+c=0','t') 4.0427 r =1.3318[ 1/2/a*(-b+(b^2-4*a*c)^(1/2))] ⎡1 2 5 ⎤ ⎡8 - 7 4⎤ [ 1/2/a*(-b-(b^2-4*a*c)^(1/2))]3. a = ⎢3 6 - 4⎥ , b = ⎢3 6 2⎥ ,观察 a 与 b 之间的⎣ ⎦ ⎣ ⎦⎡a 11 a 12 ⎤六种关系运算的结果 >> a=[1 2 3;4 5 6];8 求矩阵 A = ⎢ ⎣ 21 ⎥ 的行列式值、逆和特征根a 22 ⎦>> b=[8 –7 4;3 6 2];>> syms a11 a12 a21 a22;>> a>b >> A=[a11,a12;a21,a22] ans =>> AD=det(A) % 行列式 0 1 0 >> AI=inv(A) % 逆 11>> AE=eig(A) % 特征值 >> a>=b ans =0 1 0 1 01>> a<b ans =1 0 1 0 1>> a<=b ans =1 0 1 010 A = [ a11, a12][ a21, a22] AD =a11*a22-a12*a21 AI =[ -a22/(-a11*a22+a12*a21), a12/(-a11*a22+a12*a21)] [ a21/(-a11*a22+a12*a21), -a11/(-a11*a22+a12*a21)] AE = [1/2*a11+1/2*a22+1/2*(a11^2- 2*a11*a22+a22^2+4*a12*a21)^(1/2)] [1/2*a11+1/2*a22-1/2*(a11^2->> a==b2*a11*a22+a22^2+4*a12*a21)^(1/2)] ans =9 因式分解: x 4 - 5x 3 + 5x 2 + 5x - 60 0 0 >> syms x;0 00 >> f=x^4-5*x^3+5*x^2+5*x-6; >> a~=b ans =>> factor(f) ans =(x-1)*(x-2)*(x-3)*(x+1)4 计算多项式乘法(x 2+2x +2)(x 2+5x +4) ⎡10 f = ⎢ x 21 x ⎤⎥ ,用符号微分求 df/dx 。
Matlab上机练习一答案
![Matlab上机练习一答案](https://img.taocdn.com/s3/m/46fbb9cbac51f01dc281e53a580216fc700a53fe.png)
Matlab上机练习⼀答案Matlab 上机练习⼀班级学号姓名按要求完成题⽬,并写下指令和运⾏结果。
(不需要画图)1、求下列联⽴⽅程的解=+-+-=-+=++-=--+41025695842475412743w z y x w z x w z y x w z y x >> a=[3 4 -7 -12;5 -7 4 2;1 0 8 -5;-6 5 -2 10];>> b=[4;4;9;4];>> c=a\bc =5.22264.45701.47181.59942、设 ??++=)1(sin 35.0cos 2x x x y ,把x=0-2π间分为101点,画出以x 为横坐标,y 为纵坐标的曲线。
>> x=linspace(0,2*pi,101);>> y=cos(x)*(0.5+(1+x.^2)\3*sin(x));>> plot(x,y,'r')3、产⽣8×6阶的正态分布随机数矩阵R1, 求其各列的平均值、和、乘积。
并求该矩阵全体数的平均值。
(mean var )a=randn(8,6)mean(a)var(a)k=mean(a)k1=mean(k)i=ones(8,6)i1=i*k1i2=a-i1i3=i2.*i2g=mean(i3)g2=mean(g)或者u=reshape(a,1,48);p1=mean(u)p2=var(u)4、绘制)(222y x e x z +-=在定义域x=[-2,2],y=[-2,2]内的曲⾯。
(利⽤meshgrid )x=-2:2;y=x;[X,Y]= meshgrid(x,y);Z=X^2*exp(-(X^2+Y^2));mesh(X,Y,Z)5、求代数⽅程3x 5+4x 4+7x 3+2x 2+9x+12=0的所有根。
(利⽤roots 函数)p=[3 4 7 2 9 12];roots(p)6、把1开五次⽅,并求其全部五个根。
东南大学 数值分析上机题作业 MATLAB版
![东南大学 数值分析上机题作业 MATLAB版](https://img.taocdn.com/s3/m/69272dfd4a7302768f9939c5.png)
东南大学数值分析上机题作业MATLAB版上机作业题报告1.Chapter 11.1题目设S N= Nj=2j2−1,其精确值为11311(-- 。
22N N +1(1)编制按从大到小的顺序S N =(2)编制按从小到大的顺序S N =111,计算S N 的通用程序。
++⋯⋯+22-132-1N 2-1111,计算S N 的通用程++⋯⋯+N 2-1(N -1 2-122-1序。
(3)按两种顺序分别计算S 102, S 104, S 106, 并指出有效位数。
(编制程序时用单精度)(4)通过本次上机题,你明白了什么?1.2程序 1.3运行结果1.4结果分析按从大到小的顺序,有效位数分别为:6,4,3。
按从小到大的顺序,有效位数分别为:5,6,6。
可以看出,不同的算法造成的误差限是不同的,好的算法可以让结果更加精确。
当采用从大到小的顺序累加的算法时,误差限随着N 的增大而增大,可见在累加的过程中,误差在放大,造成结果的误差较大。
因此,采取从小到大的顺序累加得到的结果更加精确。
2.Chapter 22.1题目(1)给定初值x 0及容许误差ε,编制牛顿法解方程f(x=0的通用程序。
3(2)给定方程f (x =x-x =0, 易知其有三个根x 1*=-3, x 2*=0, x 3*=○1由牛顿方法的局部收敛性可知存在δ>0, 当x 0∈(-δ, Newton 迭代序列收敛+δ 时,于根x2*。
试确定尽可能大的δ。
○2试取若干初始值,观察当x 0∈(-∞, -1, (-1, -δ, (-δ, +δ, (δ, 1, (1, +∞ 时Newton 序列的收敛性以及收敛于哪一个根。
(3)通过本上机题,你明白了什么?2.2程序2.3运行结果(1)寻找最大的δ值。
算法为:将初值x0在从0开始不断累加搜索精度eps ,带入Newton 迭代公式,直到求得的根不再收敛于0为止,此时的x0值即为最大的sigma 值。
最新东南大学数值分析上机题matlab(前三章)
![最新东南大学数值分析上机题matlab(前三章)](https://img.taocdn.com/s3/m/4f1fcfeb3169a4517623a304.png)
数值分析上机题第一章(17题)(1)从2依次累加到N的程序function sn = sum1( n )sn=0;sn=single(sn);for i=2:nai=1/(i^2-1);sn=sn+ai;endend(2)从N依次累加到2的程序function sn = sum2( n )sn=0;sn=single(sn);for i=n:-1:2ai=1/(i^2-1);sn=sn+ai;endend(3)编制求精确值的求和函数sum0function sn = sum0( n )sn=0;sn=single(sn);sn=1/2*(3/2-1/n-1/(n+1));end按第一种顺序得到的值及有效位数如下:N=100时sn0=sum0(100);sn=sum1(100)n=fix(-log10(2*abs(sn-sn0)))得到:sn =0.7400495 n =7N=10e4时sn0=sum0(10e4);sn=sum1(10e4)n=fix(-log10(2*abs(sn-sn0)))得到:sn =0.7498521 n =3N=10e6时sn0=sum0(10e6);sn=sum1(10e6)n=fix(-log10(2*abs(sn-sn0)))得到:sn =0.7498521 n =3按第二种顺序得到的值及有效位数如下:N=100时sn0=sum0(100);sn=sum2(100)n=fix(-log10(2*abs(sn-sn0)))得到:sn =0.7400495 n =7N=10e4时sn0=sum0(10e4);sn=sum2(10e4)n=fix(-log10(2*abs(sn-sn0)))得到:sn =0.7499900 n =7N=10e6时sn0=sum0(10e6);sn=sum2(10e6)n=fix(-log10(2*abs(sn-sn0)))得到:sn =0.7499999 n =7(4)通过这道上机题,我明白了,应用计算机进行求和运算时,求和的顺序不同对结果的精度是有影响的。
数值分析上机题 舍入误差与有效数
![数值分析上机题 舍入误差与有效数](https://img.taocdn.com/s3/m/90cff658804d2b160b4ec0ea.png)
(5)通过上述分析可以看出:按从小到大的顺序计算所得的结果与真值接近,而按从大到小的顺序计算所得的结果与真值的误差较大,且有效位数较前者少。
原因:这是由于机器数在进行加法运算时,首先比较两数的阶码,将阶码较小的尾数向右移位,每移一位阶码加一,直至其阶码与另一数的阶码一致为止,且将移位后的尾数多于计算机字长的部分进行四舍五入,之后对尾数进行加减运算,最后将尾数写成规格化的形式,当从大到小的顺序进行计算式,由于越到后面数字越小,就会产生大数吃小数的情况,从而产生误差的累积,最后使计算结果的不准确。
解:(1)从大到小的matlab程序:
functions=myfun1(N)
formatlong;
k=2;
s=single(0);
fork=2:1:N
a=1/(k*k-1);
s=a+s;
end
end
(2)从小到大的matlab程序
functions=myfun2(N)
formatlong;
s=single(0);
fori=N:-1:2
a=1/(i*i-1);
s=a+s;
end
真值
有效位数
0.7400495
0.7400495
0.7400495
大小
7位
小大
7位
0.7498521
0.7499000
0.7499000
大小
4位
小大
7位
0.7498521
0.7499990
0.749999
大小
4位
小大
舍入误差与有效数
东南大学机械工程学院
设SN= ,其精确值为 )。
(1)编制按从大到小的顺序SN= + + +……+ ,计算SN通用的程序;
Matlab上机练习题及答案
![Matlab上机练习题及答案](https://img.taocdn.com/s3/m/82e6e5704b73f242336c5ff1.png)
Matlab 上机练习题及答案---------------------------------------------------------------------1、 矩阵Y= ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡3472123100451150425,给出元素1的全下标和单下标,并用函数练习全下标和单下标的转换,求出元素100的存储位置。
取出子矩阵⎥⎦⎤⎢⎣⎡21301,并求该矩阵的维数。
解:命令为:Y=[5,2,4;0,15,1;45,100,23;21,47,3] Y(2,3) Y(10)sub2ind([4 3],2,3)[i,j]=ind2sub([4 3],10)find(Y==100) sub2ind([4 3],3,2)B=Y(2:2:4,3:-2:1) 或 B=Y([2 4],[3 1]) [m n]=size(Y)---------------------------------------------------------------------2、已知矩阵A=[1 0 -1 ;2 4 1; -2 0 5],B=[0 -1 0;2 1 3;1 1 2] 求2A+B 、A 2-3B 、A*B 、B*A 、A .*B ,A/B 、A\B解:命令为:A=[1 0 -1 ;2 4 1; -2 0 5] B=[0 -1 0;2 1 3;1 1 2] E=2*A+B F=A^2-3*B G=A*B H=B*A I=A.*B J=A/B K=A\B---------------------------------------------------------------------3、利用函数产生3*4阶单位矩阵和全部元素都为8的4*4阶矩阵,并计算两者的乘积。
解:命令为: A=eye(3,4) B=8*ones(4)C=A*B---------------------------------------------------------------------4、创建矩阵a=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------7023021.5003.120498601,取出其前两列构成的矩阵b ,取出前两行构成矩阵c ,转置矩阵b 构成矩阵d ,计算a*b 、c<d ,c&d, c|d ,~c|~d解:命令为:a=[-1,0,-6,8;-9,4,0,;0,0,,-2;0,-23,0,-7] b=a(:,[1 2]) c=a([1 2],:) d=b ’ e=a*b f=c<d g=c&d h=c|d i=~c|~d---------------------------------------------------------------------5、求!201∑=n n解:命令文件为 sum=0; s=1;for n=1:20 s=n*s; sum=sum+s; end sum---------------------------------------------------------------------6、求a aa aaa aa a S n ++++=得值,其中a 是一个数字,由键盘输入,表达式中位数最多项a 的个数,也由键盘输入。
东南大学_数值分析_第六章_常微分方程数值解法
![东南大学_数值分析_第六章_常微分方程数值解法](https://img.taocdn.com/s3/m/38118bd8b14e852458fb57d5.png)
第六章 常微分方程数值解法——RK 4法、AB 4法******(学号) *****(姓名)上机题目要求见教材P307,23题。
一、算法原理题目要求采用RK 4法和AB 4法求解最简单的常微分方程初值问题(,),()y f x y a x by a η'=≤≤⎧⎨=⎩ (1)为求解式(1),采用离散化方法,就是寻求解)(x y 在区间],[b a 上的一系列点<<<<<n x x x x 321上的近似值 ,,,,21n y y y 。
记1(1,2,)i i i h x x i -=-=表示相邻两个节点的间距,称为步长。
求微分方程数值解的主要问题:(1) 如何将微分方程(,)y f x y '=离散化,并建立求其数值解的递推公式; (2) 递推公式的局部截断误差、数值数n y 与精确解)(n x y 的误差估计; (3) 递推公式的稳定性与收敛性. a) Runge-Kutta 方法基本思想:通过在1[,]i i x x +多预报几个点求斜率,并将其加权平均作为k *的近似值,以此构造更高精度的计算公式。
如果每步计算四次函数 的值,完全类似的,可以导出局部截断误差为)(5h O 的四阶Runge-Kutta 公式(RK 4):1123412132431(22),6(,),(,),221(,),22(,).n n n n n n n n n n y y k k k k k f x y h h k f x y k h k f x h y k k f x h y hk +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ (2)b) Adams 显式公式Runge-Kutta 方法是单步法,计算1+n y 时,只用到n y , 而已知信息1-n y 、2-n y 等没有被直接利用。
可以设想如果充分利用已知信息1-n y ,2-n y ,…来计算1+n y ,那么不但有可能提高精度,而且大大减少了计算量,这就是构造所谓线性多步法的基本思想。
数值分析作业-matlab上机作业
![数值分析作业-matlab上机作业](https://img.taocdn.com/s3/m/13860268a45177232f60a2ab.png)
数值分析———Matlab上机作业学院:班级:老师:姓名:学号:第二章解线性方程组的直接解法第14题【解】1、编写一个追赶法的函数输入a,b,c,d输出结果x,均为数组形式function x=Zhuiganfa(a,b,c,d)%首先说明:追赶法是适用于三对角矩阵的线性方程组求解的方法,并不适用于其他类型矩阵。
%定义三对角矩阵A的各组成单元。
方程为Ax=d%b为A的对角线元素(1~n),a为-1对角线元素(2~n),c为+1对角线元素(1~n-1)。
% A=[2 -1 0 0% -1 3 -2 0% 0 -2 4 -3% 0 0 -3 5]% a=[-1 -2 -3];c=[-1 -2 -3];b=[2 3 4 5];d=[6 1 -2 1];n=length(b);u(1)=b(1);y(1)=d(1);for i=2:nl(i)=a(i-1)/u(i-1);%先求l(i)u(i)=b(i)-c(i-1)*l(i);%再求u(i)%A=LU,Ax=LUx=d,y=Ux,%Ly=d,由于L是下三角矩阵,对角线均为1,所以可求y(i)y(i)=d(i)-l(i)*y(i-1);endx(n)=y(n)/u(n);for i=(n-1):-1:1%Ux=y,由于U是上三角矩阵,所以可求x(i)x(i)=(y(i)-c(i)*x(i+1))/u(i);end2、输入已知参数>>a=[2 2 2 2 2 2 2];>>b=[2 5 5 5 5 5 5 5];>>c=[2 2 2 2 2 2 2];>>d=[220/27 0 0 0 0 0 0 0];3、按定义格式调用函数>>x=zhuiganfa(a,b,c,d)4、输出结果x=[8.147775166909105 -4.073701092835030 2.036477565178471 -1.017492820111148 0.507254485099400 -0.250643392637350 0.119353996493976 -0.047741598597591]第15题【解】1、编写一个程序生成题目条件生成线性方程组A x=b 的系数矩阵A 和右端项量b ,分别定义矩阵A 、B 、a 、b 分别表示系数矩阵,其中1(10.1;,1,2,...,)j ij i i a x x i i j n -==+=或1(,1,2,...,)1ij a i j n i j ==+-分别构成A 、B 对应右端项量分别a 、b 。
东南大学数值分析编程作业
![东南大学数值分析编程作业](https://img.taocdn.com/s3/m/46a0d6b3c77da26925c5b071.png)
数值分析编程作业******学号:******指导老师:***目录P20.第17题:舍入误差与有效数 (3)P56.第20题:Newton迭代法 (5)P126.第39题:列主元Gauss消去法 (13)P127.第40题:逐次超松弛迭代法 (16)P195.第37题:3次样条插值函数 (22)P256.第23题:重积分的计算 (24)P307.第23题:常微分方程初值问题数值解 (26)P346.第10题:抛物方程Crank-Nicolson格式 (31)注:源程序放在对应的文件夹里,如第一题放在“P20.第17题:舍入误差与有效数”文件夹里;P20.第17题:舍入误差与有效数设2211311,--122N N+1NN j S j ==-∑其精确值为()。
(1)编制按从小到大的顺序222111+++2-13-1N -1N S =⋅⋅⋅⋅⋅,计算N S 的通用程序; (2)编制按从小到大的顺序222111+++N -1N-1-12-1N S =⋅⋅⋅⋅⋅(),计算N S 的通用程序; (3)按两种顺序分别计算210S ,410S ,610S ,并指出有效位数(编制程序时用单精度); (4)通过本上机题你明白了什么?答 (1)(2)程序:程序放在computation.m, goon.m 文件里,运行程序时只需输入命令main(3)运行结果:N=100时,S1有6位有效数字,S2都是有效数字N=10000时,S1有4位有效数字,S2都是有效数字N=1000000时,S1有3位有效数字,S2都是有效数字(4)体会:当N从小到大变化时,SN的项从大到小,反之SN的项从小到大。
从运算结果可见,第一种计算结果总是比第2、3种计算结果小,这是由于大数吃小数的原因(计算机表示一个数值是用字节表示的,当程序进行计算时,先将第一项和第二项相加,然后再加第三项……所以加到后面会由于项的值很小,不能加到这个字节里,所以就被忽略了。
(整理)东南大学信息学院年matlab上机考试题及答案
![(整理)东南大学信息学院年matlab上机考试题及答案](https://img.taocdn.com/s3/m/11c1ac69964bcf84b9d57bf8.png)
MATLAB 上机测验题(考试时间:2:20----4:20)姓名 学号考试要求:1、要求独立完成不得与他人共享,答卷雷同将做不及格处理。
2、答卷用Word 文件递交,文件名为学号+姓名.doc ,试卷写上姓名及学号。
3、答卷内容包括:(1) 程序;(2) 运行结果及其分析;(3) 图也要粘贴在文档中。
上机考题: 一、系统传递函数为1121()10.81z H z z z---+=-+,按照以下要求求解: 1) 求其极零点图,判断系统的稳定性,画出系统的频谱特性;2) 当系统输入信号为:()[5cos(0.2)2sin(0.7)]x n n n ππ=++,050n ≤≤时,画出系统的输出。
(1) 代码如下:clcclearb=[1 1 0];a=[1,-1,0.81];sys=tf(b,a,-1);figurepzmap(sys);saveas(gcf,'p1_1','bmp');figurew=0:0.1:20;freqz(b,a,w);saveas(gcf,'p1_2','bmp');%第二问开始figure;n=0:50;x=5+cos(0.2*pi*n)+2*sin(0.7*pi*n); lsim(sys,x)saveas(gcf,'p1_3','bmp');运行结果如下:图1-1图1-2系统极零图如图1-1所示,通过该离散系统的极零图可以看出,极点全部都在单位圆内,所以系统是稳定的。
系统的频谱特性曲线如图1-2所示(2)运行结果如下图所示图1-3图1-3为系统在x(n)激励下的输出,其中灰色部分为输入信号,蓝色为输出信号。
二、系统传递函数为324321130()9459750s s sH ss s s s++=++++,1)画出系统的零极点图,判断稳定性;2)给定频率范围为[0,10],步长为0.1,画出其频率响应;3)画出系统的单位脉冲响应。
数值分析实验(参考答案)
![数值分析实验(参考答案)](https://img.taocdn.com/s3/m/044a7d6d561252d380eb6eae.png)
数值分析上机实验学生姓名: 学号: 教师:实验1:1. 实验项目的性质和任务通过上机实验,使学生对病态问题、线性方程组求解和函数的数值逼近方法有一个初步理解。
2.教学内容和要求 1)对高阶多多项式201()(1)(2)(20)()k p x x x x x k ==---=-∏编程求下面方程的解 19()0p x x ε+=并绘图演示方程的解与扰动量ε的关系。
Matlab 程序:x=1:20;y=zeros(1,20); ve=zeros(1,21); plot(x,y,'o') hold on ; pause; for x=1:5ve(2)=10^(-x);e=roots(poly(1:20)+ve);plot(e,'*') hold on pause; end0510152025303540-20-15-10-5051015202)对2~20n =,生成对应的Hilbert 矩阵,计算矩阵的条件数;通过先确定解获得常向量b 的方法,确定方程组 n H x b =最后,用矩阵分解方法求解方程组,并分析计算结果。
Matlab 程序:for n=2:20; H=hilb(n); ca=cond(H,2) x=(1:n)'; b=H*x; [L,U]=lu(H); y=L\b;x1=U\yplot(x,x,'o',x1,x1,'*') hold on pause; end-500-400-300-200-100100200300400500-500-400-300-200-10001002003004005003)对函数 21()[1,1]125f x x x=∈-+的Chebyshev 点 (21)cos()1,2,...,12(1)kk x k n n π-==++编程进行Lagrange 插值,并分析插值结果。
Matlab 程序:function y=lagrangen(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i);s=0;for k=1:nL=1;for j=1:nif j~=kL=L*(z-x0(j))/(x0(k)-x0(j));endends=s+L*y0(k);endy(i)=s;endy;for n=5:20x=-1:0.01:1; y=1./(1+25*x.^2);plot(x,y,'r')hold onk=n+1:-1:1;x0=cos((2*k-1)*pi./(2*(n+1))),y0=1./(1+25*x0.^2);x=-1:0.01:1; y1=lagrangen(x0,y0,x);plot(x0,y0,'o',x,y1), pausehold off end-1-0.8-0.6-0.4-0.200.20.40.60.81-0.200.20.40.60.811.2-1-0.8-0.6-0.4-0.200.20.40.60.81-1-0.8-0.6-0.4-0.200.20.40.60.81。
东南大学信息短学期matlab作业一
![东南大学信息短学期matlab作业一](https://img.taocdn.com/s3/m/2aa88f66561252d380eb6e48.png)
实验一、>> A=[1,2,3;4,5,6;7,8,9]A =1 2 34 5 67 8 9>> B=[9,8,7;6,5,4;3,2,1],C=[4,5,6;7,8,9;1,2,3]B =9 8 76 5 43 2 1C =4 5 67 8 91 2 3>> AA =1 2 34 5 67 8 9>> BB =9 8 76 5 43 2 1>> CC =4 5 67 8 91 2 3三、>> A=[1,2,3]A =1 2 3>> B=[4,5,6]B =4 5 6>> C=A+BC =5 7 9>> D=A-BD =-3 -3 -3>> E=A.*BE =4 10 18>> F=A./BF =0.2500 0.4000 0.5000 >> G=A.^BG =1 32 729>> A=[1,2,3;4,5,6;7,8,9]A =1 2 34 5 67 8 9>> inv(A)ans =1.0e+016 *-0.4504 0.9007 -0.45040.9007 -1.8014 0.9007-0.4504 0.9007 -0.4504 >>四、2、>> A=[1,1,1;1,-2,1;1,2,3];>> b=[2;-1;-1];>> x=inv(A)*bx =3.00001.0000-2.00004、>> A=[2,3,-1;3,-2,1;1,2,1];>> b=[18;8;24];>> x=inv(A)*bx =468>>实验二二、>> k=-10:10;>> y=sin(k).*(0*(k<0)+1*(k>0));>> stem(k,y)>> k=-10:10;>> y=sin(k)+sin(pi*k);>> stem(k,y)>> k=-10:10;>> y=k.*sin(k).*(0*(k<3)+1*(k>3));>> stem(k,y)>> k=-10:10;>> y=(-1).^k+((0.5).^k).*(0*(k<0)+1*(k>0));>> stem(k,y)七、>> f1=[1,1,1,1];>> f2=[3,2,1];>> y=conv(f1,f2)y =3 5 6 6 3 1>> f1=[0,0,1,1,1];>> f2=[0,0,0.5,0.25,0.125];>> y=conv(f1,f2)y =Columns 1 through 80 0 0 0 0.5000 0.7500 0.8750 0.3750 Column 90.1250>> f1=[0,0,1,1,1];>> f2=[0,0,-0.5,0.25,-0.125];>> y=conv(f1,f2)y =Columns 1 through 80 0 0 0 -0.5000 -0.2500 -0.3750 0.1250Column 9-0.1250实验三、二、1、x=[1,2,3,4,5,6,6,5,4,3,2,1];dtft=zeros(100);for i=1:100w(i)=(i-50)/10;for k=1:12dtft(i)=dtft(i)+x(k)*exp(-j*(k-1)*w(i));endendsubplot(1,2,1);plot(w,abs(dtft),'r-');xlabel('w'); ylabel('DTFT'); title('幅频特性');subplot(1,2,2);plot(w,angle(dtft),'r-');xlabel('w'); ylabel('DTFT');title('相频特性');2、X=[1,2,3,4,5,6,6,5,4,3,2,1];F=fftshift(fft(X,1000));k=[-500:499]*2*pi/1000;subplot(1,2,1);plot(k,abs(F),'b-')xlabel('w');ylabel('DFT');title('幅频特性');subplot(1,2,2);plot(k,angle(F) ,'b-')xlabel('w');ylabel('DFT');title('相频特性');3、x=[1,2,3,4,5,6,6,5,4,3,2,1];dtft=zeros(100);for i=1:100w(i)=(i-50)/10;for k=1:12dtft(i)=dtft(i)+x(k)*exp(-j*(k-1)*w(i));endendsubplot(1,2,1);plot(w,abs(dtft),'r-');xlabel('w'); title('幅频特性'); hold;subplot(1,2,2);plot(w,angle(dtft),'r-');xlabel('w');title('相频特性'); hold;F=fftshift(fft(X,1000));k=[-500:499]*2*pi/1000;subplot(1,2,1);plot(k,abs(F),'b-')hold;subplot(1,2,2);plot(k,angle(F) ,'b-')三、1、%f1.mfunction y=fun1(x)if((-pi<x) && (x<0))y=pi+x;elseif ((0<x) && (x<pi))y=pi-x;elsey=0endfor i=1:1000g(i)=f1(2/1000*i-1);w(i)=(i-1)*0.2*pi;endfor i=1001:10000g(i)=0;w(i)=(i-1)*0.2*pi;endG=fft(g)/1000;subplot(1,2,1);plot(w(1:50),abs(G(1:50)));xlabel('w');ylabel('G');title('DFT幅度频谱1'); subplot(1,2,2);plot(w(1:50),angle(G(1:50)))xlabel('w');ylabel('Fi');title('DFT相位频谱1');2、%f2.mfunction y=fun2(x)if x<1 && x>-1y=cos(pi*x/2);elsey=0;endfor i=1:1000g(i)=f2(2/1000*i-1);w(i)=(i-1)*0.2*pi;endfor i=1001:10000g(i)=0;w(i)=(i-1)*0.2*pi;endG=fft(g)/1000;subplot(1,2,1);plot(w(1:50),abs(G(1:50)));xlabel('w');ylabel('G');title('DFT幅度频谱2'); subplot(1,2,2);plot(w(1:50),angle(G(1:50)))xlabel('w');ylabel('Fi');title('DFT相位频谱2');3、%f3.mfunction y=fun3(x)if x<0 && x>-1y=1;elseif x>0 && x<1y=-1;elsey=0endfor i=1:1000g(i)=f3(2/1000*i-1);w(i)=(i-1)*0.2*pi;endfor i=1001:10000g(i)=0;w(i)=(i-1)*0.2*pi;endG=fft(g)/1000;subplot(1,2,1);plot(w(1:100),abs(G(1:100)));xlabel('w');ylabel('G');title('DFT幅度频谱3'); subplot(1,2,2);plot(w(1:100),angle(G(1:100)))xlabel('w');ylabel('Fi');title('DFT相位频谱3');。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
数值分析上机报告第一章一、题目精确值为)11123(21+--N N 。
1) 编制按从大到小的顺序11131121222-+⋯⋯+-+-=N S N ,计算S N 的通用程序。
2) 编制按从小到大的顺序1211)1(111222-+⋯⋯+--+-=N N S N ,计算S N 的通用程序。
3) 按两种顺序分别计算64210,10,10S S S ,并指出有效位数。
(编制程序时用单精度) 4) 通过本次上机题,你明白了什么?二、通用程序三、求解结果四、结果分析可以得出,算法对误差的传播又一定的影响,在计算时选一种好的算法可以使结果更为精确。
从以上的结果可以看到从大到小的顺序导致大数吃小数的现象,容易产生较大的误差,求和运算从小数到大数算所得到的结果才比较准确。
第二章一、题目(1)给定初值0x 及容许误差ε,编制牛顿法解方程f(x)=0的通用程序。
(2)给定方程03)(3=-=x x x f ,易知其有三个根3,0,3321=*=*-=*x x xa) 由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时,Newton 迭代序列收敛于根x 2*。
试确定尽可能大的δ。
b)试取若干初始值,观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。
(3)通过本上机题,你明白了什么?二、通用程序1.运行search.m 文件结果为:The maximum delta is 0.774597即得最大的δ为0.774597,Newton 迭代序列收敛于根*2x =0的最大区间为(-0.774597,0.774597)。
2.运行Newton.m 文件在区间(,1),(1,),(,),(,1),(1,)δδδδ-∞----++∞上各输入若干个数,计算结果如下:区间(,1)-∞-上取-1000,-100,-50,-30,-10,-8,-7,-5,-3,-1.5x。
结果显示,以上初值迭代序列均收敛于-1.732051,即根*1在区间(1,)δ--即区间(-1,-0.774597)上取-0.774598,-0.8,-0.85,-0.9,-0.99,计算结果如下:计算结果显示,迭代序列局部收敛于-1.732051,即根*1x ,局部收敛于1.730251,即根*3x 。
在区间(,)δδ-即区间(-0.774597,0.774597)上,由search.m 的运行过程表明,在整个区间上均收敛于0,即根*2x 。
在区间(,1)δ即区间(0.774597,1)上取0.774598,0.8,0.85,0.9,0.99,计算结果如下:计算结果显示,迭代序列局部收敛于-1.732051,即根1x ,局部收敛于1.730251,即根3x 。
上取100,60,20,10,7,6,4,3,1.5,计算结果如下: 区间(1,)x。
结果显示,以上初值迭代序列均收敛于1.732051,即根*3综上所述:(-∞,-1)区间收敛于-1.73205,(-1,δ)区间局部收敛于 1.73205,局部收敛于-1.73205,(-δ,δ)区间收敛于0,(δ,1)区间类似于(-1,δ)区间,(1,∞)收敛于1.73205。
通过本上机题,明白了对于多根方程,Newton法求方程根时,迭代序列收敛于某一个根有一定的区间限制,在一个区间上,可能会局部收敛于不同的根。
第三章一、题目列主元Gauss 消去法对于某电路的分析,归结为求解线性方程组RI V =。
其中3113000100001335901100000931100000000107930000900030577050000074730000000030410000005002720009000229R --⎛⎫ ⎪--- ⎪ ⎪-- ⎪--- ⎪⎪=--- ⎪-- ⎪ ⎪- ⎪-- ⎪ ⎪--⎝⎭()15,27,23,0,20,12,7,7,10T T V =----(1) 编制解n 阶线性方程组Ax b =的列主元高斯消去法的通用程序;(2) 用所编程序线性方程组RI V =,并打印出解向量,保留5位有效数;二、通用程序%% 列主元Gauss 消去法求解线性方程组%%%%参数输入n=input('Please input the order of matrix A: n='); %输入线性方程组阶数n b=zeros(1,n);A=input('Input matrix A (such as a 2 order matrix:[1 2;3,4]) :'); b(1,:)=input('Input the column vector b:'); %输入行向量bb=b';C=[A,b]; %得到增广矩阵%%列主元消去得上三角矩阵for i=1:n-1 [maximum,index]=max(abs(C(i:n,i)));index=index+i-1;T=C(index,:);C(index,:)=C(i,:);C(i,:)=T;for k=i+1:n %%列主元消去if C(k,i)~=0C(k,:)=C(k,:)-C(k,i)/C(i,i)*C(i,:);endendend%% 回代求解 %%x=zeros(n,1);x(n)=C(n,n+1)/C(n,n);for i=n-1:-1:1x(i)=(C(i,n+1)-C(i,i+1:n)*x(i+1:n,1))/C(i,i);endA=C(1:n,1:n); %消元后得到的上三角矩阵disp('The upper teianguular matrix is:')for k=1:nfprintf('%f ',A(k,:));fprintf('\n');enddisp('Solution of the equations:');fprintf('%.5g\n',x); %以5位有效数字输出结果以教材第123页习题16验证通用程序的正确性。
执行程序,输入系数矩阵A 和列向量b,结果如下:结果与精确解完全一致。
三、求解结果执行程序,输入矩阵A(即题中的矩阵R)和列向量b(即题中的V),得如下结果:Please input the order of matrix A: n=9Input matrix A (such as a 2 order matrix:[1 2;3,4]):[31 -13 0 0 0 -10 0 0 0-13 35 -9 0 -11 0 0 0 00 -9 31 -10 0 0 0 0 00 0 -10 79 -30 0 0 0 -90 0 0 -30 57 -7 0 -5 00 0 0 0 -7 47 -30 0 00 0 0 0 0 -30 41 0 00 0 0 0 -5 0 0 27 -20 0 0 -9 0 0 0 -2 29]Input the column vector b:[-15 27 -23 0 -20 12 -7 7 10]31.000000 -13.000000 0.000000 0.000000 0.000000 -10.000000 0.000000 0.000000 0.000000 0.000000 29.548387 -9.000000 0.000000 -11.000000 -4.193548 0.000000 0.000000 0.000000 0.000000 0.000000 28.258734 -10.000000 -3.350437 -1.277293 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 75.461271 -31.185629 -0.451999 0.000000 0.000000 -9.000000 0.000000 0.000000 0.000000 0.000000 44.602000 -7.179695 0.000000 -5.000000 -3.577994 0.000000 0.000000 0.000000 0.000000 -0.000000 45.873193 -30.000000 -0.784718 -0.561543 0.000000 0.000000 0.000000 0.000000 -0.000000 -0.000000 21.380698 -0.513187 -0.367236 0.000000 0.000000 0.000000 0.000000 -0.000000 -0.000000 0.000000 26.413085 -2.419996 0.000000 0.000000 0.000000 0.000000 -0.000000 -0.000000 0.000000 0.000000 27.389504 Solution of the equations:-0.289230.34544-0.71281-0.22061-0.43040.15431-0.0578230.201050.29023由上述结果得:第四章一、题目二、通用程序三、求解结果1、数据输入Input n: n=10Input x:[0 1 2 3 4 5 6 7 8 9 10]Input y:[2.51 3.30 4.04 4.70 5.22 5.54 5.78 5.40 5.57 5.70 5.80] Input the derivative of y(0):0.5Input the derivative of y(n):0.22、计算结果第五章一、题目二、通用程序三、运行结果第六章一、题目二、通用程序1、RK4方法的通用程序2、AB4方法的通用程序3、AB4- AB4预测校正方法的通用程序4、带改进的AB4- AB4预测校正方法的通用程序三、结果比较四、结论。