数值分析编程 - 360文档中心

数值分析课后复习(C语言学习知识编程实现)

#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;2;t[0][0]=h*(f(a)+f(b))/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);2;t[0][k]=(t[0][k]+t[0][k-1])/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)得到近似结果为 x ≈ %lf\n\n",x,i);printf("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;初值为 (x0,y0) = ( %.8f , %.8f )\n",x,y);printf("\nfor(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);6;y=y+h*(k1+2*(k2+k3)+k4)/x=x+h;第%d次输出结果为 (x%d,y%d) = ( %.8f , %.8f )\n",i+1,i+1,i+1,x,y);printf("}}#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];}在 x = %f 处的近似值为: y = %f\n",x,y);printf("\n f(x)}#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];}方程组的解为:\n");printf("for(j=0;j<m+1;j++)printf("a%d = %f\n",j,a[j][m+1]);拟合多项式为:\n");printf("printf("P%d(x) = (%f) + (%f)x + (%f)x^2\n",m,a[0][m+1],a[1][m+1],a[2][m+1]);}else数据有误！\n");printf("}}列主元素法#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;}方程第%d行互换位置后如下\n",i+1);printf("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;}方程经%d次消元如下\n",i+1);printf("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){方程化简得\n");printf("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];}方程组的解为:\n");printf("for(j=0;j<3;j++)printf("x%d = %f\n",j+1,a[j][3]);}elseprintf("数据有误！\n");}#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;//把a矩阵转化为b矩阵//for(i=0;i<3;i++){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;}化为b矩阵如下\n");printf("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){计算结果为\n");printf("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);第%d次输出结果为(%.8f,%.8f)\n",n,x,y);printf("x0=x;y0=y;}}。

