下载文档原格式
function x=tridiagsolver(a,b)[n,n]=size(a);for i=1:nif i==1l(i)=a(i,i);y(i)=b(i)/l(i);elsel(i)=a(i,i)-a(i,i-1)*u(i-1);y(i)=(b(i)-y(i-1)*a(i,i-1))/l(i);endif i<nu(i)=a(i,i+1)/l(i);endendx(n)=y(n);for j=n-1:-1:1x(j)=y(j)-u(j)*x(j+1);end拉格朗日function yh=lagrange(x,y,xh)n=length(x);m=length(xh);yh=zeros(1,m);c1=ones(n-1,1);c2=ones(1,m);for i=1:nxp=x([1:i-1 i+1:n]);yh=yh+y(i)*prod((c1*xh-xp'*c2)./(x(i)-xp'*c2));end线性x=[x1,x2] y=[y1.y2] xh=[xh]抛物线X=[x1,x2,x3] y=[y1，y2，y3] xh=[xh]牛顿差商（输入x,y为列向量）function [p,q]=chashang(x,y)n=length(x);p(:,1)=x;p(:,2)=y;for j=3:n+1p(1:n+2-j,j)=diff(p(1:n+3-j,j-1))./(x(j-1:n)-x(1:n+2-j)); endq=p(1,2:n+1)';三次样条x=[0 1 2 3];y=[0.2 0 0.5 2.0 1.5 -1];pp=csape(x,y,'complete')[breaks,coefs,npolys,ncoefs,dim]=unmkpp(pp)x=[0.24 0.65 0.95 1.24 1.73 2.01 2.23 2.52 2.77 2.99]';y=[0.23 -0.26 -1.1 -0.45 0.27 0.1 -0.29 0.24 0.56 1]';A=[log(x) cos(x) exp(x)];Z=A\y;a0=Z(1)a1=Z(2)a2=Z(3)x=[0 0.25 0.50 0.75 1.00];y=[1.00 1.284 1.6487 2.1170 2.7183];p=polyfit(x,y,2)a2=p(1)a1=p(2)a0=p(3)复合中点function I=fmid(fun,a,b,n)h=(b-a)/n;x=linspace(a+h/2,b-h/2,n);y=feval(fun,x);I=h*sum(y);复合梯形function I=ftrapz(fun,a,b,n)h=(b-a)/n;x=linspace(a,b,n+1);y=feval(fun,x);I=h*(0.5*y(1)+0.5*y(n+1)+sum(y(2:n)));复合辛普森function I=fsimpson(fun,a,b,n)h=(b-a)/n;x=linspace(a,b,2*n+1);y=feval(fun,x);I=h/6*(y(1)+y(2*n+1)+2*sum(y(3:2:2*n-1))+4*sum(y(2:2:2*n)));雅克比迭代function [x,iter]=jacobi(A,b,tol)D=diag(diag(A));L=D-tril(A);U=D-triu(A);x=zeros(size(b));for iter=1:500x=D\(b+U*x+L*x);error=norm(b-A*x)/norm(b);if (error<tol)break;endGS迭代function [x,iter]=GS(A,b,tol)D=diag(diag(A));L=D-tril(A);U=D-triu(A);x=zeros(size(b));for iter=1:500x=(D-L)\(b+U*x);error=norm(b-A*x)/norm(b);if (error<tol)break;endendSOR迭代function [x,iter]=SOR(A,b,omega,tol)D=diag(diag(A));L=D-tril(A);U=D-triu(A);x=zeros(size(b));for iter=1:500x=(D-omega*L)\(omega*b+(1-omega)*D*x+omega*U*x); error=norm(b-A*x)/norm(b);if (error<tol)break;endend二分法function v=f(x)v=x^3-x-1;function x=erfenfa(f,a,b,tol)if nargin<4tol=1e-5;endfa=feval(f,a);fb=feval(f,b);while abs(a-b)>tolx=(a+b)/2;fx=feval(f,x);if sign(fx)==sign(fa)a=x;fa=fx;elseif sign(fx)==sign(fb)b=x;fb=fx;endend>> x=erfenfa('f',1,2)牛顿法function [x,it]=newton(f,g,x0,tol)it=0;done=0;while ~donex=x0-feval(f,x0)/feval(g,x0);it=it+1;done=(norm(x-x0)<=tol);if ~donex0=x;endendfunction r=f(x)r=polyval([1,2,10,-20],x);function r=g(x)p=polyder([1,2,10,-20]);r=polyval(p,x);牛顿下山法乘幂法function [t,y]=chengmifa(a,xinit,ep) v0=xinit;[tv,ti]=max(abs(v0));lam0=v0(ti);u0=v0/lam0;flag=0;while ~flagv1=a*u0;[tv,ti]=max(abs(v1));lam1=v1(ti);u0=v1/lam1;err=abs(lam0-lam1);if err<epflag=1endlam0=lam1;endt=lam1;y=u0;反幂法function [t,y]=fanmifa(a,xinit,ep)v0=xinit;[tv,ti]=max(abs(v0));lam0=v0(ti);u0=v0/lam0;flag=0;while ~flagv1=inv(a)*u0;[tv,ti]=max(abs(v1));lam1=v1(ti);u0=v1/lam1;err=abs(lam0-lam1);if err<epflag=1;endlam0=lam1;endt=1/lam1;y=u0;改进欧拉法function [x,y]=odeIEuler(f,y0,a,b,n) h=(b-a)/n;x=a:h:b;y(1)=y0;for i=1:nyp=y(i)+h*feval(f,x(i),y(i));yc=y(i)+h*feval(f,x(i+1),yp);y(i+1)=1/2*(yp+yc);endfunction r=f(x,y)r=y-2*x/y;。
数值分析中求解非线性方程的MATLAB求解程序1. fzero函数：fzero函数是MATLAB中最常用的求解非线性方程的函数之一、它使用了割线法、二分法和反复均值法等多种迭代算法来求解方程。
使用fzero函数可以很方便地求解单变量非线性方程和非线性方程组。
例如，要求解方程f(x) = 0，可以使用以下语法：``````2. fsolve函数：fsolve函数是MATLAB中求解多维非线性方程组的函数。
它是基于牛顿法的迭代算法来求解方程组。
使用fsolve函数可以非常方便地求解非线性方程组。
例如，要求解方程组F(x) = 0，可以使用以下语法：``````3. root函数：root函数是MATLAB中求解非线性方程组的函数之一、它采用牛顿法或拟牛顿法来求解方程组。
使用root函数可以非常方便地求解非线性方程组。
例如，要求解方程组F(x) = 0，可以使用以下语法：``````4. vpasolve函数：vpasolve函数是MATLAB中求解符号方程的函数。
它使用符号计算的方法来求解方程，可以得到精确的解。
vpasolve函数可以求解多变量非线性方程组和含有符号参数的非线性方程。
例如，要求解方程组F(x) = 0，可以使用以下语法：```x = vpasolve(F(x) == 0, x)```vpasolve函数会返回方程组的一个精确解x。
5. fsolve和lsqnonlin结合：在MATLAB中，可以将求解非线性方程转化为求解最小二乘问题的形式。
可以使用fsolve函数或lsqnonlin函数来求解最小二乘问题。
例如，要求解方程f(x) = 0，可以将其转化为最小二乘问题g(x) = min，然后使用fsolve或lsqnonlin函数来求解。
具体使用方法可以参考MATLAB官方文档。
6. Newton-Raphson法手动实现：除了使用MATLAB中的函数来求解非线性方程，还可以手动实现Newton-Raphson法来求解。
数值分析课程介绍
课程代码（学校统一编制）
课程名称数值分析
英文名称Numerical Analysis
学分： 2 修读期：大三上学期
授课对象：理工科
课程主任：洪晓英、讲师、硕士
课程简介
《Numerical Analysis》是理工科专业基础选修课。
它主要介绍各种数值方法来解决形式比较复杂的各种数学问题。
通过本课程的学习，使学生了解和掌握这门课程所涉及的各种常用的数值计算公式、数值方法的构造原理及适用范围，并通过本课学到一些现代数学的概念,为今后用计算机去有效地解决实际的科研问题及进入现代数学打下基础。
主要包括：（1）引论
（2）线性方程组的求解
（3）插值法与最小二乘法
（4）数值积分与微分
（5）常微分方程的数值解法
（6）逐次逼近法
实践教学环节（如果有）
学习计算方法的过程中，进行重要的实验（上机）是必不可少的。
通过实验一方面加深对计算方法的理解，另一方面培养学生的解决实际问题的能力。
本课程有实验（上机）的教学安排,内容以教材附录中的上机实习参考题为主，共18学时。
要求学生熟悉至少一门数学软件平台（Mathematica/Matleb/Maple）和至少一种编程语言，能够编程实现几种重要的计算方法，至少做有求解足够规模问题的大作业4－5次。
课程考核
课外作业 10%，上机实验 20%；期末闭卷考试 70%。
指定教材
计算机数值方法，施吉林，高等教育出版社，2005年3月，第2版
参考书目
【1】计算方法，易大义，浙江大学出版社，2002年6月，第2版
【2】现代数值计算方法，肖筱南，北京大学出版社，2003年7月，第1版。
拉格朗日插值#include<iostream>#include<cmath>using namespace std;int lagrange(double *px,double *py,double *pcoeff,int iNum) { double *pk,*ptemp_coeff;int i,j,k,kk;pk=new double [iNum];ptemp_coeff=new double [iNum];for(i=0;i<iNum;i++) {pcoeff[i]=0.0;ptemp_coeff[i]=0.0;}for(i=0;i<iNum;i++) {pk[i]=1.0;for(j=0;j<iNum;j++) {if(i!=j) {pk[i]*=px[i]-px[j]; }}}for(i=0;i<iNum;i++) {ptemp_coeff[0]=1.0;for(j=0,k=0;j<iNum;j++) {if(j!=i) {for(k++,kk=k;kk>0;kk--) {ptemp_coeff[kk]-=ptemp_coeff[kk-1]*px[j];}}}for(j=0;j<iNum;j++)pcoeff[j]+=py[i]*ptemp_coeff[j]/pk[i];ptemp_coeff[j]=0.0;}}delete pk,ptemp_coeff;return 0;}double lagcaculate(double *pcoeff,double dxx,int iNum) { double dsum;int i;dsum=0.0;for(i=0;i<iNum;i++) {dsum=dsum*dxx+pcoeff[i];}return dsum;}int main() {int iNum,i,iisrandom;double *px,*py,*pcoeff;cout<<"本程序求解拉格朗日插值求解函数近似值" <<endl<<"请输入已知点个数:"<<endl;cin>>iNum;cout<<"是否进行自动测试（测试函数为正弦函数）,是输入1，否输入0：";cin>>iisrandom;px=new double [iNum];py=new double [iNum];pcoeff=new double [iNum];if(iisrandom==0) {cout<<"请输入iNum个已知点的横坐标:"<<endl;for(i=0;i<iNum;i++) {cin>>px[i];}cout<<"请输入iNum个已知点的纵坐标:"<<endl;for(i=0;i<iNum;i++)cin>>py[i];}else {for(i=0;i<iNum;i++) {px[i]=i;py[i]=sin(i);}}lagrange(px,py,pcoeff, iNum);int testnum;cout<<"输入测试值数量";cin>>testnum;double *pxvalue;pxvalue=new double [testnum];double *pyvalue;pyvalue=new double [testnum];cout<<"输入"<<testnum<<"个需要测试的值"<<endl;for(i=0;i<testnum;i++) {cin>>pxvalue[i];}for(i=0;i<testnum;i++) {pyvalue[i]=lagcaculate(pcoeff,pxvalue[i],iNum);py[i]=sin(pxvalue[i]);}cout<<"拉格朗日公式求得的值"<<" "<<"标准函数值为"<<endl;for(i=0;i<testnum;i++) {cout<<pyvalue[i]<<" "<<py[i]<<endl; }delete px,py,pcoeff,pxvalue,pyvalue;return 0;}牛顿插值#include <iostream.h>#include <math.h>void main() {int n,i,j;double A[50][50];double x[50],y[50];double K=1,X=0,N=0,P;cout<<"请输入所求均差阶数:";cin>>n;for(i=0;i<=n;i++) {cout<<"请输入x"<<i<<"=";cin>>x[i];cout<<"请输入y"<<i<<"=";cin>>y[i];A[i][0]=x[i]; A[i][1]=y[i];}for(j=2;j<=n+1;j++) {for(i=1;i<=n;i++) {A[i][j]=(A[i][j-1]-A[i-1][j-1])/(A[i][0]-A[i-j+1][0]);}}for(i=0;i<=n;i++) {cout<<"输出第"<<i<<"阶均差为："<<A[i][i+1]<<endl;}cout<<"请所要代入计算的x的值：X=";cin>>X;for(i=0;i<n;i++) {K*=X-x[i];N+=A[i+1][i+2]*K;P=A[0][1]+N;}c o u t<<"所要求的函数值为：y="<<P<<endl;}复化梯形公式double f(double x){if(x==0) return 1;else return (sin(x)/x);}double FuhuaTixing(int n,double a,double b){double h = (b-a)/n;double x = a;double s = 0;for(int k=0; k< n-1; k++){x += h;s += f(x);}double T = (f(a)+s*2+f(b))*h/2;return T;}int main(){char ans='n';do{cout<<"请输入积分区间(a,b)："<<endl;double a;double b;cin>>a>>b;cout<<"请输入等分份数n： "<<endl;int n; cin>>n;cout<<"由复化梯形公式球的结果："<<FuhuaTixing(n,a,b)<<endl; cout<<"是否要继续？(y/n)";cin>>ans;}while(ans == 'y');return 0;}改进欧拉#include"stdio.h"#include"iostream"using namespace std;double x,x1,y,y1;int main(){double h;int n;cout<<endl<<"input x0=";cin>>x0;cout<<endl<<"input y0=";cin>>y0;cout<<endl<<"intput h=";cin>>h;cout<<endl<<"intput n=";cin>>n;mend_euler euler(h,n);getch();return 1;}class mend_euler{public:mend_euler(double h,int n);double f(double x,double y);private:double h;int n;};mend_euler::mend_euler(double a,int b){int i=1; h=a; n=b;while(i<=n) {x1=x0+h;y1=y0+h/2*(f(x0,y0)+f(x1,y0+h*f(x0,y0)));cout<<endl;cout<<"x1="<<x1<<" y1="<<y1<<" y="<<x0/(1+x0*x0)<<" e="<<y1-x0/(1+x0*x0) <<endl;i++;x0=x1;y0=y1;}}double mend_euler::f(double x,double y){ return 4*x*y/(1+x*x)+1; }四阶龙格库塔#include <iostream.h> #include <math.h>double f1(double t, double x, double y, double z) { return -x + 6 * y + 2 * z; }double f2(double t, double x, double y, double z) { return -y + 3 * z + 2 * sin(t); }double f3(double t, double x, double y, double z) { return -z + pow(t,2) * exp(-t) + cos(t); } double RungeKutta(double x0, double y0, double z0, double a, double b, int n, double x[], double y[], double z[], double (*f1)(double,double,double,double), double (*f2)(double,double,double,do uble), double (*f3)(double,double,double,double)){double k1, k2, k3, k4double m1, m2, m3, m4double n1, n2, n3, n4double h; double t;int count; count =0;h = (b - a) / n;x[0] = x0;y[0] = y0;z[0] = z0;for(int i = 1;i<=n;i++) {t = a + (i - 1) * h;k1 = f1(t, x[i - 1], y[i - 1], z[i - 1]);m1 = f2(t, x[i - 1], y[i - 1], z[i - 1]);n1 = f3(t, x[i - 1], y[i - 1], z[i - 1]);k2 = f1(t+h/2.0, x[i - 1] + k1 * h / 2.0, y[i - 1] + m1 * h / 2.0, z[i - 1] + n1 * h / 2.0);m2 = f2(t+h/2.0, x[i - 1] + k1 * h / 2.0, y[i - 1] + m1 * h / 2.0, z[i - 1] + n1 * h / 2.0);n2 = f3(t+h/2.0, x[i - 1] + k1 * h / 2.0, y[i - 1] + m1 * h / 2.0, z[i - 1] + n1 * h / 2.0); k3 = f1(t+h/2.0, x[i - 1] + k2 * h / 2.0, y[i - 1] + m2 * h / 2.0, z[i - 1] + n2 * h / 2.0);m3 = f2(t+h/2.0, x[i - 1] + k2 * h / 2.0, y[i - 1] + m2 * h / 2.0, z[i - 1] + n2 * h / 2.0);n3 = f3(t+h/2.0, x[i - 1] + k2 * h / 2.0, y[i - 1] + m2 * h / 2.0, z[i - 1] + n2 * h / 2.0);k4 = f1(t+h/2.0, x[i - 1] + k3 * h / 2.0, y[i - 1] + m3 * h / 2.0, z[i - 1] + n3 * h / 2.0);m4 = f2(t+h/2.0, x[i - 1] + k3 * h / 2.0, y[i - 1] + m3 * h / 2.0, z[i - 1] + n3 * h / 2.0);n4 = f3(t+h/2.0, x[i - 1] + k3 * h / 2.0, y[i - 1] + m3 * h / 2.0, z[i - 1] + n3 * h / 2.0);x[i] = x[i - 1] + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;y[i] = y[i - 1] + h * (m1 + 2 * m2 + 2 * m3 + m4) / 6.0;z[i] = z[i - 1] + h * (n1 + 2 * n2 + 2 * n3 + n4) / 6.0; ++count;cout<<"t= "<<t+0.05<<" "<<"x= "<<x[i]<<" "<<"y= "<<y[i]<<" "<<"z= "<<z[i]<<endl;}cout<<"共计算的次数为："<<count<<endl; return 0;}void main() { double x[110]; double y[110];double z[110];RungeKutta(1.0,-1.0,0.0,0.0,5.0,100,x,y,z,f1,f2,f3); }弦截法#include<iostream>#include<math.h>using namespace std;class imaroot{public:virtual float f(float x)=0;virtual float x(float x1,float x2)=0;};class root:public imaroot{public:float f(float x){float y;y=((x-5.0)*x+16.0)*x-80.0;return y;}float x(float x1,float x2){float x;x=(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1));return x;}};int main(){root a;float x,c,d;cout<<"输入任意2点横坐标(x1<x2)"<<endl;cin>>c>>d;if(a.f(c)*a.f(d)>=0){do{cout<<"f(x1)与f(x2)同号重新输入"<<endl;cin>>c>>d;}while(a.f(c)*a.f(d)>=0);}if(a.f(c)*a.f(d)<0)x=a.x(c,d);if(a.f(x)*a.f(c)>0){do{x=a.x(x,d);}while(abs(a.f(x))>=0.0001);}else{do{x=a.x(c,x);}while(abs(a.f(x))>=0.0001);}cout<<"方程根是:"<<x<<endl;return 0;}用C语言实现高斯-赛德尔迭代方法，源程序代码：#include "stdlib.h"#include "stdio.h"#include "conio.h"#include "string.h"#include "math.h"#define N 100float Table(int n,float a[N][N],float b[N]){ int i,j; float c[N][N];printf("Please input the matrix A by row!\n"); label:for(i=0;i<n;i++) { printf("Row %d:",i+1); for(j=0;j<n;j++)scanf("%f",&a[i][j]); }printf("Please input the vector b:");for(i=0;i<n;i++)scanf("%f",&b[i]);for(i=0;i<n;i++)for(j=0;j<n;j++) {if(i==j) {c[i][j]=0;continue; }c[i][j]=-a[i][j]/a[i][i]; } printf("\nThe matrix A and vector b:\n"); for(i=0;i<n;i++) { for(j=0;j<n;j++)printf("%10.5f",a[i][j]);printf("%10.5f",b[i]);printf("\n");}printf("\nThe Gauss-Seidel iterative scheme:\n");for(i=0;i<n;i++){for(j=0;j<n;j++)printf("%10.5f",c[i][j]);printf("%10.5f",b[i]/a[i][i]);printf("\n"); } }float init_vec(int n,float x[N]){I nt i;printf("Please input the initial iteration vector x:");for(i=0;i<n;i++)scanf("%f",&x[i]);printf("\nThe initial iteration vector x:\n"); for(i=0;i<n;i++) printf("%10.5f",x[i]);printf("\n");}float gs(int n,float a[N][N],float b[N],float x[N]) {int i,j,k;float tmp1,tmp2,x2[N];for(k=0;k<10001;k++){for(i=0;i<n;i++)x2[i]=x[i];for(i=0;i<n;i++){tmp1=0.0;tmp2=0.0;for(j=0;j<i;j++)tmp1+=a[i][j]*x[j];for(j=i+1;j<n;j++)tmp2+=a[i][j]*x2[j];x[i]=(b[i]-tmp1-tmp2)/a[i][i];}for(i=0,j=0;i<n;i++)if(fabs(x2[i]-x[i])<0.00001) j++;if(j==n){printf("\nThis Gauss-Seidel iterative scheme is convergent!");printf("\nNumber of iterations: %d",k+1);break; }if(k==499){printf("\nThis Jacobi iterative scheme may be not convergent!"); break; }}printf("\nThe results:\n");for(i=0;i<n;i++)printf("%12.7f",x[i]);}int main() {int n;float x[N],a[N][N],b[N];printf("Input n:");scanf("%d",&n);Table(n,a,b);init_vec(n,x);gs(n,a,b,x);getch();}高斯消去#include "Stdio.h"#include "Conio.h"/*L是矩阵的行减1,从程序上看是最外层循环的次数N 对应矩阵的行数，M对应矩阵的列数可以通过改变L、N、M来控制矩的阶数*/#define L 3#define N 4#define M 5void gauss(double a[N][M],double x[N]){int i,j,l,n,m,k=0;double temp[N];/*第一个do-while是将增广矩阵消成上三角形式*/do{n=0;for(l=k;l<L;l++)temp[n++]=a[l+1][k]/a[k][k];for(m=0,i=k;i<N;i++,m++)for(j=k;j<M;j++)a[i+1][j]-=temp[m]*a[k][j];k++;}while(k<N) ;/*第二个do-while是将矩阵消成对角形式，并且重新给k赋值*/k=L-1;do{n=0;for(l=k;l>=0;l--)temp[n++]=a[k-l][k+1]/a[k+1][k+1];for(m=0,i=k;i>=0;i--,m++)for(j=k;j<M;j++)a[k-i][j]-=temp[m]*a[k+1][j];k--;}while(k>=0) ;/*下一个for是解方程组*/for(i=0;i<N;i++)x[i]=a[i][N]/a[i][i];}void menu(){printf("\n _ _ _ _ _\n");printf(" 1.operation\n");printf(" 2.exit");printf("\n _ _ _ _ _\n");}main(){int i,j,choose;double a[N][M]={0},answer[N];clrscr();while(1){leep:menu();scanf("%d",&choose);switch(choose){case 1:printf("!!The size of Maxrix is %d * %d,each line enter %d element:\n ",N,M,M); for(i=0;i<N;i++){printf("Enter the Matrix's %d line:\n",i);for(j=0;j<N+1;j++)scanf("%lf",&a[i][j]);}printf("\nthe corss matrix is:\n_ _ _ _ _\n");gauss(a,answer);for(i=0;i<N;i++){for(j=0;j<M;j++)printf("%-2lf ",a[i][j]);putchar('\n');}printf("_ _ _ _ _\nthe solve is:\n");for(i=0;i<N;i++)printf("x%d=%lf\n",i+1,answer[i]); break;case 2:exit(0);break;default:printf("input error:\n");goto leep;}}getch();}。
数值分析⽅法计算定积分⽤C语⾔实现⼏种常⽤的数值分析⽅法计算定积分，代码如下：1 #include <stdio.h>2 #include <string.h>3 #include <stdlib.h>4 #include <math.h>56#define EPS 1.0E-14 //计算精度7#define DIV 1000 //分割区间，值越⼤，精确值越⾼8#define ITERATE 20 //⼆分区间迭代次数，区间分割为2^n，初始化应该⼩⼀点，否则会溢出910#define RECTANGLE 0 //矩形近似11#define TRAPEZOID 1 //梯形近似12#define TRAPEZOID_FORMULA 2 //递推梯形公式13#define SIMPSON_FORMULA 3 //⾟普森公式14#define BOOL_FORMULA 4 //布尔公式1516double function1(double);17double function2(double);18double function3(double);19void Integral(int, double f(double), double, double); //矩形, 梯形逼近求定积分公式20void Trapezoid_Formula(double f(double), double a, double b); //递推梯形公式21void Simpson_Formula(double f(double), double a, double b); //⾟普森公式22void Bool_Formula(double f(double), double a, double b); //布尔公式23void Formula(int, double f(double), double, double);2425double function1(double x)26 {27double y;28 y = cos(x);29return y;30 }3132double function2(double x)33 {34double y;35 y=1/x;36// y = 2+sin(2 * sqrt(x));37return y;38 }3940double function3(double x)41 {42double y;43 y = exp(x);44return y;45 }4647int main()48 {49 printf("y=cos(x) , x=[0, 1]\n");50 printf("Rectangle : "); Integral(RECTANGLE, function1, 0, 1);51 printf("Trapezoid : "); Integral(TRAPEZOID, function1, 0, 1);52 Formula(TRAPEZOID_FORMULA, function1, 0, 1);53 Formula(SIMPSON_FORMULA, function1, 0, 1);54 Formula(BOOL_FORMULA, function1, 0, 1);55 printf("\n");5657 printf("y=1/x , x=[1, 5]\n");58 printf("Rectangle : "); Integral(RECTANGLE, function2, 1, 5);59 printf("Trapezoid : "); Integral(TRAPEZOID, function2, 1, 5);60 Formula(TRAPEZOID_FORMULA, function2, 1, 5);61 Formula(SIMPSON_FORMULA, function2, 1, 5);62 Formula(BOOL_FORMULA, function2, 1, 5);63 printf("\n");6465 printf("y=exp(x) , x=[0, 1]\n");66 printf("Rectangle : "); Integral(RECTANGLE, function3, 0, 1);67 printf("Trapezoid : "); Integral(TRAPEZOID, function3, 0, 1);68 Formula(TRAPEZOID_FORMULA, function3, 0, 1);69 Formula(SIMPSON_FORMULA, function3, 0, 1);70 Formula(BOOL_FORMULA, function3, 0, 1);7172return0;73 }74//积分：分割-近似-求和-取极限75//利⽤黎曼和求解定积分76void Integral(int method, double f(double),double a,double b) //f(double)，f(x), double a，float b，区间两点77 {79double x;80double sum = 0; //矩形总⾯积81double h; //矩形宽度82double t = 0; //单个矩形⾯积8384 h = (b-a) / DIV;8586for(i=1; i <= DIV; i++)87 {88 x = a + h * i;89switch(method)90 {91case RECTANGLE :92 t = f(x) * h;93break;94case TRAPEZOID :95 t = (f(x) + f(x - h)) * h / 2;96break;97default:98break;99 }100 sum = sum + t; //各个矩形⾯积相加101 }102103 printf("the value of function is %.10f\n", sum);104 }105106//递推梯形公式107void Trapezoid_Formula(double f(double), double a, double b)//递推梯形公式108 {109 unsigned int i, j, M, J=1;110double h;111double F_2k_1 = 0;112double T[32];113double res = 0;114115 T[0] = (f(a) + f(b)) * (b-a)/2;116for(i=0; i<ITERATE; i++)117 {118 F_2k_1 = 0;119 J = 1;120121for(j=0; j<=i; j++) //区间划分122 J <<= 1; //2^J123 h = (b - a) /J; //步长124//M = J/2; //2M表⽰将积分区域划分的块数125for(j=1; j <= J; j += 2) //126 {127 F_2k_1 += f(a + j*h); //f(2k-1)累加和128 }129 T[i+1] = T[i]/2 + h * F_2k_1; //递推公式130 res = T[i+1];131132if(fabs(T[i+1] - T[i]) < EPS)133if(i < 3) continue;134else break;135 }136137 printf("Trapezoid_Formula : the value of function is %.10f\n", res);138//return res;139 }140//⾟普森公式141void Simpson_Formula(double f(double), double a, double b) //⾟普森公式142 {143 unsigned int i, j, M, J=1;144double h;145double F_2k_1 = 0;146double T[32], S[32];147double res_T = 0, res_S = 0;148149 T[0] = (f(a) + f(b)) * (b-a)/2;150for(i=0; i<ITERATE; i++)151 {152 F_2k_1 = 0;153 J = 1;154155for(j=0; j<=i; j++) //区间划分156 J <<= 1; //2^J157 h = (b - a) /J; //步长158//M = J/2; //2M表⽰将积分区域划分的块数159for(j=1; j <= J; j += 2) //160 {161 F_2k_1 += f(a + j*h); //f(2k-1)累加和163 T[i+1] = T[i]/2 + h * F_2k_1; //递推梯形公式164 res_T = T[i+1];165166if(fabs(T[i+1] - T[i]) < EPS)167if(i < 3) continue;168else break;169 }170171for(i=1; i<ITERATE; i++)172 {173 S[i] = (4 * T[i] - T[i-1]) / 3; //递推⾟普森公式174 res_S = S[i];175if(fabs(S[i] - S[i-1]) < EPS)176if(i < 3) continue;177else break;178 }179180 printf("Simpson_Formula : the value of function is %.10f\n", res_S);181//return res_S;182 }183184//布尔公式185void Bool_Formula(double f(double), double a, double b) //布尔公式186 {187 unsigned int i, j, M, J=1;188double h;189double F_2k_1 = 0;190double T[32], S[32], B[32];191double res_T = 0, res_S = 0, res_B;192193 T[0] = (f(a) + f(b)) * (b-a)/2;194for(i=0; i<ITERATE; i++)195 {196 F_2k_1 = 0;197 J = 1;198199for(j=0; j<=i; j++) //区间划分200 J <<= 1; //2^J201 h = (b - a) /J; //步长202//M = J/2; //2M表⽰将积分区域划分的块数203for(j=1; j <= J; j += 2) //204 {205 F_2k_1 += f(a + j*h); //f(2k-1)累加和206 }207 T[i+1] = T[i]/2 + h * F_2k_1; //递推公式208 res_T = T[i+1];209210if(fabs(T[i+1] - T[i]) < EPS)211if(i < 3) continue;212else break;213 }214215for(i=1; i<ITERATE; i++)216 {217 S[i] = (4 * T[i] - T[i-1]) / 3; //递推⾟普森公式218 res_S = S[i];219if(fabs(S[i] - S[i-1]) < EPS)220if(i < 3) continue;221else break;222 }223224for(i=1; i<ITERATE; i++)225 {226 B[i] = (16 * S[i] - S[i-1]) / 15; //递推布尔公式227 res_B = B[i];228if(fabs(B[i] - B[i-1]) < EPS)229if(i < 3) continue;230else break;231 }232233 printf("Bool_Formula : the value of function is %.10f\n", res_B);234//return res_B;235 }236237//采⽤⼆分区间的⽅法迭代，实际运⾏速度很慢238void Formula(int method, double f(double), double a, double b)//递推梯形公式239 {240 unsigned int i, j, M, J=1;241double h;242double F_2k_1 = 0;243double T[32], S[32], B[32];244double res_B = 0, res_S = 0, res_T = 0;247for(i=0; i<ITERATE; i++)248 {249 F_2k_1 = 0;250 J = 1;251252for(j=0; j<=i; j++) //区间划分253 J <<= 1; //2^J254 h = (b - a) /J; //步长255//M = J/2; //2M表⽰将积分区域划分的块数256for(j=1; j <= J; j += 2) //257 {258 F_2k_1 += f(a + j*h); //f(2k-1)累加和259 }260 T[i+1] = T[i]/2 + h * F_2k_1; //递推公式261 res_T = T[i+1];262263if(fabs(T[i+1] - T[i]) < EPS)264if(i < 3) continue;265else break;266 }267268switch(method)269 {270default :271case TRAPEZOID_FORMULA :272 printf("Trapezoid_Formula : the value of function is %.10f\n", res_T); 273//return res_T;274break;275case SIMPSON_FORMULA :276for(i=1; i<ITERATE; i++)277 {278 S[i] = (4 * T[i] - T[i-1]) / 3;279 res_S = S[i];280if(fabs(S[i] - S[i-1]) < EPS)281if(i < 3) continue;282else break;283 }284 printf("Simpson_Formula : the value of function is %.10f\n", res_S); 285//return res_S;286break;287case BOOL_FORMULA :288for(i=1; i<ITERATE; i++)289 {290 S[i] = (4 * T[i] - T[i-1]) / 3;291 res_S = S[i];292if(fabs(S[i] - S[i-1]) < EPS)293if(i < 3) continue;294else break;295 }296for(i=1; i<ITERATE; i++)297 {298 B[i] = (16 * S[i] - S[i-1]) / 15;299 res_B = B[i];300if(fabs(B[i] - B[i-1]) < EPS)301if(i < 3) continue;302else break;303 }304305 printf("Bool_Formula : the value of function is %.10f\n", res_B); 306//return res_B;307break;308 }309310 }测试结果：。
课程设计任务书实验方法与理论方法是推动科学技术发展的两大基本方法，但有局限性。
许多研究对象，由于空间或时间的限制，既不可能用理论精确描述，也不能用实验手段实现。
数值模拟或称为科学计算突破了实验和理论科学的局限，在科技发展中起到越来越重要的作用。
可以认为，科学计算已于实验、理论一起成为科学方法上不可或缺的三个主要手段。
计算数学的研究是科学计算的主要组成部分，而数值分析则是计算数学的核心。
数值计算是研究使用计算机来解决各种数学问题的近似计算方法与理论，其任务是提供在计算机上可解的、理论可靠的、计算复杂性低的各种常用算法。
数值分析的主要内容：1）、数值代数：求解线性和非线性方程组的解，分直接方法和间接方法两大类；2）、插值、曲线拟合和数值逼近；3）、数值微分和数值积分；4）、常微分和偏微分方程数值解法。
本文主要通过Matlab软件，对数值分析中的一些问题进行求解，如列主元Gauss消去法，Lagrange插值多项式，复化Simpson公式，Runge-Kutta方法以及数值分析在实际问题中的应用，并在求解的过程中更加熟识这门课程的主要内容，以及加强对课程知识的掌握。
在学习与设计计算方法时，从数学理论角度，学会分析方法的误差、收敛性和稳定性，保证计算方法的准确性；从实际应用的角度出发，掌握计算方法的结构与流程，能够把计算方法转换为可在计算机上直接处理的程序，保证算法的可用性。
关键词：列主元Gauss消去法；Lagrange插值；复化Simpson公式；Runge-Kutta实验一列主元Gauss消去法 (1)1.1 实验目的 (1)1.2 基本原理 (1)1.3 实验内容 (2)1.4 实验结论 (3)实验二拉格朗日插值多项式 (4)2.1 实验目的 (4)2.2 基本原理 (4)2.3 实验内容 (4)2.4 实验结论 (9)实验三复化Simpson求积公式 (10)3.1 实验目的 (10)3.2 基本原理 (10)3.3 实验内容 (10)3.4 实验结论 (12)实验四龙格-库塔(Runge-Kutta)方法 (13)4.1 实验目的 (13)4.2 基本原理 (13)4.3 实验内容 (14)4.4 实验结论 (15)实验五数值方法实际应用 (16)5.1 实验目的 (16)5.2 基本原理 (16)5.3 实验内容 (16)5.4 实验结论 (22)参考文献 (23)实验一 列主元Gauss 消去法1.1 实验目的1) 理解列主元消去法的原理；2) 熟悉列主元消去法的计算步骤，能用代码编写； 3) 解决实际问题。
注：凡“in.open("C:\\Users\\Administrator\\Desktop\\data-拉格朗日插值.txt")”等里面的路径自己合理修改。
1-1.拉格朗日插值（用一系列已知点求近似函数，并估计未知点的函数值）#include<fstream.h>#include<stdlib.h>int n; //插值节点个数double *X,*Y;//插值节点数组double x;//估计点void read(){double temp;ifstream in;in.open("C:\\Users\\Administrator\\Desktop\\data-拉格朗日插值.txt");if(!in){cout<<"can not open the file !"<<endl;exit(0);}in>>n;X=new double[n];Y=new double[n];for(int i=0;i<n;i++){in>>temp;X[i]=temp;cout<<temp<<" ";in>>temp;Y[i]=temp;cout<<temp<<"\n";}in>>temp;x=temp;cout<<"x="<<temp<<endl;in.close();}void main(){double t,y=0;read();for(int k=0;k<n;k++){t=1;for(int i=0;i<n;i++)if(k!=i)t=(x-X[i])/(X[k]-X[i])*t;y=y+t*Y[k];}cout<<"y="<<y<<endl;}1-2.埃特金逐步插值法（同上）#include<fstream.h>#include<stdlib.h>double *X,*Y; //插值节点（x，y）double x; //带估点int n; //N个插值点，插n-1次。
数 学 实 验 报 告P318习题1、给定初值问题（1）21',12(1)1y y x x xy ⎧=-≤≤⎪⎨⎪=⎩ 编制函数文件： Fun.mFunction f=fun(x,y)f=1/x.^2-y/x ------- % y’的表达式 在MA TLAB 窗口输入命令：[x,y]=ode45(‘fun’,[1,2],1) ---------- % ode45表示采用四阶、五阶Runge-Kutta 方法求解常微分方程，fun 代表导函数，[1,2]表示x 的区间范围 ，y 表示函数的初值，y(1)=1. 图形见下图1 >> x=--------% x 表示自变量Columns 1 through 181.0000 1.0250 1.0500 1.0750 1.1000 1.1250 1.1500 1.17501.2000 1.2250 1.2500 1.2750 1.3000 1.3250 1.3500 1.3750 1.4000 1.4250 Columns 19 through 361.4500 1.4750 1.5000 1.5250 1.5500 1.5750 1.6000 1.62501.6500 1.6750 1.7000 1.7250 1.7500 1.7750 1.8000 1.8250 1.8500 1.8750 Columns 37 through 411.9000 1.9250 1.9500 1.97502.0000 >> y = ---------% y 表示 Columns 1 through 181.0000 0.9997 0.9988 0.9975 0.9957 0.9936 0.9911 0.98830.9853 0.9820 0.9785 0.9749 0.9710 0.9671 0.9630 0.9589 0.9546 0.9503 Columns 19 through 360.9459 0.9415 0.9370 0.9325 0.9279 0.9233 0.9188 0.91420.9096 0.9050 0.9004 0.8958 0.8912 0.8866 0.8821 0.8776 0.8731 0.8686Columns 37 through 410.8641 0.8597 0.8553 0.8509 0.84661 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.920.840.860.880.90.920.940.960.9811.02图1 y-x 的图形（2）2'50502,011(0)3y y x x x y ⎧=-++≤≤⎪⎨=⎪⎩编制函数文件 Fun.mFunction f=fun(x,y) 在MA TLAB 输入命令[x,y]=ode45(‘fun’,[0,1],1/3) 图形见下图2 >> x'=Columns 1 through 180 0.0010 0.0020 0.0030 0.0040 0.0081 0.0122 0.01630.0204 0.0244 0.0283 0.0322 0.0361 0.0401 0.0440 0.0479 0.0518 0.0558Columns 19 through 360.0597 0.0637 0.0676 0.0716 0.0756 0.0796 0.08360.0877 0.0919 0.0960 0.1002 0.1047 0.1093 0.11380.1183 0.1237 0.1291 0.1345Columns 37 through 540.1399 0.1462 0.1525 0.1589 0.1652 0.1725 0.17970.1870 0.1943 0.2023 0.2103 0.2183 0.2262 0.23480.2433 0.2518 0.2604 0.2693Columns 55 through 720.2783 0.2873 0.2963 0.3056 0.3150 0.3244 0.33380.3435 0.3533 0.3631 0.3729 0.3830 0.3931 0.40330.4134 0.4239 0.4344 0.4449Columns 73 through 900.4553 0.4662 0.4770 0.4878 0.4986 0.5098 0.52090.5321 0.5432 0.5546 0.5661 0.5775 0.5890 0.60070.6124 0.6242 0.6359 0.6479Columns 91 through 1080.6599 0.6719 0.6839 0.6962 0.7085 0.7207 0.73300.7455 0.7580 0.7706 0.7831 0.7958 0.8086 0.82130.8341 0.8470 0.8600 0.8730Columns 109 through 1210.8860 0.8991 0.9123 0.9255 0.9387 0.9521 0.96540.9788 0.9922 0.9941 0.9961 0.9980 1.0000>> y =Columns 1 through 180.3333 0.3170 0.3015 0.2867 0.2727 0.2220 0.18090.1475 0.1204 0.0992 0.0818 0.0677 0.0561 0.04660.0389 0.0327 0.0277 0.0236Columns 19 through 360.0204 0.0179 0.0159 0.0144 0.0133 0.0126 0.01210.0118 0.0118 0.0120 0.0123 0.0127 0.0134 0.01410.0149 0.0160 0.0172 0.0185Columns 37 through 540.0199 0.0216 0.0234 0.0254 0.0274 0.0298 0.03230.0350 0.0378 0.0409 0.0442 0.0476 0.0512 0.05510.0592 0.0634 0.0678 0.0725Columns 55 through 720.0774 0.0825 0.0878 0.0934 0.0992 0.1052 0.11140.1180 0.1248 0.1318 0.1391 0.1467 0.1545 0.1626 0.1709 0.1797 0.1886 0.1979Columns 73 through 900.2074 0.2173 0.2275 0.2380 0.2487 0.2598 0.27130.2831 0.2951 0.3076 0.3204 0.3336 0.3470 0.3608 0.3750 0.3896 0.4045 0.4197Columns 91 through 1080.4354 0.4515 0.4679 0.4846 0.5018 0.5195 0.53750.5557 0.5745 0.5938 0.6134 0.6333 0.6536 0.6746 0.6959 0.7174 0.7394 0.7621Columns 109 through 1210.7851 0.8083 0.8321 0.8566 0.8814 0.9063 0.93180.9581 0.9847 0.9886 0.9924 0.9963 1.000200.10.20.30.40.50.60.70.80.910.20.40.60.811.21.4图2 y-x 的图形收获：通过上课所学和以上练习，了解以及掌握了表示用二阶、三阶和用四阶、五阶的Runge-Kutta 法解常微分方程的函数ode()、ode45（）的用法，可以使用MATLAB 解答简单的微分方程问题。
一、目录一、目录 (1)二、设计内容 (2)1.Gauss列主元素消去法解线性方程组的算法设计 (2)1.1 模块设计 (2)1.2主要程序代码 (2)1.3 运行结果 (5)2.Romberg算法 (6)2.1模块设计 (6)2.2 主要程序代码 (6)2.3 运行结果 (7)三、小结 (8)四、参考文献 (8)二、设计内容1.Gauss列主元素消去法解线性方程组的算法设计1.1 模块设计求任意阶（n<10）线性方程组的能力，主要包括一下模块：（1）在源程序中定义求解的方程阶数；（2）输入方程的系数矩阵和常数矩阵；（3）获取增广矩阵（4）Gauss列主元素消去：找主元、消元（分解）、回代、输出结果。
1.2主要程序代码#include <iostream.h>#include <iomanip.h>#include <math.h>#define n 4class Guss{private:float a[n][n];float b[n];float Array[n][n+1];public:Guss() //构造矩阵{cout<<"请输入"<<n<<"*"<<n<<"维系数矩阵:"<<endl;for(int i = 0;i<n;i++) // 输入系数矩阵{for(int j = 0;j<n;j++)cin>>a[i][j];}cout<<"请输入"<<1<<"*"<<n<<"维常数矩阵："<<endl;for(i = 0;i<n;i++) // 输入常数矩阵{cin>>b[i];}/* for(int i1 = 0;i1<n;i1++){for(int j1 =0;j1<n+1;j1++){cout<<a[i1][j1]<<setw(15);}cout<<b[i1];cout<<endl;}*/}void zgjz() //获得增广矩阵{// float c[n][n+1]; //定义增广矩阵// float *p;// p = &a[0][0];for(int i = 0;i<n;i++) //有系数矩阵和常数矩阵获得增广矩阵205 {for(int j = 0;j<n+1;j++){if(j == n)Array[i][j] = b[i];elseArray[i][j] = a[i][j];}}// return &c[0][0];// return *p; //返回增广矩阵cout<<"增广矩阵为："<<endl;for(int i1 = 0;i1<n;i1++){for(int j1 =0;j1<n+1;j1++){cout<<Array[i1][j1]<<setw(15);}cout<<endl;}}void Gjfc() //Guss列主元素消元法{float max,u,z; //定义三个float变量，max标记最大值，u用来交换，z用来计算最后的结果int flag; //flag用来标记下标for(int i = 0;i<n-1;i++){max = Array[i][i];flag = i;for(int k = i+1;k<n;k++){if(fabs(max) < fabs(Array[k][i])) //找主元{max = Array[k][i];flag = k;}}if(flag != i){for(int j = i;j<n+1;j++){u = Array[i][j];Array[i][j] = Array[flag][j];Array[flag][j] = u;}}for(int x = i+1;x<n;x++) //消元{u = Array[x][i]/max;for(int y = i;y<n+1;y++){Array[x][y] = Array[x][y]-u*Array[i][y];}}}for(int x=n-1;x>=0;x--) //回代{z=0;for(int y=x+1;y<n;y++)z=z+Array[x][y]*Array[y][n];Array[x][n]=(Array[x][n]-z)/(Array[x][x]);//计算结果}cout<<"Guss消元所得方程组的解为："<<endl;for(i = 0;i<n;i++) //输出结果cout<<Array[i][n]<<setw(15);cout<<endl;}};//类的结束int main(){Guss G;// float Array[n][n+1];G.zgjz();G.Gjfc();return 0;}1.3 运行结果2.Romberg算法2.1模块设计利用递推的复合梯形公式，并结合外推加速公式的Romberg算法。
Matlab 数值分析求高斯勒让德积分function [ result ] = gslrdjf(f, a,b,n ,GaussP,GaussA)%-----------------------------------------------------------------------------------%高斯勒让德积分(gslrdjf.m)%输入被积分函数式f，积分区域a到b，高斯点GaussP，高斯系数GaussA%直接输出积分结果%提醒：超精确的高斯点和高斯系数可由LZ.m子程序单独计算%提醒：精度要求不高时，可以自己手动输入。
GaussP,GaussA均为数组%使用范例：% 0.先输入必要参数：% a=0;%函数求积分的区间下限% b=pi；% f=@(x)exp(x)*cos(x);% n=5;% [GaussP,GaussA]=LZ(n);% 1.复化高斯勒让德积分方法，结果可以为精确解% 例如：% h=(b-a)/1000; %积分区间划分为1000份，也可以弄成一万份，一亿份。
% A=0; %预先开一个空集% for i=1:1000 %循环累加这1000份% A=A+gslrdjf(f, a+(i-1)*h,a+i*h ,5,GaussP, GaussA);% end% A %输出精确解结果% 2.高斯勒让德积分方法，结果虽然不可能是精确解，但是精度也不会差。
% 例如：% A=gslrdjf(f,a,b,n,GaussP, GaussA)%---------------------------------------------------------------------------------------p=GaussP; %获取高斯点数组Xkq=GaussA; %获取高斯点对应系数AkA=(b-a)/2; %积分区间转化为[-1,1]B=(b+a)/2;result=0; %预先开一个空集for i=1:n+1 %循环实现累加result=result+q(i)*f(A*p(i)+B); %高斯求积公式endresult=A*result; %补上积分区间转化后应有的系数end。
数值分析上机实验报告《数值分析》上机实验报告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实现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=，进行迭代。
课程设计任务书摘要实验方法与理论方法是推动科学技术发展的两大基本方法，但有局限性。
许多研究对象，由于空间或时间的限制，既不可能用理论精确描述，也不能用实验手段实现。
数值模拟或称为科学计算突破了实验和理论科学的局限，在科技发展中起到越来越重要的作用。
可以认为，科学计算已于实验、理论一起成为科学方法上不可或缺的三个主要手段。
计算数学的研究是科学计算的主要组成部分，而数值分析则是计算数学的核心。
数值计算是研究使用计算机来解决各种数学问题的近似计算方法与理论，其任务是提供在计算机上可解的、理论可靠的、计算复杂性低的各种常用算法。
数值分析的主要内容：1）、数值代数：求解线性和非线性方程组的解，分直接方法和间接方法两大类；2）、插值、曲线拟合和数值逼近；3）、数值微分和数值积分；4）、常微分和偏微分方程数值解法。
本文主要通过Matlab软件，对数值分析中的一些问题进行求解，如列主元Gauss 消去法，Lagrange插值多项式，复化Simpson公式，Runge-Kutta方法以及数值分析在实际问题中的应用，并在求解的过程中更加熟识这门课程的主要内容，以及加强对课程知识的掌握。
在学习与设计计算方法时，从数学理论角度，学会分析方法的误差、收敛性和稳定性，保证计算方法的准确性；从实际应用的角度出发，掌握计算方法的结构与流程，能够把计算方法转换为可在计算机上直接处理的程序，保证算法的可用性。
关键词：列主元Gauss消去法；Lagrange插值；复化Simpson公式；Runge-Kutta目录实验一列主元Gauss消去法 01.1 实验目的 01.2 基本原理 01.3 实验内容 (1)1.4 实验结论 (2)实验二拉格朗日插值多项式 (3)2.1 实验目的 (3)2.2 基本原理 (3)2.3 实验内容 (3)2.4 实验结论 (8)实验三复化Simpson求积公式 (9)3.1 实验目的 (9)3.2 基本原理 (9)3.3 实验内容 (9)3.4 实验结论 (11)实验四龙格-库塔(Runge-Kutta)方法 (12)4.1 实验目的 (12)4.2 基本原理 (12)4.3 实验内容 (13)4.4 实验结论 (14)实验五数值方法实际应用 (15)5.1 实验目的 (15)5.2 基本原理 (15)5.3 实验内容 (15)5.4 实验结论 (21)参考文献 (22)实验一 列主元Gauss 消去法1.1 实验目的1) 理解列主元消去法的原理；2) 熟悉列主元消去法的计算步骤，能用代码编写； 3) 解决实际问题。
Fractor是一个基于Python的开源软件，用于执行复杂的数学计算和数值分析。
它使用Python语言编写，并使用NumPy和SciPy等科学计算库进行数学运算。
以下是使用Fractor进行数学计算和数值分析的示例代码：import numpy as npfrom fractor import Fractor# 定义一个Fractor对象f = Fractor()# 定义一个数值数组x = np.linspace(-10, 10, 100)# 定义一个函数def f(x):return np.sin(x) + np.cos(x)# 使用Fractor对象计算函数的值y = f(x)# 将计算结果保存到变量y中f.result('y', y)# 使用Fractor对象进行数值积分和微分运算integral = f.integrate('y', 0, 10) # 对y在[0,10]区间进行积分derivative = f.derivative('y') # 对y进行微分# 输出结果print('Integral:', integral)print('Derivative:', derivative)在上面的代码中，我们首先导入了NumPy和Fractor库。
然后，我们定义了一个Fractor对象，并使用NumPy的linspace函数定义了一个数值数组x。
接下来，我们定义了一个函数f，该函数接受一个数值参数x，并返回sin(x) +cos(x)的值。
然后，我们使用Fractor对象f计算了函数f在数值数组x上的值，并将结果保存到变量y中。
最后，我们使用Fractor对象f对变量y进行数值积分和微分运算，并将结果输出到控制台。
需要注意的是，Fractor还提供了许多其他功能，如矩阵运算、特征值计算、曲线拟合等等。
如果您需要更多关于Fractor的详细信息和使用示例，请查阅Fractor的官方文档。
实验一：拉格朗日插值多项式命名（源程序.cpp及工作区.dsw）：lagrange问题：4//Lagrange.cpp#include <stdio.h>#include <conio.h>#define N 4int checkvalid(double x[], int n);void printLag (double x[], double y[], double varx, int n);double Lagrange(double x[], double y[], double varx, int n);void main (){double x[N+1] = {0.4, 0.55, 0.8, 0.9, 1};double y[N+1] = {0.41075, 0.57815, 0.88811, 1.02652, 1.17520};double varx = 0.5;if (checkvalid(x, N) == 1){printf("\n\n插值结果: P(%f)=%f\n", varx, Lagrange(x, y, varx, N));}else{printf("结点必须互异");}getch();}int checkvalid (double x[], int n){int i,j;for (i = 0; i < n; i++){for (j = i + 1; j < n+1; j++){if (x[i] == x[j])//若出现两个相同的结点，返回-1{return -1;}}}return 1;}double Lagrange (double x[], double y[], double varx, int n) {double fenmu;double fenzi;double result = 0;int i,j;printf("Ln(x) =\n");for (i = 0; i < n+1; i++){fenmu = 1;for (j = 0; j < n+1; j++){if (i != j){fenmu = fenmu * (x[i] - x[j]);}}printf("\t%f", y[i] / fenmu);fenzi = 1;for (j = 0; j < n+1; j++){if (i != j){printf("*(x-%f)", x[j]);fenzi = fenzi * (varx - x[j]);}}if (i != n){printf("+\n");}result += y[i] / fenmu * fenzi;}return result;}运行及结果显示：实验二：牛顿插值多项式命名（源程序.cpp及工作区.dsw）：newton_cz问题：4//newton_cz.cpp#include<stdio.h>#include<iostream.h>#include<conio.h>#define N 4int checkvaild(double x[],int n){int i,j;for(i=0;i<n+1;i++){for(j=i+1;j<n+1;j++)if(x[i]==x[j])return -1;}return 1;}void chashang(double x[],double y[],double f[][N+1]){int i,j,t=0;for(i=0;i<N+1;i++){f[i][0]=y[i];//f[0][0]=y[0],f[1][0]=y[1];f[2][0]=y[2];f[3][0]=y[3];f[4][0]=y[4] }for(j=1;j<N+1;j++)// 阶数j{for(i=0;i<N-j+1;i++) //差商个数if[i][j]=(f[i+1][j-1]-f[i][j-1])/(x[t+i+1]-x[i]);//一阶为f[0][1]~f[3][1];二阶为f[0][2]~f[2][2]//三阶为f[0][3]~f[1][3];四阶为f[0][4]t++;}}double compvalue(double t[][N+1],double x[],double varx){int i,j,r=0;double sum=t[0][0],m[N+1]={sum,1,1,1,1};printf("i\tXi\t F(Xi)\t 1阶\t 2阶\t\t3阶\t 4阶\n");printf("--------------------------------------------------------------------------------");for(i=0;i<N+1;i++){r=i;printf("x%d\t%f ",i,x[i]);for(j=0;j<=i;j++){printf("%f ",t[r][j]);r--;}printf("\n");}printf("--------------------------------------------------------------------------------");/**/ printf("N(x) =\n");printf(" %f\n",t[0][0]);for(i=1;i<N+1;i++){for(j=0;j<i;j++)m[i]=m[i]*(varx-x[j]);m[i]=t[0][i]*m[i];sum=sum+m[i];}for(i=1;i<N+1;i++){ printf(" +%f*",t[0][i]);for(j=0;j<i;j++)printf("(x-%f)",x[j]);printf("\n");}return sum;}void main(){double varx,f[N+1][N+1];double x[N+1]={0.4,0.55,0.8,0.9,1};double y[N+1]={0.41075,0.57815,0.88811,1.02652,1.17520};checkvaild(x,N);chashang(x,y,f);varx=0.5;if(checkvaild(x,N)==1){chashang(x,y,f);printf("\n\n牛顿插值结果: P(%f)=%f\n",varx,compvalue(f,x,varx));}elseprintf("输入的插值节点的x值必须互异!");getch();} 运行及结果显示：实验三：自适应梯形公式命名（源程序.cpp 及工作区.dsw ）：autotrap ][n T问题：计算⎰+=10214dx x π的近似值，使得误差（EPS ）不超过610- //autotrap.cpp #include<stdio.h> #include<conio.h> #include<math.h> #define EPS 1e-6 double f(double x) { return 4/(1+x*x); } double AutoTrap(double(*f)(double),double a,double b,double eps) { double t2,t1,sum=0.0; int i,k; t1=(b-a)*(f(a)+f(b))/2; printf("T(%4d)=%f\n",1,t1); for(k=1;k<11;k++) { for(i=1;i<=pow(2,k-1);i++) sum+=f(a+(2*i-1)*(b-a)/pow(2,k)); t2=t1/2+(b-a)*sum/pow(2,k); printf("T(%4d)=%f\n",(int)pow(2,k),t2); if(fabs(t2-t1)<EPS) break; else t1=t2; sum=0.0; } return sum; } void main(){ double s;s=AutoTrap(f,0.0,1.0,EPS);getch();} 运行及结果显示：实验四：龙贝格算法命名（源程序.cpp 及工作区.dsw ）：romberg问题：求⎰+=10214dx x π的近似值，要求误差不超过610-//romberg.cpp #include <stdio.h> #include <conio.h> #include <math.h> #define N 20 #define EPS 1e-6 double f(double x){return 4/(1+x*x);} double Romberg(double a,double b,double (*f)(double),double eps) { double T[N][N],p,h;int k=1,i,j,m=0;T[0][0]=(b-a)/2*(f(a)+f(b));do{p=0;h=(b-a)/pow(2,k-1);for(i=1;i<=pow(2,k-1);i++)p=p+f(a+(2*i-1)*h/2);T[0][k]=T[0][k-1]/2+p*h/2;for(i=1;i<=k;i++){j=k-i;T[i][j]=(pow(4,i)*T[i-1][j+1]-T[i-1][j])/(pow(4,i)-1); }k++; p=fabs(T[k-1][0]-T[k-2][0]);}while(p>=EPS);k--; while(m<=k){for(i=0;i<=m;i++) printf("T(%d,%d)=%f ",i,m-i,T[i][m-i]);m++;printf("\n");}return T[k][0]; } void main() {printf("\n\n 积分结果 = %f",Romberg(0,1,f,EPS)); getch(); } 运行及结果显示：实验五：牛顿切线法问题：求方程01)(3=--=x x x f 在5.1=x ，6.0=x 附近的根（精度31021-⨯=） //newton_qxf.cpp#include <math.h>#include<conio.h>#include <stdio.h>#define MaxK 1000#define EPS 0.5e-3double f(double x){return x*x*x-x-1;}double f1(double x){return 3*x*x-1;}int newton(double (*f)(double), double (*f1)(double), double &x0, double eps) {double xk, xk1;int count = 0;printf("k\txk\n");printf("-----------------------\n");xk = x0;printf("%d\t%f\n", count, xk);do{if((*f1)(xk)==0)return 2;xk1 = xk - (*f)(xk) / (*f1)(xk);if (fabs(xk - xk1) < eps){count++;xk = xk1;printf("%d\t%f\n", count, xk);x0 = xk1;return 1;}count++;xk = xk1;printf("%d\t%f\n", count, xk);}while(count < MaxK);return 0;}void main(){for(int i=0;i<2;i++){double x=0.6;if(i==1)x=1.5;printf("x0初值为%f:\n",x);if (newton(f, f1, x, EPS) == 1){printf("-----------------------\n");printf("the root is x=%f\n\n\n", x);}else{printf("the method is fail!");}}getch();}实验六：牛顿下山法命名（源程序.cpp 及工作区.dsw ）：newton_downhill 问题：求方程01)(3=--=x x x f 在5.1=x ，6.0=x 附近的根（精度31021-⨯=） //newton_downhill.cpp #include <stdio.h> #include <conio.h> #include <math.h> #include <stdlib.h> #define Et 1e-3//下山因子下界 #define E1 1e-3//根的误差限 #define E2 1e-3//残量精度 double f(double x) { return x * x * x - x - 1; } double f1(double x) { return 3 * x * x - 1; } void errormess(int b){char *mess;switch(b){case -1:mess = "f(x)的导数为0!";break;case -2:mess = "下山因子已越界，下山处理失败";break;default:mess = "其他类型错误!";}printf("the method has faild! because %s", mess);}int Newton(double (*f)(double), double (*f1)(double), double &x0) {int k = 0;double t;double xk, xk1;double fxk, fxk_temp;printf("k t xk f(xk)\n");printf("----------------------------------------------------------\n");printf("%-20d", k);xk = x0;printf("%-15f", x0);fxk = (*f)(xk);printf("%-20f", fxk);printf("\n");for (k = 1; ; k++){t = 1;while(1){printf("%-10d", k);printf("%-10f", t);if((*f1)(xk) != 0){xk1 = xk - t * (*f)(xk) / (*f1)(xk);}else{return -1;}printf("%-15f", xk1);fxk_temp = (*f)(xk1);printf("%-20f", fxk_temp);if(fabs(fxk_temp) >fabs(fxk)){t = t / 2;printf("\n");if (t < Et){return -2;}}else{printf("下山成功\n");break;}}if (fabs(xk-xk1)<E1){x0 = xk1;return 1; } xk = xk1; } }void main() {int b;double x0; x0 = 0.6;b = Newton(f, f1, x0); if (b == 1) printf("\nthe root x=%f\n", x0); else errormess(b); getch(); }运行及结果显示：实验七：埃特金加速算法命名（源程序.cpp 及工作区.dsw ）：aitken问题：求方程01)(3=--=x x x f 在5.1=x ，6.0=x 附近的根（精度31021-⨯=） //aitken.cpp#include <math.h> #include <stdio.h> #include <conio.h> #define MaxK 100 #define EPS 0.5e-3double g(double x) {return x * x * x - 1; }int aitken (double (*g)(double), double &x, double eps) {int k;double xk = x, yk, zk, xk1;printf("k xk yk zk xk+1\n");printf("-------------------------------------------------------------------\n"); for (k = 0;k<MaxK; k++) { yk = (*g)(xk); zk = (*g)(yk);xk1 = xk - (yk - xk) * (yk - xk) / (zk - 2 * yk + xk); printf("%-10d%-15f%-15f%-15f%-15f\n", k, xk, yk, zk, xk1); if (fabs(yk-xk)<=eps) { x = xk1; return k + 1; } xk = xk1; }return -1; }void main () {double x = 1.5; int k;k = aitken(g, x, EPS); if (k == -1) printf("迭代次数越界!\n"); else printf("\n 经k=%d 次迭代，所得方程根为:x=%f\n", k, x); getch(); }运行及结果显示：实验八：正割法问题：求方程01)(3=--=x x x f 在5.1=x 附近的根（精度0.5e-8） //ZhengGe.cpp #include <math.h> #include <stdio.h> #include <conio.h>#define MaxK 1000 #define EPS 0.5e-8double f(double x) {return x*x*x-x-1;}int ZhengGe(double (*f)(double), double x0, double &x1, double eps) {printf("k xk f(xk)\n");printf("---------------------------------------------\n");int k;double xk0, xk, xk1;xk0 = x0;printf("%-10d%-15f%-15f\n", 0, x0, (*f)(x0));xk = x1;for (k=1;k<MaxK;k++){if((*f)(xk)-(*f)(xk0)==0)return -1;xk1 = xk - (*f)(xk) / ( (*f)(xk) - (*f)(xk0) ) * (xk - xk0);printf("%-10d%-15f%-15f\n", k, xk, (*f)(xk));if (fabs(xk1 - xk)<=eps){printf("%-10d%-15f%-15f\n\n", k + 1, xk1, (*f)(xk1));printf("---------------------------------------------\n\n");x1 = xk1;return 1;}xk0 = xk;xk = xk1;}return -2;}void main (){double x0 = 1, x1 = 2;if (ZhengGe(f, x0, x1, EPS) == 1){printf("the root is x = %f\n", x1);}else{printf("the method is fail!");}getch();}实验九：高斯列选主元算法命名（源程序.cpp 及工作区.dsw ）：colpivot问题：求解方程组并计算系数矩阵行列式值 ⎪⎩⎪⎨⎧=-+=+-=-+2240532321321321x x x x x x x x x//colpivot.cpp#include <math.h> #include <stdio.h> #include <conio.h> #define N 3static double aa[N][N]={{1,2,-1},{1,-1,5},{4,1,-2}}; static double bb[N+1]={3,0,2};void main() {int i,j;double a[N+1][N+1],b[N+1],x[N+1],det;double gaussl(double a[][N+1],double b[],double x[]); for(i=1;i<=N;i++) { for(j=1;j<=N;j++) a[i][j]=aa[i-1][j-1]; b[i]=bb[i-1]; }det=gaussl(a,b,x); if(det!=0) { printf("\n 方程组的解为:"); for(i=1;i<=N;i++) printf(" x[%d]=%f",i,x[i]); printf("\n\n 系数矩阵的行列式值＝%f",det); }else printf("\n\n 系数矩阵奇异阵，高斯方法失败 !:"); getch(); }double gaussl(double a[][N+1],double b[],double x[])//a传入增广矩阵(有效元素为a[1,1]...a[n,n+1]),x负责传入迭代初值并返回解向量//(有效值为x[1]...x[n]);返回值：系数矩阵行列式值detA{double det=1.0,F,m,temp;int i,j,k,r;void disp(double a[][N+1],double b[]);printf("消元前增广矩阵:\n");disp(a,b);for(k=1;k<N;k++){temp=a[k][k];r=k;for(i=k+1;i<=N;i++){if(fabs(temp)<fabs(a[i][k])){temp=a[i][k];r=i;}//按列选主元，即确定ik}if(a[r][k]==0)return 0;//如果aik,k=0,则A为奇异矩阵，停止计算if(r!=k){for(j=k;j<=N;j++){a[k][j]+=a[r][j];a[r][j]=a[k][j]-a[r][j];a[k][j]-=a[r][j];}b[k]+=b[r];b[r]=b[k]-b[r];b[k]-=b[r];det=-det;//如果ik!=k,则交换[A,b]第ik行与第k行元素printf("交换第%d, %d行:\n",k,r);disp(a,b);}for(i=k+1;i<=N;i++){m=a[i][k]/a[k][k];for(j=1;j<=k;j++)a[i][j]=0.0;for(j=k+1;j<=N;j++)a[i][j]-=m*a[k][j];b[i]-=m*b[k];}printf("第%d次消元:\n",k);disp(a,b);det*=a[k][k];} //大FOR循环结束x[N]=b[N]/a[N][N];for(i=N-1;i>0;i--){F=0.0;for(j=i+1;j<=N;j++)F+=a[i][j]*x[j];x[i]=(b[i]-F)/a[i][i];}det*=a[N][N];return det;}void disp(double a[][N+1],double b[])//显示选主元及消元运算中间增广矩阵结果 {int i,j;for(i=1;i<=N;i++) {for(j=1;j<=N;j++)printf("%10f\t",a[i][j]); printf("%10f\n",b[i]); } }运行及结果显示：实验十：高斯全主元消去算法 命名（源程序.cpp 及工作区.dsw ）：fullpivot问题：利用全主元消去法求解方程组⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----9555.04525.015.0201.0321x x x //fullpivot.cpp#include <math.h> #include <stdio.h> #include <conio.h> #define N 3 double x[N+1];static double aa[N][N]={{0.01,2,-0.5},{-1,-0.5,2},{5,-4,0.5}}; static double bb[N]={-5,5,9};void main(){int i,j;double a[N+1][N+1],b[N+1],x[N+1],det;int t[N+1];//引入列交换保存x向量的各分量的位置顺序double gaussl(double a[][N+1],double b[],double x[],int t[]);for(i=1;i<=N;i++){for(j=1;j<=N;j++) a[i][j]=aa[i-1][j-1];b[i]=bb[i-1];}for(i=1;i<=N;i++)t[i]=i;det=gaussl(a,b,x,t);if(det!=0){ printf("\n方程组的解为:");for(i=N;i>=1;i--){printf(" x[%d]=%f",t[i],x[i]);if(i>1)printf(" -->");}printf("\n\n系数矩阵的行列式值＝%f",det);}else printf("\n\n系数矩阵奇异阵，高斯方法失败!:");getch();}double gaussl(double a[][N+1],double b[],double x[],int t[])//a传入增广矩阵(有效元素为a[1,1]...a[n,n+1]),x负责传入迭代初值并返回解向量//(有效值为x[1]...x[n]);返回值：系数矩阵行列式值detA{double det=1.0,F,m,temp;int i,j,k,r,s;void disp(double a[][N+1],double b[],int x[]);// printf("消元前增广矩阵:\n");//disp(a,b,t);for(k=1;k<N;k++){temp=a[k][k];r=k;s=k;for(i=k;i<=N;i++)for(j=k;j<=N;j++)if(fabs(temp)<fabs(a[i][j])){temp=a[i][j];r=i;s=j;}//选主元，选取ik,jkif(a[r][s]==0)return 0;if(r!=k){for(j=k;j<=N;j++){a[k][j]+=a[r][j];a[r][j]=a[k][j]-a[r][j];a[k][j]-=a[r][j];}b[k]+=b[r];b[r]=b[k]-b[r];b[k]-=b[r];det=-det;printf("交换第%d, %d行:\n",k,r);disp(a,b,t);}if(s!=k){for(i=0;i<=N;i++){a[i][k]+=a[i][s];a[i][s]=a[i][k]-a[i][s];a[i][k]-=a[i][s];}t[k]+=t[s];t[s]=t[k]=t[s];t[k]-=t[s];det=-det;printf("交换第%d, %d列:\n",k,s);disp(a,b,t);}for(i=k+1;i<=N;i++){m=a[i][k]/a[k][k];for(j=1;j<=k;j++)a[i][j]=0.0;for(j=k+1;j<=N;j++)a[i][j]-=m*a[k][j];b[i]-=m*b[k];} //消元计算printf("第%d次消元:\n",k);disp(a,b,t);det*=a[k][k];} //大FOR循环结束x[N]=b[N]/a[N][N];for(i=N-1;i>0;i--){F=0.0;for(j=i+1;j<=N;j++)F+=a[i][j]*x[j];x[i]=(b[i]-F)/a[i][i];}//回代计算det*=a[N][N];return det;}void disp(double a[][N+1],double b[],int x[])//显示选主元及消元运算中间增广矩阵结果{int i,j;for(i=1;i<=N;i++){for(j=1;j<=N;j++)printf("%f\t",a[i][j]);printf("| x%d = %f\n",x[i],b[i]);}}运行及结果显示：实验十一：LU 分解算法命名（源程序.cpp 及工作区.dsw ）：LU问题：利用LU 法求解方程组 ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡472111312613321x x x //LU.cpp#include<math.h> #include<stdio.h> #include<conio.h> #define N 3static aa[N][N]={{3,1,6},{2,1,3},{1,1,1}}; static bb[N]={2,7,4};void main() {int i,j;double a[N+1][N+1],b[N+1];void solve(double a[][N+1],double b[N+1]); int LU(double a[][N+1]); for(i=1;i<=N;i++) { for(j=1;j<=N;j++) a[i][j]=aa[i-1][j-1]; b[i]=bb[i-1]; }if(LU(a)==1){printf("矩阵L如下\n");for(i=1;i<=N;i++){for(j=1;j<=i;j++)if(i==j)printf("%12d ",1);else printf("%12f",a[i][j]);printf("\n");}printf("\n矩阵U如下\n");for(i=1;i<=N;i++){for(j=1;j<=N;j++)if(i<=j)printf("%12f",a[i][j]);else printf("%12s","");printf("\n");}solve(a,b);for(i=1;i<=N;i++)printf("x[%d]=%f ",i,b[i]);printf("\n");}else printf("\nthe LU method failed!\n");getch();}int LU(double a[][N+1])//对N阶矩A阵进行LU分解，结果L、U存放在A的相应位置{int i,j,k,s;double m,n;for(i=2;i<=N;i++)a[i][1]/=a[1][1];for(k=2;k<=N;k++){for(j=k;j<=N;j++){m=0;for(s=1;s<k;s++)m+=a[k][s]*a[s][j];a[k][j]-=m;}if(a[k][k]==0){printf("a[%d][%d]=%d ",k,k,0);return -1;//存在元素akk为0}for(i=k+1;i<=N;i++){n=0;for(s=1;s<k;s++)n+=a[i][s]*a[s][k];a[i][k]=(a[i][k]-n)/a[k][k];}}return 1;//正常结束}void solve(double a[][N+1],double b[N+1])//利用分解的LU求x//回代求解，L和U元素均在A矩阵相应位置；b存放常数列，返回时存放方程组的解{double y[N+1],F;int i,j;y[1]=b[1];for(i=2;i<=N;i++){F=0.0;for(j=1;j<i;j++)F+=a[i][j]*y[j];y[i]=b[i]-F;}b[N]=y[N]/a[N][N];for(i=N-1;i>0;i--){F=0.0;for(j=N;j>i;j--)F+=a[i][j]*b[j];b[i]=(y[i]-F)/a[i][i];}}运行及结果显示：实验十二：Guass-Sediel 算法命名（源程序.cpp 及工作区.dsw ）：GS问题：利用G －S 法求解方程组 ⎪⎩⎪⎨⎧=+--=-+-=--2.453.82102.7210321321321x x x x x x x x x （精度为41021-⨯）//Guass_sediel.cpp#include "iostream.h"#include "math.h"#include "stdio.h"#include<conio.h>#define N 3 //方程的阶数#define MaxK 100 //最大迭代次数#define EPS 0.5e-4 //精度控制static double aa[N][N]={{10,-1,-2},{-1,10,-2},{-1,-1,5}};static double bb[N]={7.2,8.3,4.2};void main(){int i,j;double x[N+1];double a[N+1][N+1],b[N+1];int GaussSeidel(double a[][N+1],double b[],double x[]);for(i=1;i<=N;i++){for(j=1;j<=N;j++)a[i][j]=aa[i-1][j-1];b[i]=bb[i-1];x[i]=0;}if(GaussSeidel(a,b,x)==1){printf("the roots is:");for(i=1;i<=N;i++)printf(" x[%d]=%f ",i,x[i]);printf("\n");}else printf("\nthe G-S method failed!\n");getch();}int GaussSeidel(double a[][N+1],double b[],double x[])//a传入系数矩阵(有效元素为a[1,1]...a[n,n]),b为方程组右边常数列，x传入迭代初值并返回解向量//x**(k+#x)=Gx**(k)+f{int k=1,i,j;double m,max,t[N+1];while(true){printf("k=%d:",k);max=0.0;for(i=1;i<=N;i++){if(a[i][i]==0)return -1;//存在元素akk为0m=0.0;t[i]=x[i];for(j=1;j<=N;j++)if(j!=i)m+=a[i][j]*x[j];x[i]=(b[i]-m)/a[i][i];if(max<fabs(x[i]-t[i]))max=fabs(x[i]-t[i]);printf(" x[%d]=%f ",i,x[i]);}printf("\n");if(max<EPS)return 1;//正常结束elsek++;if(k>MaxK)return -2;//迭代次数越界}}运行及结果显示：实验十三：SOR 算法命名（源程序.cpp 及工作区.dsw ）：sor问题：利用SOR 法求解方程组⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----111141111411114111144321x x x x （精度为51021-⨯） //SOR.cpp#include "math.h"#include "stdio.h"#include "conio.h"#define N 4#define MaxK 1000#define EPS 0.5e-5static double aa[N][N]={{-4,1,1,1},{1,-4,1,1},{1,1,-4,1},{1,1,1,-4}}; static double bb[N]={1,1,1,1};void main(){int i,j;double x[N+1];double a[N+1][N+1],b[N+1];int Sor(double a[][N+1],double b[],double w,double x[]);for(i=1;i<=N;i++){for(j=1;j<=N;j++)a[i][j]=aa[i-1][j-1];b[i]=bb[i-1];x[i]=0;}if(Sor(a,b,1.3,x)==1){printf("the roots is:");for(i=1;i<=N;i++)printf(" x[%d]=%f ",i,x[i]);printf("\n");}else printf("\nthe G-S method failed!\n");getch();}int Sor(double a[][N+1],double b[],double w,double x[]) {int i,j,k;double loft,comb;for(k=0;k<=N;k++)x[k]=0;for(i=1;i<=MaxK && fabs(comb)>EPS;i++){printf("k=%d:\t",i);comb=0;for(k=1;k<=N;k++){loft=b[k];for(j=1;j<N+1;j++)loft-=a[k][j]*x[j];if(a[k][k]==0)return -1;else loft=w*loft/a[k][k];x[k]=x[k]+loft;printf("x[%d]=%10f \t",k,x[k]);if(fabs(loft)>comb)comb=fabs(loft);}printf("\n");}if(i>MaxK)return -2;elsereturn 1;}运行结果（w=1）=。
