数值分析——编程作业..docx
数值分析作业(C++)
数值分析上机题Chapter 1（1）从大到小顺序#include<iostream>#include<math.h>#include<string>using namespace std;void main(){int N;cout<<"Please input N:"<<endl;cin>>N;cout<<"The number you input is:N="<<N<<endl;float sn=0;for(int i=2;i<=N;++i)sn+=1.0/(i*i-1);cout<<"The result you want is:"<<endl<<"Sn="<<sn<<endl; }（2）从小到大的顺序#include<iostream>#include<math.h>#include<string>using namespace std;void main(){int N;cout<<"Please input N:"<<endl;cin>>N;cout<<"The number you input is:N="<<N<<endl;float sn=0;for(int i=N;i>1;--i)sn+=1.0/(i*i-1);cout<<"The result you want is:"<<endl<<"Sn="<<sn<<endl; }（3）结果：N 100 10000 1000000(1) 0.740049 0.749852 -14.2546(2) 0.74005 0.7499 -14.2551 Chapter 2(1) 通用程序：#include <iostream>#include <string>#include <math.h>using namespace std;double h;//允许误差double x[100];double x0;//初值double xl;//所求结果double f(double x) //原函数{double a;//函数值a=f(x);return a;}double f1(double x)//导函数{double b=0.01,c[100];d;for(int i=0; ;++i){c[i]=(f(x+b)-f(x))/b;//导函数if(i>=1){if(fabs(c[i]-c[i-1])<=h){d=c[i];break;}}b/=10;}return d;}void main(){for(int k=0; ;++i){x[0]=x0;//初值x[k+1]=x[k]-f(x[k])/f1(x[k]);if(fabs(x[k+1]-x[k])<=h)break;}cout<<"xl="<<x[k+1]<<endl;//输出结果}（2）#include <iostream>#include <string>#include <math.h>using namespace std;int k;double h=0.001;//允许误差double x[100];double x0;//初值double xl;//所求结果double f(double x) //原函数{double a;//函数值a=x*x*x/3-x;return a;}double f1(double x)//导函数{double b=0.01,c[100],d;for(int i=0; ;++i){c[i]=(f(x+b)-f(x))/b;//导函数if(i>=1){if(fabs(c[i]-c[i-1])<=h){d=c[i];break;}}b/=10;}return d;}int main(){for(int i=1;i<=1000;++i){x0=i/1000.0;for(k=0; ;++k){x[0]=x0;x[k+1]=x[k]-f(x[k])/f1(x[k]);if(fabs(x[k+1]-x[k])<=h)break;}if(fabs(x[k+1])>=1)break;}cout<<"when deta="<<x0<<": "<<"xl*="<<x[k+1]<<endl;//输出结果cout<<"This is the biggest deta,or the result will not be zero."<<endl;return 0;}Chapter 3(1).通用程序：#include <iostream>#include <string>#include <math.h>#define n 3//定义阶次using namespace std;void main(){int i,j;double a[n][n+1],temp,x[n];cout<<"请输入需求解矩阵:"<<endl;//输入求解矩阵for(i=0;i<n;++i)for(j=0;j<=n;++j)cin>>a[i][j];for(j=0;j<n;++j)//列主元{ int k=j;double h=a[j][j];//列首元int b=j;for(i=j+1;i<n;++i)//列最大元素并标记{if(fabs(a[i][k])>fabs(h)){h=a[i][k];b=i;}}for(int m=j;m<=n;++m)//交换主元{temp=a[k][m];a[k][m]=a[b][m];a[b][m]=temp;}for(int c=j+1;c<n;++c)//化为上三角矩阵for(int z=j;z<=n;++z){double yz=a[c][j]/h;a[c][z]=a[c][z]-a[j][z]*yz;}}//求解过程x[n-1]=a[n-1][n]/a[n-1][n-1];cout<<"x"<<"["<<n-1<<"]="<<x[n-1]<<endl;for(int e=n-2;e>=0;--e){ double l=0;for(int w=e+1;w<n;++w){l=l+a[e][w]*x[w];}x[e]=(a[e][n]-l)/a[e][e];cout<<"x"<<"["<<e<<"]="<<x[e]<<endl;}}x0=-0.23645,x1=0.47743,x2=-0.76698,x3=-0.077649,x4=-0.30792,x5=0.14634,x6=-0.17 073,x7=0.28480,x8=0.34483;Chapter 437.(1)通用程序程序在使用前要先设定n的值#include<iostream>#include<math.h>#include<string>using namespace std;#define n 10double x[n+1],y[n+1],h[n],u[n],u_u[n],d[n+1];double a[n+1][n+1],m[n+1];double F0,Fn,arr;double b[n+1][n+1];double f(double &a){double function;for(int i=0;i<n+1;++i)if(a==x[i])function=y[i];return function;}double f1(double &a,double &b){if(a!=b)return (f(b)-f(a))/(b-a);else{if(a==x[0])return F0;elsereturn Fn;}}double f2(double &a,double &b,double &c){return (f1(b,c)-f1(a,b))/(c-a);}double s(double x0){double sx;if(x0<=x[0]||x0>=x[n]){cout<<"不?在¨²定¡§义°?域®¨°内¨²"<<endl;return 0;}else{int j;for(int i=0;i<n+1;++i)if(x0>=x[i]&&x0<=x[i+1]) j=i;sx=y[j]+(f1(x[j],x[j+1])-(m[j]/3.0+m[j+1]/6.0)*h[j])*(x0-x[j])+(m[j]/2.0)*(x0-x[j])*(x0-x [j])+((m[j+1]-m[j])/6*h[j])*(x0-x[j])*(x0-x[j])*(x0-x[j]);return sx;}}int main(){cout<<"please input x:";for(int i1=0;i1<n+1;++i1)cin>>x[i1];cout<<endl<<endl;cout<<"please input y:";for(int i2=0;i2<n+1;++i2)cin>>y[i2];cout<<endl<<endl;cout<<"please input F0,Fn:";cin>>F0>>Fn;cout<<endl<<endl;for(int i3=0;i3<n;++i3)h[i3]=x[i3+1]-x[i3];for(int i4=1;i4<n;++i4){u[i4]=h[i4-1]/(h[i4-1]+h[i4]);u_u[i4]=1-u[i4];}for(int i=0;i<n+1;++i)for(int j=0;j<n+1;++j){if(i==j)a[i][j]=2;else if(j==i-1){if(i!=n)a[i][j]=u[i];elsea[i][j]=1;}else if(j==i+1){if(i!=0){a[i][j]=u_u[i];}elsea[i][j]=1;}elsea[i][j]=0;cout<<a[i][j]<<" ";if(j==n)cout<<endl<<endl;}cout<<"d[i]的Ì?值¦Ì：êo"<<endl;for(int i5=0;i5<n+1;++i5){if(i5==0)d[i5]=6*f2(x[0],x[0],x[1]);else if(i5==n)d[i5]=6*f2(x[n-1],x[n],x[n]);elsed[i5]=6*f2(x[i5-1],x[i5],x[i5+1]); cout<<d[i5]<<" ";}cout<<endl<<endl;cout<<"系¦Ì数ºy矩?阵¨®：êo"<<endl<<endl; for(int i=0;i<n+1;++i)for(int j=0;j<=n+1;++j){if(j<n+1)b[i][j]=a[i][j];elseb[i][j]=d[i];cout<<b[i][j]<<" ";if(j==n+1)cout<<endl;}cout<<endl<<endl;for (int i=1;i<=n;i++){ arr=b[i][i-1];for (int j=0;j<=n+1;j++)b[i][j]=b[i][j]-a[i-1][j]*arr/b[i-1][i-1]; }m[n]=b[n][n+1]/b[n][n];for(int i=n-1;i>=0;i--) m[i]=(b[i][n+1]-b[i][i+1]*m[i+1])/b[i][i];for(int i=0;i<n+1;++i)cout<<m[i]<<" ";cout<<endl<<endl;for(int j=0;j<n;++j){cout<<"当Ì¡Àx在¨²区?间?["<<x[j]<<","<<x[j+1]<<"]时º¡À:";cout<<"S(x)="<<y[j]<<"+"<<(f1(x[j],x[j+1])-(m[j]/3.0+m[j+1]/6.0)*h[j]) <<"(x-"<<x[j]<<")+"<<m[j]/2.0<<"(x-"<<x[j]<<")^2+"<<(m[j+1]-m[j])/6*h[j]<<"(x-"<<x[j]<<")^3"<<endl<<endl;}for(int i=0;i<10;++i)cout<<"s("<<i<<"+0.5)="<<s(i+0.5)<<" "<<endl;return 0;}（2）此时n为10，将程序中的n设为10，输入数据运行，得到结果：Chapter 523（1）通用程序程序在运行之前要先设定积分区间，即a,b,c,d的值，本题为了编程方便，已经把数据写入程序中。
数值分析编程作业(jacobi,sor)
编程作业一、第三章数值求解正方形区域的poisson 方程边值问题。
{−(∂2u ∂x 2+∂2u ∂y 2)=f(x ，y)=1,0<x,y <1u (0,y )=u (1,y )=0u (x,0)=u (x,1)=x(1−x)写成五点差分格式得到线性方程组为4u i,j −u i−1,j −u i+1,j −u i,j−1−u i,j+1=h 2f ij ,1≤i,j ≤N u 0,j =u 1,j =0,u i,0=u i,1=(i −1)∗h ∗(1−(i −1)∗h)写成矩阵形式Au=f 。
其中v 1=(u 1,1,u 2,1,…,u N,1)T,v 2=(u 1,2,u 2,2,…,u N,2)T,……,v 1=(u 1,N ,u 2,N ,…,u N,N )Tb 1=h 2(f 1,1,f 2,1,f 3,1,…,f N,1)T,b 2=h 2(f 1,2,f 2,2,f 3,2,…,f N,2)T,…,b N =h 2(f 1,N ,f 2,N ,f 3,N ,…,f N,N )Th =1N +1,取N =99，h =0.01.f i,j =2,i,j =1,2,…,N1、用jacobi 迭代法求解方程组 程序：function[u,k]=jaco(n) h=1/(n+1);f(2:n+1,2:n+1)=2*h*h; %h2*f u=zeros(n+2,n+2); v=zeros(n+2,n+2);for i=2:n+1 %边界条件u(i,1)=(i-1)*h*(1-(i-1)*h); u(i,n+2)= (i-1)*h*(1-(i-1)*h); ende=0.000000001;for k=1:1000 %迭代求解 er=0;for i=2:n+1for j=2:n+1v(i,j)=(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1)+f(i,j))/4; er=er+abs(v(i,j)-u(i,j)); %估计误差1122NN A I I A A I IA -⎛⎫ ⎪-⎪= ⎪- ⎪-⎝⎭1122N N v b v b u f v b ⎛⎫⎛⎫ ⎪ ⎪ ⎪⎪== ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭4114114ii A -⎛⎫ ⎪- ⎪= ⎪- ⎪-⎝⎭endendif er/n^2<e,break;end %判断是否达到计算精度，如果达到则退出循环u(2:n+1,2:n+1)=v(2:n+1,2:n+1);end计算结果：[u,k]=jaco(9)u=0 0 0 0 0 0 0 0 0 0 0 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0 0 0 0 0 0 0 0 0 0 0 k =3182、块jacobi迭代法求解程序：function [u,k]=kjac(n)% xsbj：用块Jacobi迭代法求解线性方程组A*u=f% u：方程组的解； k：迭代次数； n：非边界点数；% f：线性方程组A*u=f的右端矩阵f；% a：方程组系数矩阵的下对角线元素； b：方程组系数矩阵的主对角线元素；% c：方程组系数矩阵的上对角线元素； d：追赶法所求方程的右端向量；% e：允许误差界； er：迭代误差；h=1/(n+1);f(2:n+1,2:n+1)=2*h*h;a=-1*ones(1,n);b=4*ones(1,n);c=-1*ones(1,n);u=zeros(n+2,n+2);v=zeros(n+2,n+2);e=0.000000001;for i=2:n+1u(i,1)=(i-1)*h*(1-(i-1)*h);u(i,n+2)=(i-1)*h*(1-(i-1)*h);endfor i=2:n+1f(i,1)=2*h*h+(i-1)*h*(1-(i-1)*h);f(i,n+2)=2*h*h+(i-1)*h*(1-(i-1)*h);endfor k=1:2000er=0;for j=2:n+1d(1:n)=f(2:n+1,j)+u(2:n+1,j-1)+u(2:n+1,j+1) ;x=zg(a,b,c,d); % 用追赶法求解v(2:n+1,j)=x';er=er+norm(v(:,j)-u(:,j),1);endif er/n^2<e,break;endu(2:n+1,2:n+1)=v(2:n+1,2:n+1);end计算结果：[u,k,t]=kjac(9)u =0 0 0 0 0 0 0 0 0 0 0 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0 0 0 0 0 0 0 0 0 0 0 k =1703、SOR方法求解方程组程序：function [u,k,t]=sor(w,n)h=1/(n+1);f(2:n+1,2:n+1)=2*h*h;u=zeros(n+2,n+2);for i=2:n+1u(i,1)=(i-1)*h*(1-(i-1)*h);u(i,n+2)=(i-1)*h*(1-(i-1)*h);ende=0.000000001;ticfor k=1:1000er=0;for j=2:n+1for i=2:n+1du(i,j)=w*(-4*u(i,j)+u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1)+f(i,j))/4;u(i,j)=u(i,j)+du(i,j);er=er+abs(du(i,j));endendif er/n^2<e,break;endend计算结果：[u,k]=sor(1.5,9)u=0 0 0 0 0 0 0 0 0 0 0 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2500 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2400 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.2100 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.1600 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0.0900 0 0 0 0 0 0 0 0 0 0 0 k =46二、第五章插值比较程序：x=(1200:400:3600);y=(1200:400:3600);xx=(1200:100:3600);yy=(1200:100:3600);z=[1130 1250 1280 1230 1040 900 5001320 1450 1420 1400 1300 700 9001390 1500 1500 1400 900 1100 10601500 1200 1100 1350 1450 1200 11501500 1200 1100 1550 1600 1550 13801500 1550 1600 1550 1600 1600 16001480 1500 1550 1510 1430 1300 1200];[xx,yy]=meshgrid(xx,yy);zz=Interp2(x,y,z,xx,yy);mesh(xx,yy,zz);结果：三、计算二重积分%首先以0.01*0.01的方格把求积区间等分%每个小区间采用Simpson公式来求积h=0.01; %小正方形的边长（也可以是矩形）Xa=1.5;Xb=2;Nx=(Xb-Xa)/h; %Nx等分Yc=0;Yd=3;Ny=(Yd-Yc)/h; %Ny等分x(1,1:Nx+1)=(1.5:h:2);y(1,1:Ny+1)=(0:h:3);A=0;for v=1:Nxfor w=1:Ny[S]=Simpson(x(v),x(v+1),y(w),y(w+1));A=A+S;endenddouble A采用Simpson公式计算矩形区域的二重积分function [S]=Simpson(a,b,c,d) %注意一般要求a<x<b, c<y<d C=[1/6 4/6 1/6]; %牛顿科特斯系数x=[a 0.5*(a+b) b];y=[c 0.5*(c+d) d];S=0;for i=0:2he=0;for k=0:2he=he+C(k+1)*log(x(i+1)+2*y(k+1));endS=S+C(i+1)*he;endS=(b-a)*(d-c)*S;计算结果：A=65。
数值分析编程
习题一：已知A=⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛-------------------------------00001.40446782.4213474.2238417.0784317.1037585.0010103.1123124.3031743.2446782.4719334.30123789.3103456.1121314.271828.0336993.0112384.1067813.3213474.2123789.37138465.14103458.0111223.3101023.2056781.1784165.0742382.1238417.0103456.1103458.0789365.9431637.0741856.3836742.1563849.1718719.0784317.1121314.2111223.3431637.0897918.19101011.4031454.2189736.2113584.0037585.071828.0101023.2741856.3101011.4108437.27123848.3012345.1112336.1010103.1336993.0056781.1741856.3031454.2123848.3567914.15125432.3061074.1123124.3112348.1784165.0836742.1189736.2012345.1125432.3141823.19115237.2031743.2067813.3742382.1718719.0113584.0112336.1061074.1115237.238412.12b=(2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392);1． 用Household 变换，把A 化为三对角阵B 。
数值分析——编程作业
《数值分析》实验报告第二章：解线性方程组的直接方法2、试用MATLAB软件编程实现追赶法求解三对角方程组的算法，并考虑梯形电阻电路问题，电路如下：其中电路中的各个电流{1i ，2i ，…，8i }须满足下列线性方程组：R V i i =- 22 210 252321=-+-i i i 0 252 432=-+-i i i 0 252 543=-+-i i i 0 252 654=-+-i i i 0 252 765=-+-i i i 0 252 876=-+-i i i 052 87=+-i i设V 220=V ，Ω=27R ，运用求各段电路的电流量。
解：1481.827220≈=R V 上述方程组可用矩阵表示为：⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--------------00000001481.8522520000002520000002520000002520000002520000002520000002287654321i i i i i i i iMatLab 程序： %赋初值；a=[0 -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=[8.1481 0 0 0 0 0 0 0]; %三对角方程的追赶法for i=2:8%“追”的过程； a(i)=a(i)/b(i-1);b(i)=b(i)-c(i-1)*a(i); d(i)=d(i)-a(i)*d(i-1); end;d(8)=d(8)/b(8);%“赶”的过程； for i=7:-1:1d(i)=(d(i)-c(i)*d(i+1))/b(i); end;x=d; x程序运行结果：x =8.1477 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477即(A) 0477.01194.02506.05073.00175.10365.20737.41477.887654321⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=i i i i i i i i I 。
《数值分析》课程设计—作业实验一...
《数值分析》课程设计—作业实验一1.1 水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。
由于旅途的颠簸，大家都很疲惫，很快就入睡了。
第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。
第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题。
解：一、问题分析：对于本题，比较简单，我们只需要判断原来椰子的个数及每个人私藏了一份之后剩下的是否能被5除余1，直到最后分完。
二、问题求解：通过matlab 建立M 文件，有如下程序：或者对于第一个程序，n 取2000；对于第二个程序，n 取20001，就能得到我们想要的结果，即原先一共有15621个椰子，最终平均每人得4092个椰子。
n=input('input n:');forx=1:n p=5*x+1;for k=1:5 p=5*p/4+1;end if p==fix(p) break ; end enddisp([x,p]) input n:20001023 15621function fentao(n)a=cat(1,7);for j=n:-1:1 a(1)=j;i=1; while i<7a(i+1)=4*(a(i)-1)/5; i=i+1;endif a(7)==fix(a(7)) a, end end end>> fentao(20001) a =15621 12496 9996 7996 6396 5116 4092（本文档内的有些运行结果，限于篇幅，使文档结构更和谐、紧凑，已做相关的改动，程序代码没变）1.2 当0,1,2,,100n = 时，选择稳定的算法计算积分1d 10n xn xe I x e--=+⎰.解：一、问题分析：由1d 10n xn xe I x e--=+⎰知： 110101==+⎰dx I I以及： )1(110101011)1(1nnxxnxxn n n endx edx ee eI I ----+-+-==++=+⎰⎰得递推关系：⎪⎩⎪⎨⎧--=-=-+n n n I e n I I I 10)1(1101101，但是通过仔细观察就能知道上述递推公式每一步都将误差放大十倍，即使初始误差很小，但是误差的传播会逐步扩大，也就是说用它构造的算法是不稳定的，因此我们改进上述递推公式（算法）如下：⎪⎪⎩⎪⎪⎨⎧--=-=+-))1(1(101)1(101110n n n I e n I I I通过比较不难得出该误差是逐步缩小的，即算法是稳定的。
(完整word版)同济大学数值分析matlab编程题汇编
MATLAB 编程题库 1.下面的数据表近似地满足函数21cxbax y ++=，请适当变换成为线性最小二乘问题，编程求最好的系数c b a ,,，并在同一个图上画出所有数据和函数图像.625.0718.0801.0823.0802.0687.0606.0356.0995.0628.0544.0008.0213.0362.0586.0931.0ii y x ----解：x=[-0.931 -0.586 -0.362 -0.213 0.008 0.544 0.628 0.995]'; y=[0.356 0.606 0.687 0.802 0.823 0.801 0.718 0.625]'; A=[x ones(8,1) -x.^2.*y]; z=A\y;a=z(1); b=z(2); c=z(3); xh=-1:0.1:1;yh=(a.*xh+b)./(1+c.*xh.^2); plot(x,y,'r+',xh,yh,'b*')2.若在Matlab工作目录下已经有如下两个函数文件，写一个割线法程序，求出这两个函数10 的近似根，并写出调用方式：精度为10解：>> edit gexianfa.mfunction [x iter]=gexianfa(f,x0,x1,tol)iter=0;while(norm(x1-x0)>tol)iter=iter+1;x=x1-feval(f,x1).*(x1-x0)./(feval(f,x1)-feval(f,x0));x0=x1;x1=x;end>> edit f.mfunction v=f(x)v=x.*log(x)-1;>> edit g.mfunction z=g(y)z=y.^5+y-1;>> [x1 iter1]=gexianfa('f',1,3,1e-10)x1 =1.7632iter1 =6>> [x2 iter2]=gexianfa('g',0,1,1e-10)x2 =0.7549iter2 =83.使用GS 迭代求解下述线性代数方程组：123123123521242103103x x x x x x x x x ì++=-ïïïï-++=íïïï-+=ïî解：>> edit gsdiedai.mfunction [x iter]=gsdiedai(A,x0,b,tol) D=diag(diag(A)); L=D-tril(A); U=D-triu(A); iter=0; x=x0;while((norm(b-A*x)./norm(b))>tol) iter=iter+1; x0=x;x=(D-L)\(U*x0+b); end>> A=[5 2 1;-1 4 2;1 -3 10]; >> b=[-12 10 3]'; >>tol=1e-4; >>x0=[0 0 0]';>> [x iter]=gsdiedai(A,x0,b,tol); >>x x =-3.0910 1.2372 0.9802 >>iter iter = 64.用四阶Range-kutta 方法求解下述常微分方程初值问题（取步长h=0.01）,(1)2x dy y e xy dx y ìïï=++ïíïï=ïî解：>> edit ksf2.mfunction v=ksf2(x,y) v=y+exp(x)+x.*y;>> a=1;b=2;h=0.01; >> n=(b-a)./h; >> x=[1:0.01:2]; >>y(1)=2;>>fori=2:(n+1)k1=h*ksf2(x(i-1),y(i-1));k2=h*ksf2(x(i-1)+0.5*h,y(i-1)+0.5*k1); k3=h*ksf2(x(i-1)+0.5*h,y(i-1)+0.5*k2); k4=h*ksf2(x(i-1)+h,y(i-1)+k3); y(i)=y(i-1)+(k1+2*k2+2*k3+k4)./6; end >>y调用函数方法>> edit Rangekutta.mfunction [x y]=Rangekutta(f,a,b,h,y0) x=[a:h:b]; n=(b-a)/h; y(1)=y0; fori=2:(n+1)k1=h*(feval(f,x(i-1),y(i-1)));k2=h*(feval(f,x(i-1)+0.5*h,y(i-1)+0.5*k1)); k3=h*(feval(f,x(i-1)+0.5*h,y(i-1)+0.5*k2)); k4=h*(feval(f,x(i-1)+h,y(i-1)+k3)); y(i)=y(i-1)+(k1+2*k2+2*k3+k4)./6; end>> [x y]=Rangekutta('ksf2',1,2,0.01,2); >>y5.取0.2h =，请编写Matlab 程序，分别用欧拉方法、改进欧拉方法在12x ≤≤上求解初值问题。
(完整版)数值分析第一次作业
问题1：20.给定数据如下表：试求三次样条插值S(x)，并满足条件 (1)S`(0.25)=1.0000，S`(0.53)=0.6868； （2）S ’’(0.25)=S ’’(0.53)=0。
分析：本问题是已知五个点，由这五个点求一三次样条插值函数。
边界条件有两种，（1）是已知一阶倒数，（2）是已知自然边界条件。
对于第一种边界(已知边界的一阶倒数值)，可写出下面的矩阵方程。
⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡432104321034322110d M M M M M 200020000020022d d d d λμμλμλμλ其中μj =j1-j 1-j h h h +，λi=j1-j j h h h +，dj=6f[x j-1,x j ,x j+1]， μn =1，λ0=1对于第一种边界条件d 0=0h 6（f[x 0,x 1]-f 0`），d n =1-n h 6（f`n-f `[x n-1,x n ]） 解：由matlab 计算得：由此得矩阵形式的线性方程组为：⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡ 2.1150-2.4286-3.2667-4.3143-5.5200-M M M M M 25714.00001204286.000004000.026000.0006429.023571.0001243210解得 M 0=-2.0286;M 1=-1.4627;M 2= -1.0333; M 3= -0.8058; M 4=-0.6546S(x)=⎪⎪⎩⎪⎪⎨⎧∈-+-+-∈-+-+-∈-+-+-∈-+-+-]53.0,45.0[x 5.40x 9.1087x 35.03956.8.450-x 1.3637-x .5301.67881- ]45.0,39.0[x 9.30x 11.188x 54.010.418793.0-x 2.2384-x .450(2.87040-]39.0,30.0[x 03.0x 6.9544x 9.30 6.107503.0-x 1.9136-x .3902.708779-]30.0,25.0[x 5.20x 10.9662x 0.3010.01695.20-x 4.8758-x .3006.76209-33333333），（）（）（）（），（）（）（）），（）（）（）（），（）（）（）（Matlab 程序代码如下：function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。
(完整word版)数值分析作业(C语言编程实现)(word文档良心出品)
#include <stdio.h>#include <math.h>double f(double x){double ans;ans=exp(x);return ans;}void main(){double a=1,b=3,error=0.0001,t[20][20],h,c;int i,j,k,m,n;h=b-a;t[0][0]=h*(f(a)+f(b))/2;k=1;while(1){t[0][k]=0;m=1;for(j=0;j<k-1;j++)m=m*2;for(i=1;i<=m;i++)t[0][k]=t[0][k]+h*f(a+(i-0.5)*h);t[0][k]=(t[0][k]+t[0][k-1])/2;for(j=1;j<=k;j++){ c=1;for(n=0;n<j;n++)c=c*4;t[j][k-j]=(c*t[j-1][k-j+1]-t[j-1][k-j])/(c-1);}if(fabs(t[k][0]-t[k-1][0])<error){ printf("\n积分结果I ≈%lf\n",t[k][0]);break;}else{ h=h/2;k++;}}}#include <stdio.h>#include <math.h>double f(double t){double ans;ans=pow(cos(t),1.0/3);return ans;}void main(){double x=0,eslong=0.000001,x0;int N=20,i;printf("\n近似初值x0 = %lf\n",x);for(i=0;i<N;i++){x0=x;x=f(x);printf(" x%d = %lf\n",i+1,x);if(fabs(x-x0)<eslong)break;}if(fabs(x-x0)<eslong)printf("得到近似结果为x ≈%lf\n\n",x,i);elseprintf("迭代失败\n");}#include <stdio.h>#include <math.h>double a=0,b=1,x,y=0,h=0.1,k1,k2,k3,k4;int i,N;double f(double t,double s){double ans;ans=1+t*t;return ans;}void main(){N=(b-a)/h;x=a;printf("\n 初值为(x0,y0) = ( %.8f , %.8f )\n",x,y);for(i=0;i<N;i++){k1=f(x,y);k2=f(x+h/2,y+h*k1/2);k3=f(x+h/2,y+h*k2/2);k4=f(x+h,y+h*k3);y=y+h*(k1+2*(k2+k3)+k4)/6;x=x+h;printf(" 第%d次输出结果为(x%d,y%d) = ( %.8f , %.8f )\n",i+1,i+1,i+1,x,y);}}#include <stdio.h>void main(){double datax[4]={1.2,2.9,4.6,5.8},datay[10]={14.84,33.71,58.36,79.24},l[3],x=1.5,y;int i,j;y=0;for(i=0;i<=3;i++){l[i]=1;for(j=0;j<i;j++)l[i]=(x-datax[j])/(datax[i]-datax[j])*l[i];for(j=i+1;j<=3;j++)l[i]=(x-datax[j])/(datax[i]-datax[j])*l[i];y=y+datay[i]*l[i];}printf("\n f(x)在x = %f 处的近似值为: y = %f\n",x,y);}#include <stdio.h>void main(){double datay[9]={11.7,14.87,21.44,31.39,44.73,61.46,81.57,105.11,131.91};int m=2,i,j,k;double p,data[9][4],a[3][4],datax[9]={1.2,2.3,3.4,4.5,5.6,6.7,7.8,8.9,10.0};for(i=0;i<9;i++)for(j=1;j<2*m+1;j++){data[i][j]=1;for(k=0;k<j;k++)data[i][j]=data[i][j]*datax[i];}for(i=0;i<m+1;i++){for(j=0;j<m+1;j++){a[i][j]=0;for(k=0;k<9;k++)a[i][j]=a[i][j]+data[k][i+j];}}a[0][0]=9;a[0][m+1]=0;for(i=0;i<9;i++)a[0][m+1]=a[0][m+1]+datay[i];for(i=1;i<m+1;i++){ a[i][m+1]=0;for(j=0;j<9;j++){p=datay[j];for(k=0;k<i;k++)p=p*datax[j];a[i][m+1]=a[i][m+1]+p;}} //生成m+1行，m+2列增广矩阵// for(i=0;i<m+1;i++) //显示方程组//for(j=0;j<m+2;j++){if(j!=m+1){ printf("(%f)a%d ",a[i][j],j);if(j!=m)printf("+ ");}elseprintf("= %f \n",a[i][j]);}for(i=0;i<m;i++) //高斯消去法//{if(a[i][i]!=0){ for(j=i+1;j<m+1;j++){ a[j][i]=a[j][i]/a[i][i];for(k=i+1;k<m+2;k++)a[j][k]=a[j][k]-a[i][k]*a[j][i];}}elsebreak;}if(a[m][m]!=0&&i==m){ a[m][m+1]=a[m][m+1]/a[m][m];for(i=2;i<=m+1;i++){ for(j=1;j<i;j++)a[m+1-i][m+1]=a[m+1-i][m+1]-a[m+1-i][m+1-j]*a[m+1-j][m+1];a[m+1-i][m+1]=a[m+1-i][m+1]/a[m+1-i][m+1-i];}printf("方程组的解为:\n");for(j=0;j<m+1;j++)printf("a%d = %f\n",j,a[j][m+1]);printf("拟合多项式为:\n");printf("P%d(x) = (%f) + (%f)x + (%f)x^2\n",m,a[0][m+1],a[1][m+1],a[2][m+1]);}elseprintf("数据有误！\n");}}列主元素法#include <stdio.h>#include <math.h>void main(){double a[3][4]={1,-2,-1,3,-2,10,-3,15,-1,-2,5,10},mov,comp;int i,j,k,nrow;for(i=0;i<2;i++){comp=fabs(a[i][i]);for(k=i;k<3;k++) //比较绝对值大小并进行主元列交换// if(fabs(a[k][i])>=comp){nrow=k;comp=fabs(a[k][i]);}for(j=0;j<=3;j++){mov=a[i][j];a[i][j]=a[nrow][j];a[nrow][j]=mov;}printf("方程第%d行互换位置后如下\n",i+1);for(j=0;j<3;j++)printf("(%f)x1 + (%f)x2 + (%f)x3 = %f\n",a[j][0],a[j][1],a[j][2],a[j][3]);if(a[i][i]!=0){for(j=i+1;j<3;j++){a[j][i]=a[j][i]/a[i][i];for(k=i+1;k<=3;k++)a[j][k]=a[j][k]-a[i][k]*a[j][i];a[j][i]=0;}printf("方程经%d次消元如下\n",i+1);for(j=0;j<3;j++)printf("(%f)x1 + (%f)x2 + (%f)x3 = %f\n",a[j][0],a[j][1],a[j][2],a[j][3]);}elsebreak;}if(a[2][2]!=0&&i==2){printf("方程化简得\n");for(i=0;i<3;i++)printf("(%f)x1 + (%f)x2 + (%f)x3 = %f\n",a[i][0],a[i][1],a[i][2],a[i][3]);a[2][3]=a[2][3]/a[2][2];for(i=2;i<=3;i++){for(j=1;j<i;j++)a[3-i][3]=a[3-i][3]-a[3-i][3-j]*a[3-j][3];a[3-i][3]=a[3-i][3]/a[3-i][3-i];}printf("方程组的解为:\n");for(j=0;j<3;j++)printf("x%d = %f\n",j+1,a[j][3]);}elseprintf("数据有误！\n");}Jacobi迭代法#include <stdio.h>#include <math.h>void main(){double a[3][7]={{1,-2,-1,3},{-2,10,-3,15},{-1,-2,5,10}},error=0.000001,norm;int N=423,i,j,k;a[0][4]=0,a[1][4]=0,a[2][4]=0;for(i=0;i<3;i++) //把a矩阵转化为b矩阵//{a[i][6]=a[i][i];for(j=0;j<3;j++){a[i][j]=-a[i][j]/a[i][6];}a[i][3]=a[i][3]/a[i][6];a[i][i]=0;}printf("化为b矩阵如下\n");for(i=0;i<3;i++){printf("%f %f %f %f\n",a[i][0],a[i][1],a[i][2],a[i][3]);}for(i=1;i<N;i++){for(j=0;j<3;j++){a[j][5]=0;for(k=0;k<3;k++){a[j][5]=a[k][4]*a[j][k]+a[j][5];}a[j][5]=a[j][5]+a[j][3];}norm=0;for(k=0;k<3;k++)norm=norm+fabs(a[k][4]-a[k][5]);if(norm<error)break;elsefor(k=0;k<3;k++)a[k][4]=a[k][5];}if(norm<error){printf("计算结果为\n");for(i=0;i<3;i++){printf(" x%d = %.3f\n",i+1,a[i][5]);}}elseprintf("迭代失败\n");}题目1#include "stdio.h"#include "math.h"double f(double x){double ans;ans=exp(x);return(ans);}void main(){double a=-1,b=1,error=0.0001,m=1,h,T0,T,F;int k;h=(b-a)/2;T0=h*(f(a)+f(b));while(1){F=0;for(k=1;k<=pow(2.0,m-1);k++)F=F+f(a+(2*k-1)*h);T=T0/2+h*F;if(fabs(T-T0)<error)break;m++;h=h/2;T0=T;}printf("积分结果为I ≈%f\n",T);}#include "stdio.h"double f(double t,double s){double ans;ans=1+t*t;return(ans);}void main(){double a=0,b=1,h=0.2,x0=0,y0=0,x,k1,k2,k3,y;int N,n;N=(b-a)/h;for(n=1;n<=N;n++){x=x0+h;k1=f(x0,y0);k2=f(x0+h/2,y0+h/2*k1);k3=f(x0+h,y0-h*k1+2*h*k2);y=y0+h/6*(k1+4*k2+k3);printf("第%d次输出结果为(%.8f,%.8f)\n",n,x,y);x0=x;y0=y;}}。
北航研究生数值分析编程大作业1
数值分析大作业一、算法设计方案1、矩阵初始化矩阵[]501501⨯=ij a A 的下半带宽r=2,上半带宽s=2，设置矩阵[][]5011++s r C ，在矩阵C 中检索矩阵A 中的带内元素ij a 的方法是：j s j i ij c a ,1++-=。
这样所需要的存储单元数大大减少，从而极大提高了运算效率。
2、利用幂法求出5011λλ，幂法迭代格式：0111111nk k k k kk T k k k u R y u u Ay y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止。
首先对于矩阵A 利用幂法迭代求出一个λ，然后求出矩阵B ，其中I A B λ-=(I 为单位矩阵），对矩阵B 进行幂法迭代，求出λ'，之后令λλλ+'=''，比较的大小与λλ''，大者为501λ，小者为1λ。
3、利用反幂法求出ik s λλ,反幂法迭代格式：0111111nk k k k kk T k k k u R y u Au y y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止，1s k λβ=。
每迭代一次都要求解一次线性方程组1-=k k y Au ，求解过程为：（1）作分解LU A =对于n k ,...,2,1=执行[][]s k n r k k k i c c c c c n s k k k j c cc c k s ks k t k s k r i t t s t i k s k i k s k i js j t k s j r k t t s t k j s j k j s j k <+++=-=++=-=+++----=++-++-++-++----=++-++-++-∑∑);,min(,...,2,1/)(:),min(,...,1,:,1,11),,1max(,1,1,1,11),,1max(,1,1,1（2）求解y Ux b Ly ==,（数组b 先是存放原方程组右端向量，后来存放中间向量y))1,...,2,1(/)(:/:),...,3,2(:,1),min(1.1.11),1max(,1--=-===-=+++-++-+--=++-∑∑n n i c x c b x c b x n i b c b b i s t n s i i t t s t i i i ns n n ti r i t t s t i i i使用反幂法，直接可以求得矩阵按模最小的特征值s λ。
数值分析编程题c语言.doc
数值分析编程题c语言.word 资料上机实习题一一、题目：已知A与b12.38412,2.115237,-一、题目：已知A与b12.38412,2.115237,：在高斯方法已求出x(m),x(m-1)的基础上，组合新的序列，从而加快收敛速度。
其算式：其中ω是超松弛因子，当ω1时，可以加快收敛速度3 、用消去法求解用追赶消去法求Bx=b的方法：, , , -其算式：其中ω是超松弛因子，当ω1时，可以加快收敛速度3 、用消去法求解用追赶消去法求Bx=b的方法：, , , ：y1(0.025000)=0.025000,y2(0.025000)=0.151599,y3(0.025000)=8.33542 6y1(0.045000)=0.045000,y2(0.045000)=0.312860,y3(0.045000)=7.5362 80y1(0.085000)=0.085000,y2(0.085000)=0.560581,y3(0.085000)=4.946 353y1(0.100000)=0.100000,y2(0.100000)=0.628881,y3(0.100000)=4.18 0993四、问题讨论：此程序用高阶方程组的四阶古典的Runge-Kutta方法求解初值问题，控制循环输出语句时，应用时不能用t=0.025（其他判断语句类似），因t是双精度格式的，只能控制t在某一个很小的范围内。
Runge_Kutta方法的优点:1. 精度高,不必用别的方法求开始几点的函数值.2. 可根据f’(t,y)变化的情况与需要的精度自动修改步长.3. 方法稳定.4. 程序简单,存储量少.Runge_Kutta方法的缺点:每步要计算函数值f(t,y)四次,在f(t,y)，较复杂时,工作量大, 且每一步缺乏可靠的检查.习题六用对分法求B(即A)的小于23且最接近23的特征值的近似值（误差小于0.1），再用反幂法求B的该特征值的更精确近似值及相应的特征向量。
(完整word版)matlab数值分析例题
1、 在MATLAB 中用Jacobi 迭代法讨论线性方程组,1231231234748212515x x x x x x x x x -+=⎧⎪-+=-⎨⎪-++=⎩（1）给出Jacobi 迭代法的迭代方程，并判定Jacobi 迭代法求解此方程组是否收敛。
（2）若收敛，编程求解该线性方程组.解（1）：A=[4 -1 1;4 —8 1;-2 1 5］ %线性方程组系数矩阵A =4 -1 1 4 -8 1 —2 1 5>> D=diag(diag(A)）D =4 0 0 0 —8 0 0 0 5〉〉 L=—tril （A,-1) ％ A 的下三角矩阵L =0 0 0 —4 0 0 2 —1 0〉〉U=-triu(A，1）% A的上三角矩阵U =0 1 —10 0 —10 0 0B=inv（D）*（L+U）% B为雅可比迭代矩阵B =0 0.2500 —0。
25000.5000 0 0.12500。
4000 —0.2000 0〉〉r=eigs(B,1）%B的谱半径r =0。
3347 〈1Jacobi迭代法收敛。
(2）在matlab上编写程序如下:A=[4 —1 1;4 -8 1;—2 1 5]；〉〉b=［7 —21 15］'；>〉x0=［0 0 0]’;〉〉［x,k］=jacobi（A,b，x0，1e—7)x =2。
00004.00003。
0000k =17附jacobi迭代法的matlab程序如下：function [x,k］=jacobi（A,b，x0，eps)% 采用Jacobi迭代法求Ax=b的解％A为系数矩阵％b为常数向量％x0为迭代初始向量％eps为解的精度控制max1= 300; %默认最多迭代300，超过300次给出警告D=diag（diag（A）)；％求A的对角矩阵L=-tril(A，—1); %求A的下三角阵U=—triu（A,1); %求A的上三角阵B=D\(L+U);f=D\b;x=B＊x0+f；k=1；%迭代次数while norm（x-x0）>=epsx0=x;x=B＊x0+f；k=k+1;if（k〉=max1)disp(’迭代超过300次，方程组可能不收敛’）；return；endend2、设有某实验数据如下：（1)在MATLAB中作图观察离散点的结构,用多项式拟合的方法拟合一个合适的多项式函数；（2）在MATLAB中作出离散点和拟合曲线图。
数值分析上机题参考答案.docx
如有帮助欢迎下载支持数值分析上机题姓名：陈作添学号： 040816习题 120．（上机题）舍入误差与有效数N11 3 1 1设S N，其精确值为 。
22 2 NN 1j 2j1（1）编制按从大到小的顺序111 ，计算 S 的通用程序。
S N1 321N 21 N22（2）编制按从小到大的顺序111，计算 S 的通用程序。
S N1(N 1)2 122 1NN 2（3）按两种顺序分别计算S 102 ， S 104 ， S 106 ，并指出有效位数。
（编制程序时用单精度）（4）通过本上机题，你明白了什么？按从大到小的顺序计算 S N 的通用程序为： 按从小到大的顺序计算 S N 的通用程序为：#include<iostream.h> #include<iostream.h> float sum(float N) float sum(float N) {{float j,s,sum=0; float j,s,sum=0; for(j=2;j<=N;j++) for(j=N;j>=2;j--) {{s=1/(j*j-1); s=1/(j*j-1); sum+=s;sum+=s;}}return sum;return sum;}}从大到小的顺序的值从小到大的顺序的值精确值有效位数从大到小从小到大0.7400490.740050.74004965 S 1020.7498520.74990.74994 4S 1040.7498520.7499990.74999936S 106通过本上机题， 看出按两种不同的顺序计算的结果是不相同的，按从大到小的顺序计算的值与精确值有较大的误差， 而按从小到大的顺序计算的值与精确值吻合。
从大到小的顺序计算得到的结果的有效位数少。
计算机在进行数值计算时会出现“大数吃小数”的现象，导致计算结果的精度有所降低，我们在计算机中进行同号数的加法时，采用绝对值较小者先加的算法，其结果的相对误差较小。
数值分析编程实例-1
数值计算方法上机题目1、利用Steffensen 方法和Muller (抛物线)方法计算下列方程2260()cos()x x f x e x -=++-=在区间13[,]上的实根，要求精度为510ε-=，并比较迭代次数。
Steffensen 法:function A=steffensen(f,x,ep) %f 为迭代函数 %x 为初始值 %ep 为精度 A=[];A(1,1)=x; for i=1:1000y=feval(f,A(i,1)); z=feval(f,y);A(i,2)=A(i,1)-(y-A(i,1)).^2/(z-2*y+A(i,1)); if abs(A(i,2)-A(i,1))<ep break endA(i+1,1)=A(i,2); end 执行:f=inline('log(6-2*cos(x)-2^(-x))'); y=1;ep=10.^-5;steffensen(f,y,ep) 结果: ans =1.0000 1.8682 1.8682 1.8295 1.8295 1.8294 1.8294 1.8294 Muller 法:function[A]=muller(f,x0,x1,x2,ep) %f 为原函数%x0为输入初值最小的那个 %x1为输入初值最大的那个 %x2为输入初值中间的那个 %ep 为输入精度 A=[];A(1,1)=x0; A(1,2)=x1; A(1,3)=x2; for i=1:1000y=(A(i,3)-A(i,2))/(A(i,2)-A(i,1)); x=1+y;a=feval(f,A(i,1))*y.^2-feval(f,A(i,2))*y*x+feval(f,A(i,3))*y; b=feval(f,A(i,1))*y.^2-feval(f,A(i,2))*x.^2+feval(f,A(i,3))*(y+x); c=feval(f,A(i,3))*x;h=-2*c/(b+sign(b)*sqrt(b.^2-4*a*c)); A(i,4)=A(i,3)+h*(A(i,3)-A(i,2)); if(h*(A(i,3)-A(i,2))<ep) break endA(i+1,1)=A(i,2); A(i+1,2)=A(i,3); A(i+1,3)=A(i,4); end 执行: x0=3; x1=1; x2=1.8; ep=10.^-5;f=inline('2^(-x)+exp(x)+2*cos(x)-6'); muller(f,x0,x1,x2,ep) 结果: ans =3.0000 1.0000 1.8000 1.8218 1.0000 1.8000 1.8218 1.8294 1.8000 1.8218 1.8294 1.82942、用Gauss 列主元消去法求解下列方程组123454062120130266319263600255222651621x x x x x ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦function X=gauss(A,B)%A为系数矩阵%B为值矩阵C=[A B];n=length(B);for i=1:n[y,c]=max(C(i:n,i));D=C(i,:);C(i,:)=C(c+i-1,:);C(c+i-1,:)=D;if C(i,i)==0breakelseC(i+1:n,i)=C(i+1:n,i)/C(i,i);C(i+1:n,i+1:n+1)=C(i+1:n,i+1:n+1)-C(i+1:n,i)*C(i,i+1:n+1); endendC(n,n+1)=C(n,n+1)/C(n,n);for i=n-1:-1:1C(1:i,n+1)=C(1:i,n+1)-C(1:i,i+1)*C(i+1,n+1);C(i,n+1)=C(i,n+1)/C(i,i);endX=C(1:n,n+1);执行:A=[4 0 6 0 2;0 1 3 0 2;6 3 19 2 6;0 0 2 5 -5;2 2 6 -5 16];B=[12 6 36 2 21]';gauss(A,B)结果:ans =1.00001.00001.00001.00001.00003、已知函数表f的近似值。
同济大学数值分析matlab编程题汇编.doc
同济大学数值分析matlab编程题汇编.MATLAB编程题库1.下面的数据表近似地满足函数,请适当变换成为线性最小二乘问题,编程求最好的系数,并在同一个图上画出所有数据和函数图像。
解:x=[-x=[:文件一文件二函数v=f(x)v=x . * log(x)-1；函数z=g(y)z=y. y-1；解以下内容:编辑gex AFA。
m函数[x ITER]=gex AFA(f，x0，x1，tol)ITER=0；而(标准(x1-编辑gex AFA。
m函数[x ITER )=gex AFA(f，x0，x1，tol)ITER=0；同时（标准(x1:解以下内容:编辑gsdiedai。
m函数[x iter]=gsdiedai(A，x0，b，tol)D=diag(diag(A))；函数[x iter]=gsdiedai(A，x0，b，tol)D=diag(diag(A))；L=D:编辑ksf 2。
m函数v=ksf2(x，y)v=y exp(x)x . * y；a=1；b=2；h=0.01n=(b-a)./h .x=[1:0.01:2]；y(1)-省略部分-0.5000 1.0000 0.5000-1.0000 1.0000 UU=2.0000 3.0000 4.0000 0-0.5000 7.0000 0 0-1.0000 x=林空间(0，1，11)；x ' ans=0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.90001.0000 x=[1 2 3 4]；y=[6 11 18 27]；p=polyfit(x，y，2)p=1.0000 2.0000 3.0000 diag(1(4，1)，1)ans=0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 012 .编程实现求解满足下列条件的区间[-1，2]上的三次样条函数S(x)，并画出此样条函数的图形: Xi-1 0 1 2 f(Xi)-1 0 1 0f(Xi)' 0-1函数splx=[-1 0 1 2]y=[0-1 0 1 0-1]PP=csape(x，y，'完整')[断点coefs、npolys、ncoefs、dim]=unmkpp(PP)xh=-1:0.133602 if-1=xh=0 yh=coefs(1，1)*(xh 1).3系数(1，2)*(xh 1).2系数(1，3)*(xh 1)系数(1，4)否则，如果0=xh=1 yh=系数(2，1)*(xh).3系数(2，2)*(xh).2系数(2，3)*(xh)系数(2，4)否则，如果1=xh=2 yh=系数(3，1)*(xh-1).3系数(3，2)*(xh-1).2系数(3，3)*(xh-1)系数(3，4)否则返回图(xh，yh，' r ')13 .二分法程序如果nargintol x=(a b)/2 fx=feval(f，x)如果符号(外汇)==符号(fa) a=x fa=fx elseif符号(外汇)==符号b=x FB=FX否则返回结束。
数值分析作业
1. 迭代法、初始值与收敛性(210x x --=)A 、迭代式为121-=+nn x x 时，程序为： clc;clear;x(1)=1.5;n=30;for i=2:30if i<=nx(i)=x(i-1)^2-1;endend结果为：B 、迭代式为n n x x /111+=+时，程序为：clc;clear;x(1)=1.5;n=30;for i=2:30if i<=nx(i)=x(i-1)^2-1;endend结果为：C 、迭代式为11+=+n n x x 时，程序为：clc;clear;x(1)=1.5;n=30;for i=2:30if i<=nx(i)=(x(i-1)+1)^0.5;endend结果为：由上述程序迭代结果可以看出，第一种方法不收敛，第二种和第三种迭代方法的收敛性较好，且第三种的收敛速度更快，收敛性更好。
第二种迭代方法与初始值的选取相关性不大，选取不同的初始值时，迭代次数相差不大。
第三种迭代方法与初始值选取相关性较大，当初始值较接近不动点时，迭代次数会减少。
2、分别用牛顿法、弦割法、抛物线法计算x3+4x2-10=0在[1，2]上的根。
（1）牛顿法编程调试，输入初始值x=1.5执行，计算显示的结果如下图所示：#include<stdio.h>#include<math.h>#include<stdlib.h>#include<time.h>#include<windows.h>void main(){double x,y,start,end,freq;LARGE_INTEGER t1,t2,tc;QueryPerformanceFrequency(&tc);int i=1;scanf("%lf",&x);QueryPerformanceCounter(&t1);y=x-(x*x*x+4*x*x-10)/(3*x*x+8*x);do{x=y;y=x-(x*x*x+4*x*x-10)/(3*x*x+8*x);printf("x=%f\n",x);i++;}while(fabs(x-y)>=0.0000001);printf("x=%f i=%d\n",x,i);QueryPerformanceCounter(&t2);freq=(double)tc.QuadPart;start=(double)t1.QuadPart;end=(double)t2.QuadPart;printf("cpu time: %f ",(end-start)/freq);system("pause");}（2）弦割法编程调试，输入初始值x1=1，x2=2执行，计算显示的结果如下图所示：#include<stdio.h>#include<math.h>#include<stdlib.h>#include<time.h>#include<windows.h>double f(double x){return(x+4)*x*x-10;}void main(){double x,x1,x2,start,end,freq;LARGE_INTEGER t1,t2,tc;QueryPerformanceFrequency(&tc);int i=1;scanf("%lf,%lf",&x1,&x2);QueryPerformanceCounter(&t1);x=(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1));do{x1=x;x=(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1));printf("x=%f\n",x);i++;}while(fabs(x-x1)>=0.0000001);printf("x=%f i=%d\n",x,i);QueryPerformanceCounter(&t2);freq=(double)tc.QuadPart;start=(double)t1.QuadPart;end=(double)t2.QuadPart;printf("cpu time: %f ",(end-start)/freq);system("pause");}（3）抛物线法编程调试，输入初始值x=1，y=1.5，z=2执行，计算显示的结果如下图所示：#include<stdio.h>#include<math.h>#include<stdlib.h>#include<time.h>#include<windows.h>double fx(double x){return x*x*x+4*x*x-10 ;}void main(){double x,y,z,a,b,c,l3,d3,l4,w,start,end,freq;LARGE_INTEGER t1,t2,tc;int i=1;QueryPerformanceFrequency(&tc);scanf("%lf,%lf,%lf",&x,&y,&z);QueryPerformanceCounter(&t1);loop:l3=(z-y)/(y-x);d3=1+l3;a=fx(x)*l3*l3-fx(y)*l3*d3+fx(z)*l3;b=fx(x)*l3*l3-fx(y)*d3*d3+fx(z)*(l3+d3);c=fx(z)*d3;l4=-2*c/(b+(b/fabs(b))*pow(b*b-4*a*c,0.5));w=z+l4*(z-y);printf("%d %f\n",i,w);if(fabs(l4*(w-z))>=0.0001){x=y;y=z;z=w;i++;goto loop;}QueryPerformanceCounter(&t2);freq=(double)tc.QuadPart;start=(double)t1.QuadPart;end=(double)t2.QuadPart;printf("cpu time: %f ",(end-start)*1000/freq);system("pause");}3、分别用牛顿法、弦割法、抛物线法计算e x+10x-2=0在(1，2)上的根。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《数值分析》实验报告
第二章=解线性方程组的直接方法
2、试用MATLAB软件编程实现追赶法求解三对角方程组的算法,
并考虑梯形电阻电路问题，电路如下：
其中电路中的各个电流f2?…，珀｝须满足下列线性方程组:
一耳十务一巧=0
一巧4■笃一為=0
一迢+呂一逼=0
一耳+兀_%=0
一爲+5% 一巧=0
-2f6+5f7-2f8 =0
-2i7+5r g=0
设F = 220V,氏= 27Q,运用求各段电路的电流量。
解：v!R = — «8_1481
27
上述方程组可用矩阵表示为：
・2_2000000h81481
_25_200000h0
0_25-20000h0
00_25_2000h0
000-25_200h0
0000_25-20S0
00000-25_2h0
000000-25?8.0
認翟蚁蛊
X
-p u a) A .H )q /(i p n .H )p T (.rH )p
T
二—二丄岂4
晒口®〈狠* A8)q
、
(8)PH(8)p
-p u a)
-U
—.H )
P
*-H )
F
(.H )
P A
.H )
P
丄)
y (.H )q A .H )q
Q
丄)q/u)y(2
晒口§前厂<
00
勺：。
4
坦祖咽鱼磋£紙夜
lll <
=0 0 0 0 0
0 0
T
畧
oo v
p 二甲
d — d —
甲甲甲
d —
v 。
二g
g g g
g g
g d v q
氏
—甲
甲甲甲甲
d — O V B
曲恥 qedelN
寸 61—<10905
30p ^
s
£0-
S S I
黑理q。
脅
u)
旺=味蛋二
迭代初始向量取评=(0, 0, 0, 0, 0)T o
解：实验步骤及程序、结果
取要求达到的精度f=i(r8。
以下程序中的丘均表示迭代次数。
(1) Jacobi迭代法MatLab源程序。
format long
4[10丄2,3,4;1,9厂1,2厂3;2厂1,7,3厂5;3,2,3,12厂1;4,-3厂
5厂1,15]; b=[12r27,14,-17,12];
x0=[0,0,0,0,0];xl=x0;
Nmax=1000;
k=0;
for i=l:5
sum=0;
for j=l:5
sum=sum+A(i,j) *x0 (j);
end;
end;
Nl(i)=(b(i)
end;
while abs(norm(x 1 -x0,inf))>le-8 &k<Nmax x0=xl;
for i=l:5
sum=0;
for j=l:5
sum=sum+A(i,j) *x0 (j);
end;
end; xl(i)=(b 诡
end;
k=k+l;
end;
xl
k
输岀结果=
xl =
1.0000000 -
2.0000000 2.8446139 -1.3335823 0.
k=88
迭代了88次。
(2)Gauss-Seidel 迭代法MatLab 源程序。
4[10丄2,3,4;1,9厂1,2,-3;2厂1,7,3厂5;323,12厂1;4,-3厂5厂
1,1习;
b=[12,-27,14,-17,12J;
xO=[O,O,O,O,O];xl=xO;
Nmax=1000;
k=l;sum=O;
for j=2:5
sum=sum+A( 1 ,j) *x0 (j);
end;
for i=2:4
sum=0;
forj=l:(i-l)
sum=sum+A(i,j) *x0 (j);
end;
xl(i)=b
sum=0;
forj=(i+l): 了
sum=sum+A(i,j) *x0 (j);
end;
xl(i)=(x 1 ①
end;
sum=0;
forj=l:4
sum=sum+A(5 ,j) *x0 (j); end;
xl(5)=(b(5)-sum)/A(5,5);
while abs(norm(x 1 -xO,inf))>le-8 & k<Nmax xO=xl;
sum=0;
for j=2:5
sum=sum+A( 1 ,j) *x0 (j);
end;
for i=2:4
sum=0;
for 尸1©1)
sum=sum+A(i,j) *x0 (j);
end;
xl(i)=b
sum=0;
for 尸(i+l):了
sum=sum+A(i,j) *x0 ①；
end;
xl(i)=(x 1 ①-sum)/A(i,i)；
end;。