数值分析作业(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的值，本题为了编程方便，已经把数据写入程序中。

(Fortran编程)数值分析之Newton法求解非线性方程组

BUP(:)=AB(:,N+1)
!调用上三角方程组的回代方法
call UP_TRI(AUP,BUP,X,N)
end subroutine LINEQ
subroutine UP_TRI(A,B,X,N)
!上三角方程组的回代方法
imp阵
open(unit=11,file="jieguo.txt")
write(11,100)
100 format(/,T6,"Newton法计算非线性方程组迭代序列：",/)
X=(/1.6,1.2/)
TOL=1D-9
do I=1,ITMAX
ELMAX=DABS(AB(K,K))
ID_MAX=K
!这段程序不是为了赋值最大元素给ELMAX,而是为了找出最大元素对应的标号
do I=K+1,N
if(dabs(AB(I,K))>ELMAX) then
ELMAX=AB(I,K)
module M_GAUSS
!高斯列主元消去法模块
contains
subroutine LINEQ(A,B,X,N)
!高斯列主元消去法
implicit real*8(A-Z)
integer::I,K,N
integer::ID_MAX !主元素标号
real*8::A(N,N),B(N),X(N)
!调用Newton法中的函数
call SOLVE()
end program newton
end do
X(I)=X(I)/A(I,I)
end do
end subroutine UP_TRI

混凝土结构M-N关系图计算机数值分析C++编程

练习2：一、题目钢筋混凝土矩形截面：b=300mm ，h=600mm ，h 0=560mm ，a s ’=25mm ，a s =40mm ，A s ’=157mm2，A s =804mm2，f y ’=280MPa ，f y =280MPa ，E s =200GPa ，E c =25.5GPa ，f c =13.4MPa ，f t =1.54MPa ，ε0=0.002，εcu =0.0038，εs u ≤10%=0.10。

.利用数值方法计算截面的M~N 关系，并附简化计算结果N u 。

b=300mm2Φ10h=600mm4Φ16二、简单分析：本次作业是在上一次作业的基础上继续进行其他计算。

主要任务是利用计算机软件来计算特定截面偏心受压情况下的纵向压力与弯矩的关系。

可参考第一次作业的程序，做适当的修改即可。

混凝土应力—应变曲线采用的是R üsch 建议的曲线。

曲线的上升段采用抛物线形式，下降段为斜直线。

R üsch 建议的曲线：当0εε≤时，'2002[()]c cf εεσεε=-当0cu εεε<≤时，'c c f σ=根据《高等钢筋混凝土结构学》提供的公式：由平衡条件：0N =∑，0M =∑可知，''0()()d x c ci s s s sN b y y A A σεσσ=+-⎰0'''00000()()()(+)d ()x s c ci i s s s N e y b y h x y y A h a σεσ+=-+-⎰对于离纵向力较远钢筋应力的取值可参照以下情况：平衡破坏：即受压边缘混凝土应变恰好达到极限应变时，受拉钢筋刚好达到屈服强度280MPa 。

受拉破坏：荷载偏心较大时，钢筋先屈服（达到280MPa ），经过一个过程后，混凝土达到极限压应变。

受压破坏：荷载偏心较小时，构件产生受压破坏。

受压破坏是指受压较大一侧的混凝土达到极限压应变，而离纵向力较远一侧的钢筋可能受拉或者受压但都不屈服。

数值分析课程设计（最终版）

数值分析课程设计（最终版）本⽂主要通过Matlab 软件，对数值分析中的LU 分解法、最⼩⼆乘法、复化Simpon 积分、Runge-Kutta ⽅法进⾏编程，并利⽤这些⽅法在MATLAB 中对⼀些问题进⾏求解，并得出结论。

实验⼀线性⽅程组数值解法中，本⽂选取LU 分解法，并选取数据于《数值分析》教材第5章第153页例5进⾏实验。

所谓LU 分解法就是将⾼斯消去法改写为紧凑形式，可以直接从矩阵A 的元素得到计算L 、U 元素的递推公式，⽽不需要任何步骤。

⽤此⽅法得到L 、U 矩阵，从⽽计算Y 、X 。

实验⼆插值法和数据拟合中，本⽂选取最⼩⼆乘拟合⽅法进⾏实验，数据来源于我们课堂学习该章节时的课件中的多项式拟合例⼦进⾏实验。

最⼩⼆乘拟合是⼀种数学上的近似和优化，利⽤已知的数据得出⼀条直线或者曲线，使之在坐标系上与已知数据之间的距离的平⽅和最⼩。

利⽤excel 的⾃带函数可以较为⽅便的拟合线性的数据分析。

实验三数值积分中，本⽂选取复化Simpon 积分⽅法进⾏实验，通过将复化Simpson 公式编译成MATLAB 语⾔求积分∫e ;x dx 10完成实验过程的同时，也对复化Simpon 积分章节的知识进⾏了巩固。

实验四常微分⽅程数值解，本⽂选取Runge-Kutta ⽅法进⾏实验，通过实验了解Runge-Kutta 法的收敛性与稳定性同时学会了学会⽤Matlab 编程实现Runge-Kutta 法解常微分⽅程，并在实验的过程中意识到尽管我们熟知的四种⽅法，事实上，在求解微分⽅程初值问题，四阶法是单步长中最优秀的⽅法，通常都是⽤该⽅法求解的实际问题，计算效果⽐较理想的。

实验五数值⽅法实际应⽤，本⽂采⽤最⼩⼆乘法拟合我国2001年到2015年的⼈⼝增长模型，并预测2020年我国⼈⼝数量。

关键词：Matlab ；LU 分解法；最⼩⼆乘法；复化Simpon 积分；Runge-Kutta⼀．LU分解法 (1)1.1实验⽬的 (1)1.2基本原理 (1)1.3实验内容 (2)1.4数据来源 (3)1.5实验结论 (3)⼆．Lagrange插值 (4)2.1实验⽬的 (4)2.2基本原理 (5)2.3实验内容 (5)2.4数据来源 (6)2.5实验结论 (6)三．复化simpon积分 (7)3.1实验⽬的 (7)3.2基本原理 (7)3.3实验内容 (7)3.4数据来源 (8)3.5实验结论 (8)四．Runge-Kutta⽅法 (9)4.1实验⽬的 (9)4.2基本原理 (9)4.3实验内容 (10)4.4数据来源 (11)4.5实验结论 (11)五．数值⽅法实际应⽤ (11)5.1实验⽬的 (11)5.2基本原理 (12)5.3实验内容 (12)5.4数据来源 (13)5.5实验结论 (13)总结 (16)参考⽂献 (17)⼀．LU 分解法1.1实验⽬的[1] 了解LU 分解法的基本原理和⽅法；[2] 通过实例掌握⽤MATLAB 求线性⽅程组数值解的⽅法； [3] 编程实现LU 分解1.2基本原理对于矩阵A ，若存在⼀个单位下三⾓矩阵L 和⼀个上三⾓U ，使得A =LU （1.1）。

数值分析编程报告

数值分析编程报告思源1201班12274005刘家玮列主消元法解方程一、解题方法的理论依据列主消元法和高斯消元法接方程的基本思想是一致的，唯一不同的是：在把方程组的增广矩阵化成上三角矩阵的时候，高斯消元法用的是高斯消元法，即不交换仍和行列的秩序，一次一行一行消元。

而列主消元法是在高斯消元法的基础上，增加了一个步骤，即在消元前看主元在该列中的绝对值是否为最大，若最大，则进行消元，若不是，则进行行交换，保证每次消元时主元的绝对值在该行中最大。

此目的是为了控制消元时引起的误差。

计算公式：x i =1/a ii (i)[b i (i)-∑+=ni j 1j (i)ij x a ],i=n,n-1,…,1a 为方程的系数矩阵，b 为方程右边的矩阵公式中以n 阶为例。

二、程序设计思路及计算程序 2 设计思路2、计算程序 #include"stdio.h" #include"math.h" main() {int i,j,k,m,h,max;double x[9],l,x0=0,t=0，sum=0; 定义变量 doubleb[9]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392},A[9][9]={{12.38412 ,2.115237 ,-1.061704,1.112336 ,-0.113584,0.718719 ,1.742382 ,3.067813 , -2.031743},{2.115237 ,19.141823,-3.125432,-1.012345,2.189736 ,1.563849 ,-0.784165 ,1.112348 ,3.123124 },{-1.061074,-3.125432,15.567914,3.123848 ,2.031454 ,1.836742 ,-1.056781 ,0.336993 ,-1.010103 },{1.112336 ,-1.012345,3.123848 ,27.108437,4.101011 ,-3.741856,2.101023 ,-0.71828 ,-0.037585 },{-0.113584,2.189736 ,2.031454 ,4.010111 ,19.897918,0.431637 ,-3.111223 ,2.121314 ,1.784317 } ,{0.718719 ,1.563849 ,1.836742 ,-3.741856,0.431637 ,9.789365 ,-0.103458 ,-1.103456,0.238417 } ,{1.742382 ,-0.784165,-1.056781,2.101023 ,-3.111223,-0.103458,14.7138465,3.123789 ,-2.21347 4},{3.067813 ,1.112348 ,0.336993 ,-0.71828 ,2.121314 ,-1.103456,3.123789 ,30.719334,4.446782 },{-2.031743,3.123124 ,-1.010103,-0.037585,1.784317 ,0.238417 ,-2.213474 ,4.446782 ,40.00001 } };for(j=0;j<=7;j++){max=j;for(m=j;m<=7;m++) 找出最大绝对值行元素所在的行标{if(fabs(A[max][j])-fabs(A[m+1][j])>=0) h=max;else max=h=m+1;}for(m=j;m<=8;m++){t=A[j][m];sum=b[j];A[j][m]=A[h][m]; 行交换，保证消元时主元的绝对值在该列为最大A[h][m]=t;b[j]=b[h];b[h]=sum;}for(i=j+1;i<=8;i++){l=A[i][j]/A[j][j]; 算出消元系数for(k=j;k<=8;k++) 消元A[i][k]=A[i][k]-l*A[j][k];b[i]=b[i]-l*b[j];}}for(i=8;i>=0;i--) 求出方程组的解{for(j=i+1;j<=8;j++)x0=x0+A[i][j]*x[j];x[i]=1.0/A[i][i]*(b[i]-x0);x0=0.0;}printf("The result is:\n"); 打印出结果for(i=0;i<=8;i++)printf("%lf\n",x[i]);}三、计算结果列主消元法的结果松弛法迭代9次时的值松弛法迭代20次时的值四、问题讨论由于本程序循环中数组的下标对应比较复杂，在调试的过程中飞了不少劲，最后通过手动一层层的推导才得出了正确的对应关系的。

数值分析MATLAB编程——数值积分法

数值分析MATLAB编程——数值积分法1、调用函数--f.Mfunction y=f(x)%------------------------------------------------------------函数1 y=sqrt(4-sin(x)*sin(x));%------------------------------------------------------------函数2 %y=sin(x)/x;%if x==0% y=0;%end%------------------------------------------------------------函数3 %y=exp(x)/(4+x*x);%------------------------------------------------------------函数4 %y=(log(1+x))/(1+x*x);2、复合梯形公式--tixing.M%复合梯形公式clear alla=input('请输入积分下限:');b=input('请输入积分上限:');n=input('区间n等分：');h=(b-a)/n;x=a:h:b;T=0;for k=1:n;T=0.5*h*(f(x(k))+f(x(k+1)))+T;endT=vpa(T,8)3、复合Simpson公式--simpson.M%复合Simpson公式clear alla=input('请输入积分下限:');b=input('请输入积分上限:');n=input('区间n等分：');h=(b-a)/n;x=a:h:b;S=0;for k=1:n;xx=(x(k)+x(k+1))/2;S=(1/6)*h*(f(x(k))+4*f(xx)+f(x(k+1)))+S;endS=vpa(S,8)4、Romberg算法--romberg.M%Romberg算法clear alla=input('请输入积分下限:');b=input('请输入积分上限:');n=input('区间n等分：');num=0:n;R=[num'];h=b-a;T=h*(f(a)+f(b))/2;t(1)=T;for i=2:n+1;u=h/2;H=0;x=a+u;while x<b;H=H+f(x);x=x+h;endt(i)=(T+h*H)/2;T=t(i);h=u;endR=[R,t'];for i=2:n+1for j=n+1:-1:1if j>=it(j)=(4^(i-1)*t(j)-t(j-1))/(4^(i-1)-1);elset(j)=0;endendR=[R,t'];endR=vpa(R,8)R(n,n)5、变步长算法（以复化梯形公式为例）--tixing2.M%复合梯形公式，确定最佳步长format longclear alla=input('请输入积分下限:');b=input('请输入积分上限:');eps=input('请输入误差：');k=1;T1=(b-a)*(f(a)+f(b))/2;T2=(T1+(b-a)*(f((a+b)/2)))/2; while abs((T1-T2)/3)>=epsM=0;n=2^k;h=(b-a)/n;T1=T2;x=a:h:b;for i=1:n;xx=(x(i)+x(i+1))/2;M=M+f(xx);endT2=(T1+h*M)/2;k=k+1;endT=vpa(T2,8)n=2^k。

数值分析.pdf

1、证明1 - x - sin x = 0 在[0,1]内有一个根，使用二分法求误差不大于 21 ⨯10-4 的根要迭代多少次并输出每一步的迭代解和迭代误差。

解：因为()()1*2/2/+-=-≤-k k k k a b a b x x ， ()4110*212/01-+≤-k ，所以经过计算之后有14=k ，即需要迭代14次。

运用matlab 编程法程序如下：function T1a = 0;b = 1;fprintf('n || a || b || c \n') for k=1:14c = (a+b)/2;fa = 1-a-sin(a);fb = 1-b-sin(b);fc = 1-c-sin(c);fprintf('%d || %f || %f || %f \n',k,a,b,c);if abs(fc)<(1/2)*10^(-4) r=c; sprintf('the root is: %d' , r); elseif fa*fc<0 b=c;elseif fb*fc<0 a=c;endendroot = (a+b)/2运行结果如下：n || a || b || c1 || 0.000000 || 1.000000 || 0.5000002 || 0.500000 || 1.000000 || 0.7500003 || 0.500000 || 0.750000 || 0.6250004 || 0.500000 || 0.625000 || 0.5625005 || 0.500000 || 0.562500 || 0.5312506 || 0.500000 || 0.531250 || 0.5156257 || 0.500000 || 0.515625 || 0.5078138 || 0.507813 || 0.515625 || 0.5117199 || 0.507813 || 0.511719 || 0.50976610 || 0.509766 || 0.511719 || 0.51074211 || 0.510742 || 0.511719 || 0.51123012 || 0.510742 || 0.511230 || 0.51098613 || 0.510742 || 0.511230 || 0.51098614 || 0.510742 || 0.511230 || 0.510986root =0.5110二求解方程x e x -=的根，要求5.00=x ，分别用简单迭代法，迭代法的加速方法：)(1)(11k k k k x x p p x f x --==-++-，以及艾特金方法求解，要求误差应满足｜k k x x -+1｜510-≤ 解：简单迭代法：程序如下：x0 = 0.5;epsilon = 10^(-5); e = 1;fprintf('k || xk ||error \n'); k = 1;while e>epsilonx1 = exp(-x0);e = abs(x1-x0);x0=x1;k=k+1;fprintf('%d || %f ||%f\n',k,x1,e); endx0k || xk ||error2 || 0.606531 ||0.1065313 || 0.545239 ||0.0612914 || 0.579703 ||0.0344645 || 0.560065 ||0.0196386 || 0.571172 ||0.0111087 || 0.564863 ||0.0063098 || 0.568438 ||0.0035759 || 0.566409 ||0.00202910 || 0.567560 ||0.00115011 || 0.566907 ||0.00065212 || 0.567277 ||0.00037013 || 0.567067 ||0.00021014 || 0.567186 ||0.00011915 || 0.567119 ||0.00006716 || 0.567157 ||0.00003817 || 0.567135 ||0.00002218 || 0.567148 ||0.00001219 || 0.567141 ||0.000007x0 =0.5671迭代法的加速方法：程序如下：x0 = 0.5;epsilon = 10^(-5);e = 1;fprintf('k || xk ||error \n');k = 1;while e>epsilonx11=exp(-x0);x1=x11-p*(x11-x0)/(p-1);e = abs(x0-x1);x0=x1;k=k+1;fprintf('%d || %f ||%f\n',k,x1,e);endx0结果如下：k || xk ||error2 || 0.566582 ||0.0665823 || 0.567132 ||0.0005504 || 0.567143 ||0.0000115 || 0.567143 ||0.000000x0 = 0.5671艾特金加速迭代法：程序如下：x0 = 0.5;epsilon = 10^(-5);e = 1;fprintf('k || xk ||error \n');k = 1;while e>epsilonx1=exp(-x0);x2=exp(-x1);x00=x0-(x1-x0)^2/(x0-2*x1+x2);e = abs(x0-x00);x0=x00;k=k+1;fprintf('%d || %f ||%f\n',k,x1,e);endx0k || xk ||error2 || 0.606531 ||0.0676243 || 0.566871 ||0.0004814 || 0.567143 ||0.000000x0 =0.56713、分别用单点弦割法和双点弦割法求f(x)=x3+2x2+10x-20=0的根，要求x k+1- x k<10-6。

浅谈《数值分析》课程教学中的编程能力培养问题

楚雄师范学院学报
２０１３年第６期
３Ｏ
ｒ．１｝＝１：ｍａｘ１
４０
ｒＪ＝１：Ｎ
５０
６０
ｘ（ｊ）＝（Ｂ（ｊ）一ａ（ｊ，［１：一１ √＋１： Ⅳ］）Ｐ（［１：一１ √ ＋１： Ⅳ］））／ａ（ｊ，Ｊ）；
目标。
＝
＝
关键词：数值分析；编程能力；培养
中图分类号：０２４１文章标识码：Ａ文章编号：１６７１— ７４０６（２０１３）０６—００１１一Ｏ６
６１６２
《数值分析》课程是信息与计算科学专业中计算数学方向的核心课程，《信息与计算科学专业教学规范（讨论稿）》 … 中就指出该专业的业务培养要求之一为：具备熟练应用计算机（包括常
，
二
二：：：二＝
二二）
其中Ｊ＝１，２， …，Ｎ，ｋ＝１，２， …。教材［２］在第２１２页给出了雅可比迭代的程序（行号是为了叙述方便加上的，程序中没有）如
下：
１０ｆｕｎｃｔｉｏｎＸ＝ｊａｃｏｂｉ（Ａ，Ｂ，Ｐ，ｄｅｌｔａ，ｍａｘ１）
ｆ（／ｅｆｔ＂＜ｄｅｌｔａ）Ｉ（ＦｅｌｅＩ＇ｒ＂＜ｄｅｌｔａ）

北航研究生数值分析编程大作业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 λ。

东南大学数值分析编程作业

数值分析编程作业******学号：******指导老师：***目录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种计算结果小，这是由于大数吃小数的原因（计算机表示一个数值是用字节表示的，当程序进行计算时，先将第一项和第二项相加，然后再加第三项……所以加到后面会由于项的值很小，不能加到这个字节里，所以就被忽略了。

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

数值分析上机实验报告《数值分析》上机实验报告1.用Newton 法求方程 X 7-X 4+14=0在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。

1.1 理论依据：设函数在有限区间[a ，b]上二阶导数存在，且满足条件{}αϕ上的惟一解在区间平方收敛于方程所生的迭代序列迭代过程由则对任意初始近似值达到的一个中使是其中上不变号在区间],[0)(3,2,1,0,)(')()(],,[x |))(),((|,|,)(||)(|.4;0)(.3],[)(.20)()(.110......b a x f x k x f x f x x x Newton b a b f a f mir b a c x f ab c f x f b a x f b f x f k k k k k k ==-==∈≤-≠>+令)9.1()9.1(0)8(4233642)(0)16(71127)(0)9.1(,0)1.0(,1428)(3225333647>⋅''<-=-=''<-=-='<>+-=f f x x x x x f x x x x x f f f x x x f故以1.9为起点⎪⎩⎪⎨⎧='-=+9.1)()(01x x f x f x x k k k k 如此一次一次的迭代，逼近x 的真实根。

当前后两个的差<=ε时，就认为求出了近似的根。

本程序用Newton 法求代数方程(最高次数不大于10)在（a,b ）区间的根。

1.2 C语言程序原代码：#include<stdio.h>#include<math.h>main(){double x2,f,f1;double x1=1.9; //取初值为1.9do{x2=x1;f=pow(x2,7)-28*pow(x2,4)+14;f1=7*pow(x2,6)-4*28*pow(x2,3);x1=x2-f/f1;}while(fabs(x1-x2)>=0.00001||x1<0.1); //限制循环次数printf("计算结果：x=%f\n",x1);}1.3 运行结果：1.4 MATLAB上机程序function y=Newton(f,df,x0,eps,M)d=0;for k=1:Mif feval(df,x0)==0d=2;breakelsex1=x0-feval(f,x0)/feval(df,x0);ende=abs(x1-x0);x0=x1;if e<=eps&&abs(feval(f,x1))<=epsd=1;breakendendif d==1y=x1;elseif d==0y='迭代M次失败';elsey= '奇异'endfunction y=df(x)y=7*x^6-28*4*x^3;Endfunction y=f(x)y=x^7-28*x^4+14;End>> x0=1.9;>> eps=0.00001;>> M=100;>> x=Newton('f','df',x0,eps,M);>> vpa(x,7)1.5 问题讨论：1.使用此方法求方解，用误差来控制循环迭代次数，可以在误差允许的范围内得到比较理想的计算结果。

数值分析在生活中的应用举例及Matlab实现

一、最小二乘法，用MATLAB实现1. 数值实例下面给定的是乌鲁木齐最近1个月早晨7：00左右（新疆时间）的天气预报所得到的温度，按照数据找出任意次曲线拟合方程和它的图像。

下面用MATLAB编程对上述数据进行最小二乘拟合。

下面用MATLAB编程对上述数据进行最小二乘拟合2、程序代码x=[1:1:30];y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1];a1=polyfit(x,y,3) %三次多项式拟合%a2= polyfit(x,y,9) %九次多项式拟合%a3= polyfit(x,y,15) %十五次多项式拟合%b1=polyval(a1,x)b2=polyval(a2,x)b3=polyval(a3,x)r1= sum((y-b1).^2) %三次多项式误差平方和%r2= sum((y-b2).^2) %九次次多项式误差平方和%r3= sum((y-b3).^2) %十五次多项式误差平方和%plot(x,y,'*') %用*画出x,y图像%hold onplot(x,b1, 'r') %用红色线画出x,b1图像%hold onplot(x,b2, 'g') %用绿色线画出x,b2图像%hold onplot(x,b3, 'b:o') %用蓝色o线画出x,b3图像%3、数值结果不同次数多项式拟合误差平方和为：r1=67.6659r2=20.1060r3=3.7952r1、r2、r3分别表示三次、九次、十五次多项式误差平方和。

4、拟合曲线如下图二、线性方程组的求解（高斯-塞德尔迭代算法）1、实例：求解线性方程组（见书P233页）⎪⎪⎩⎪⎪⎨⎧=++=-+=+-3612363311420238321321321x x x x x x x x x 记A x=b, 其中⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=363320,,12361114238321b x A x x x任取初始值()()Tx0000=，进行迭代。

fortran语言的基本概念

fortran语言的基本概念Fortran语言的基本概念概述Fortran（Formula Translation）是一种用于科学计算和数值分析的高级编程语言。

它是最早的编程语言之一，由IBM公司在20世纪50年代开发。

Fortran语言以其高效、可靠和快速的特点，长期被广泛应用于科学和工程计算领域。

基本结构Fortran程序以程序单元作为基本的执行单位，程序单元分为主程序和子程序两种类型。

主程序•主程序是Fortran程序的入口点，每个Fortran程序必须包含一个主程序。

•主程序由PROGRAM关键字定义，后接程序名称。

•主程序内包含了程序的执行逻辑和控制流程，通过调用子程序完成具体的计算任务。

•子程序是主程序的辅助部分，用于定义重复使用的计算任务或功能模块。

•子程序由SUBROUTINE或FUNCTION关键字定义。

•SUBROUTINE用于定义过程（Procedure），只执行特定的计算任务，不返回结果。

•FUNCTION也是定义过程，但它还可以返回一个值作为结果。

变量与数据类型Fortran语言使用变量来存储和操作数据，变量需要先进行声明。

Fortran提供了多种数据类型，用于存储不同类型的数据。

声明变量•变量声明使用DIMENSION或INTEGER/REAL/COMPLEX/LOGICAL等关键字。

•变量声明语句通常出现在主程序的开头或子程序的参数列表中。

数据类型•INTEGER用于表示整数类型的数据。

•REAL用于表示浮点数类型的数据。

•COMPLEX用于表示复数类型的数据。

•LOGICAL用于表示逻辑类型的数据，可取值为TRUE或FALSE。

Fortran语言提供了多种控制流程语句，用于在程序中实现条件判断和循环操作。

条件判断•IF-THEN语句用于根据条件执行不同的代码块。

•IF-THEN-ELSE语句可根据条件选择执行不同的代码块。

循环•DO语句用于实现循环操作，可指定循环的起始值、结束值和步进值。

数值分析实验报告

《数值分析》实验报告学院：计算机科学与软件学院姓名：XXX班级：计算机XX班学号：XXXXXX实验一：舍入误差与数值稳定性实验目的：1、通过上机编程，复习巩固以前所学程序设计语言；2、通过上机计算，了解舍入误差所引起的数值不稳定性。

3、通过上机计算，了解运算次序对计算结果的影响，从而尽量避免大数吃小数的现象。

实验内容：用两种不同的顺序计算644834.11000012≈∑=-n n ，分析其误差的变化。

实验流程图：实验源程序:#include <stdio.h>#include <math.h>void main(){ int i;float s1=0,s2=0,d1,d2;for (i=1;i<=10000;i++)s1=s1+1.0f/(i*i);for (i=10000;i>=1;i--)s2=s2+1.0f/(i*i);d1=(float)(fabs(1.644834-s1));d2=(float)(fabs(1.644834-s2));printf("正向求和结果为%f\n 误差为%f\n\n",s1,d1);printf("反向求和结果为%f\n 误差为%f\n\n",s2,d2);if(d1<d2)printf("正向求和误差小于负向求和误差\n");else if(d1==d2)printf("正向求和误差等于负向求和误差\n"); elseprintf("正向求和误差大于负向求和误差\n");}实验结果：实验分析：第一次做数值实验，又一次使用C语言编程，没有了刚学习C语言的艰难，能够将实验步骤转换成流程图并编写出完整的实验代码，在经过多次调试、改正后得到正确的程序和结果。

这个实验较简单，计算误差时如果输入数据有误差，而在计算过程中舍入误差不增长，则称此算法是稳定的，否则称此算法是数值不稳定的，减少运算次数可以减小舍入误差。

数值分析上机实验报告

一、实验目的通过本次上机实验，掌握数值分析中常用的算法，如二分法、牛顿法、不动点迭代法、弦截法等，并能够运用这些算法解决实际问题。

同时，提高编程能力，加深对数值分析理论知识的理解。

二、实验环境1. 操作系统：Windows 102. 编程语言：MATLAB3. 实验工具：MATLAB数值分析工具箱三、实验内容1. 二分法求方程根二分法是一种常用的求方程根的方法，适用于连续函数。

其基本思想是：从区间[a, b]中选取中点c，判断f(c)的符号，若f(c)与f(a)同号，则新的区间为[a, c]，否则为[c, b]。

重复此过程，直至满足精度要求。

2. 牛顿法求方程根牛顿法是一种迭代法，适用于可导函数。

其基本思想是：利用函数在某点的导数值，求出函数在该点的切线方程，切线与x轴的交点即为方程的近似根。

3. 不动点迭代法求方程根不动点迭代法是一种迭代法，适用于具有不动点的函数。

其基本思想是：从初始值x0开始，不断迭代函数g(x)的值，直至满足精度要求。

4. 弦截法求方程根弦截法是一种线性近似方法，适用于可导函数。

其基本思想是：利用两点间的直线近似代替曲线，求出直线与x轴的交点作为方程的近似根。

四、实验步骤1. 二分法求方程根（1）编写二分法函数：function [root, error] = bisection(a, b, tol)（2）输入初始区间[a, b]和精度要求tol（3）调用函数计算根：[root, error] = bisection(a, b, tol)2. 牛顿法求方程根（1）编写牛顿法函数：function [root, error] = newton(f, df, x0, tol)（2）输入函数f、导数df、初始值x0和精度要求tol（3）调用函数计算根：[root, error] = newton(f, df, x0, tol)3. 不动点迭代法求方程根（1）编写不动点迭代法函数：function [root, error] = fixed_point(g, x0, tol)（2）输入函数g、初始值x0和精度要求tol（3）调用函数计算根：[root, error] = fixed_point(g, x0, tol)4. 弦截法求方程根（1）编写弦截法函数：function [root, error] = secant(f, x0, x1, tol)（2）输入函数f、初始值x0和x1，以及精度要求tol（3）调用函数计算根：[root, error] = secant(f, x0, x1, tol)五、实验结果与分析1. 二分法求方程根以方程f(x) = x^2 - 2 = 0为例，输入初始区间[a, b]为[1, 3]，精度要求tol 为1e-6。

同济大学数值分析matlab编程

同济⼤学数值分析matlab编程MATLAB 编程题库1.下⾯的数据表近似地满⾜函数21cx bax 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 ⼯作⽬录下已经有如下两个函数⽂件，写⼀个割线法程序，求出这两个函数精度为1010-的近似根，并写出调⽤⽅式：>> edit gexianfa.mfunction [x iter]=gexianfa(f,x0,x1,tol) iter=0;x=x1;while(abs(feval(f,x))>tol) iter=iter+1;x=x1-feval(f,x1).*(x1-x0)./(feval(f,x1)-feval(f,x0)); x0=x1;x1=x; end>> edit f.m function v=f(x) v=x.*log(x)-1;>> edit g.m function z=g(y) z=y.^5+y-1;>> [x1 iter1]=gexianfa('f',1,3,1e-10) x1 =1.7632 iter1 = 6>> [x2 iter2]=gexianfa('g',0,1,1e-10) x2 =0.7549 iter2 = 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;>> 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;>>for i=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; for i=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 ≤≤上求解初值问题。

同济大学数值分析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否则返回结束。