计算方法上机作业
计算方法上机报告_董晓壮
计算方法A上机报告学院（系）：电气工程学院学生姓名：陶然学号：授课老师：完成日期：2019年12月03日西安交通大学Xi'an Jiaotong University目录1 QR分解法求解线性方程组 (2)1.1 算法原理 (2)1.1.1 基于吉文斯变换的QR分解 (2)1.1.2 基于豪斯霍尔德变换的QR分解 (3)1.2 程序流程图 (4)1.2.1 基于吉文斯变换的QR分解流程图 (4)1.2.2 基于豪斯霍尔德变换的QR分解流程图 (5)1.3 程序使用说明 (5)1.3.1 基于吉文斯变换的QR分解程序说明 (5)1.3.2 基于豪斯霍尔德变换的QR分解程序说明 (7)1.4 算例计算结果 (8)2 共轭梯度法求解线性方程组 (10)2.1 算法原理 (10)2.2 程序流程图 (10)2.3 程序使用说明 (11)2.4 算例计算结果 (12)3 三次样条插值 (14)3.1 算法原理 (14)3.2 程序流程图 (16)3.3 程序使用说明 (17)3.4 算例计算结果 (19)4 龙贝格积分 (21)4.1 算法原理 (21)4.2 程序流程图 (22)4.3 程序使用说明 (23)4.4 算例计算结果 (24)结论 (26)1 QR 分解法求解线性方程组1.1 算法原理矩阵的QR 分解是指，可以将矩阵A 分解为一个正交矩阵Q 和一个上三角矩阵R 的乘积，实际中，QR 分解经常被用来解决线性最小二乘问题，分解情况如图1.1所示。
=⨯图1.1 QR 分解示意图本次上机学习主要进行了两个最基本的正交变换—吉文斯（Givens ）变换和豪斯霍尔德（Householder ）变换，并由此导出矩阵的QR 分解以及求解线性方程组的的方法。
1.1.1 基于吉文斯变换的QR 分解吉文斯矩阵也称初等旋转阵，如式（1.1）所示，它把n 阶单位矩阵I 的第,i j 行的对角元改为c ，将第i 行第j 列的元素改为s ，第j 行第i 列的元素改为s −后形成的矩阵。
计算方法大作业——龙贝格积分
3 龙贝格积分3.1 算法原理及程序框图龙贝格积分法是在复化梯形求积公式、复化辛普森求积公式和复化科茨求积公式关系的基础上，构造出的一种精度更高的数值积分方法。
对于复化梯形求积公式而言，近似积分为()2221[]41n n n n I f T T T T ≈+-=-.(11) 对于复化辛普森求积公式和复化科茨求积公式而言，也有类似的关系，如公式(12)和公式(13)。
()22221[]41n n n n I f S S S S ≈+-=- (12)()22231[]41n n n n I f C C C C ≈+-=- (13)通过对公式(11)~(13)做进一步分析，可得到公式(14)和公式(15)。
()22141n n n n S T T T =+--(14)()222141n n n n C S S S =+-- (15)根据公式(14)和公式(15)表现出来的规律，令龙贝格积分为()223141n n n n R C C C =+-- (16)其截断误差为c R h 8f (8)(η)，已经具有很高的精度。
龙贝格积分法是将区间[a , b ]逐次分半进行计算，因此，对已知函数f (x )在区间[a , b ]上的龙贝格积分法的计算公式的算法如下，程序框图如图13所示。
(1) 计算T 1：[]1()()2b aT f a f b -=+；(2) 逐次计算T 2k +1：()1211221121,0,1,2222kk k k k i b a b a T T f a i k +++=--⎛⎫=++-= ⎪⎝⎭∑；(3) 逐次计算S 2k 、C 2k 和R 2k ：()()()11111122222222232222141141141kk k k k k k k k k k k S T T T C S S S R C C C ++++++⎧=+-⎪-⎪⎪=+-⎨-⎪⎪=+-⎪-⎩；(4) 若122k k R R ε+-<，则取[]12k I f R +≈；否则，继续计算，直到满足精度为止。
计算方法上机作业——龙格库塔法
计算方法上机作业——龙格库塔法龙格库塔法（Runge-Kutta method）是一种常用于求解常微分方程（Ordinary Differential Equation，ODE）的数值解法。
它是由德国数学家卡尔·龙格（Carl Runge）和马丁·威尔海姆·库塔（Martin Wilhelm Kutta）在20世纪初提出的。
龙格库塔法的基本思想是通过数值逼近来计算微分方程的近似解。
在讲解龙格库塔法之前，我们先来简单回顾一下ODE的一阶常微分方程的基本形式：y′(y)=y(y,y)其中，y(y,y)是已知函数。
龙格库塔法的核心是使用差分逼近计算函数的斜率。
假设我们要求解的方程为：y′(y)=y(y,y),y(y)=y₀所需计算的点为y₀,y₁,...,yy，对应的函数值为y₀,y₁,...,yy，其中y是步长的个数。
龙格库塔法通过递推关系式来计算估计值，并不断更新当前点的函数值。
接下来以龙格库塔法的经典四阶形式为例进行说明。
该方法的基本方程如下：yy+1=yy+(y₁+2y₂+2y₃+y₄)/6y₁=ℎy(yy,yy)y₂=ℎy(yy+ℎ/2,yy+y₁/2)y₃=ℎy(yy+ℎ/2,yy+y₂/2)y₄=ℎy(yy+ℎ,yy+y₃)其中y表示当前步骤，ℎ表示步长，yy表示当前点的函数值，y₁,y₂,y₃和y₄则表示对应的斜率。
使用龙格库塔法，我们可以通过不断递归计算来求得指定区间（例如[y,y]）上的函数值。
具体步骤如下：1.确定求解区间[y,y]和初始点(y₀,y₀)以及步长ℎ。
2.初始化：设置yy=y₀，yy=y₀。
3.对所有y=0,...,y−1：计算y₁,y₂,y₃和y₄，根据上述递推关系式。
根据递推关系式计算yy+1更新当前点的函数值，即yy+1=y(yy+1)。
更新当前点的y值，即yy+1=yy+ℎ。
4.返回结果：最终求得的函数值。
需要注意的是，选择适当的步长对最终结果的精度和计算效率都有重要影响。
数值计算方法上机实验报告
数值计算方法上机实验报告
一、实验目的
本次实验的主要目的是熟悉和掌握数值计算方法，学习梯度下降法的
原理和实际应用，熟悉Python语言的编程基础知识，掌握Python语言的
基本语法。
二、设计思路
本次实验主要使用的python语言，利用python下的numpy，matplotlib这两个工具，来实现数值计算和可视化的任务。
1. 首先了解numpy的基本使用方法，学习numpy的矩阵操作，以及numpy提供的常见算法，如矩阵分解、特征值分解等。
2. 在了解numpy的基本操作后，可以学习matplotlib库中的可视化
技术，掌握如何将生成的数据以图表的形式展示出来。
3. 接下来就是要学习梯度下降法，首先了解梯度下降法的主要原理，以及具体的实际应用，用python实现梯度下降法给出的算法框架，最终
可以达到所期望的优化结果。
三、实验步骤
1. 熟悉Python语言的基本语法。
首先是熟悉Python语言的基本语法，学习如何使用Python实现变量
定义，控制语句，函数定义，类使用，以及面向对象编程的基本概念。
2. 学习numpy库的使用方法。
其次是学习numpy库的使用方法，学习如何使用numpy库构建矩阵，学习numpy库的向量，矩阵操作，以及numpy库提供的常见算法，如矩阵分解，特征值分解等。
3. 学习matplotlib库的使用方法。
西安电子科技大学出版社计算方法上机答案
西安电子科技大学出版社《计算方法》任传祥等编著第九章计算方法上机参考答案实验一，算法一#include <stdio.h>#include <math.h>double I0=log(6)/log(5),I1;int n=1;main (){while(1){I1=1.0/(n)-I0*5.0;printf("%d %lf\n", n,I1);if(n>=20)break;elseI0=I1;n++;}}实验一，算法二#include <stdio.h>#include <math.h>double I0=(1/105.0+1/126.0)/2,I1;int n=20;main (){printf("%d %lf\n", n,I0);while(1){I1=1.0/(5.0*n)-I0/5.0;printf("%d %lf\n", n-1,I1);if(n<2)break;elseI0=I1;n--;}}实验二，二分法#include <stdio.h>#include <math.h>#define esp 1e-3double f(double x);main (){double a=1,b=2,x;while(fabs(b-a)>esp){x=(a+b)/2;printf("x=%lf\n",x);if(f(x)==0)break;elseif(f(x)*f(a)<0)b=x;elsea=x;}}double f(double x){return pow(x,3)-x-1;}实验二，牛顿迭代法#include<stdio.h>#include<math.h>double f(double x);double f1(double x);#define esp 1e-3void main(){double x0 = 1.5, x1;x1 = x0 - f(x0) / f1(x0);printf("x=%lf\n", x1);x0 = x1;x1 = x0 - f(x0) / f1(x0);printf("x=%lf\n", x1);while (fabs(x1 - x0)>esp){x0 = x1;x1 = x0 - f(x0) / f1(x0);printf("x=%lf\n", x1);} }double f(double x){return pow(x, 3) - x - 1;} double f1(double x){return 3 * x*x - 1;}弦割法#include<stdio.h>#include<math.h>double f(double x);#define esp 1e-3void main(){double x0 = 1.5, x1=2.0,x2;do{ x2=x1 - (x1-x0)*f(x1) /(f(x1)-f(x0));x0=x1;x1=x2;printf("x=%lf\n", x1);}while (fabs(x1 - x0)>esp);{printf("x=%lf\n", x1);}}double f(double x){return pow(x, 3) - x - 1;}实验3#include <stdio.h>/*列主元高斯消去法*/#include <math.h>float x[3],temp,max;float A[3][4]={10,-2,-1,3,-2,10,-1,15,-1,-2,5,10},c[3][4]={10,-2,-1,3,-2,10,-1,15,-1,-2,5,10}; int n=3,i,k,j,m;void main(){for(i=0;i<n;i++){max=A[i][i];k=i;for(j=j+1;j<n;j++){{max=fabs(A[j][i]);k=j;}}if(k!=i){for(j=i+1;j<=n;j++){temp=A[i][j];A[i][j]=A[k][j];A[k][j]=temp;}}for(j=i+1;j<n;j++)for(m=i+1;m<=n;m++){c[j][m]=c[j][m]+(-c[j][i]/c[i][i])*c[i][m];}}for(i=n-1;i>=0;i--){temp=0.0;for(j=n-1;j>=i+1;j--)temp=temp+c[i][j]*x[j];x[i]=(c[i][n]-temp)/c[i][i];}printf("x[1]=%f\nx[2]=%f\nx[3]=%f\n",x[0],x[1],x[2]);实验四，拉格朗日插值#include<stdio.h>int n=5,i,j;double l,L=0,X=0.5;main(){double x[5]={0.4,0.55,0.65,0.8,0.9};doubley[5]={0.41075,0.57815,0.69675,0.88811,1.02652}; for(i=0;i<n;i++){l=y[i];for(j=0;j<n;j++){if(j!=i)l=l*(X-x[j])/(x[i]-x[j]); } L=L+l;}printf("%lf\n",L);return 0;} X=0.5 X=0.7 X=0.85牛顿插值法#include<stdio.h>#include<math.h>main(){double x[5]={0.4,0.55,0.65,0.8,0.9};doubley[5]={0.41075,0.57815,0.69675,0.88811,1.02652};int n=5,i,j;double z;printf("input z\n");scanf("%lf",&z);double a[5][5];for(i=0;i<5;i++)a[i][0]=y[i];for(i=1;i<5;i++)for(j=i;j<5;j++)a[j][i]=(a[j][i-1]-a[j-1][i-1])/(x[j]-x[j-i]);double N=a[0][0],temp=1.0;for(i=1;i<n;i++){temp=temp*(z-x[i-1]);N=N+a[i][i]*temp;}printf("N=%lf\n",N);return 0;}实验五曲线拟合#include <stdio.h>#include <math.h>float x[5]={1,2,3,4,5};float y[5]={7,11,17,27,40};float A[2][3],c[2][3];float z[2],temp,max;int i,j,k,m;int n=2;void main(){for(i=0;i<5;i++){c[0][0]=A[0][0]+=1;c[0][1]=A[0][1]+=x[i];c[0][2]=A[0][2]+=y[i];c[1][0]=A[1][0]+=x[i];c[1][1]=A[1][1]+=x[i]*x[i];c[1][2]=A[1][2]+=x[i]*y[i];}/* for(i=0;i<2;i++){printf(" %lf %lf %lf\n",A[i][0],A[i][1],A[i ][2]);}*/for(i=0;i<n;i++){max=A[i][i];k=i;for(j=j+1;j<n;j++){if(fabs(A[j][i])>max){max=fabs(A[j][i]);k=j;}} if(k!=i){for(j=i+1;j<=n;j++){temp=A[i][j];A[i][j]=A[k][j];A[k][j]=temp;}}for(j=i+1;j<n;j++)for(m=i+1;m<=n;m++){c[j][m]=c[j][m]+(-c[j][i]/c[i][i])*c[i][m];}}for(i=n-1;i>=0;i--){temp=0.0;for(j=n-1;j>=i+1;j--)temp=temp+c[i][j]*z[j];z[i]=(c[i][n]-temp)/c[i][i];}printf("a=%f\nxb=%f\n",z[0],z[1]); }实验六数值积分/*梯形*/#include<stdio.h>#include<math.h> double f(double x); main(){double x[10],y[10];double h,b=1,a=0,I;int n,i;printf("n\n");scanf("%d",&n);h=(b-a)/n;for(i=0;i<=n;i++){x[i]=a+(i*h);y[i]=f(x[i]);}I=f(a)+f(b);for(i=1;i<=n-1;i++){I=I+2*y[i];}I=(h/2)*I;printf("%lf",I);}double f(double x){double f;f=1.0/(1.0+(x*x));return(f);}/*辛普森*/#include<stdio.h>#include<math.h>double f(double x);main(){double x[30],y[30];double h,b=1,a=0,I;int n,i;printf("n\n");scanf("%d",&n);//点乘2扩展h=(b-a)/n;x[10]=1;y[10]=f(x[10]);for(i=0;i<n;i++){x[2*i]=a+(i*h);y[2*i]=f(x[2*i]);x[2*i+1]=a+(i+(1.0/2.0))*h;y[(2*i)+1]=f(x[(2*i)+1]);}I=f(a)+f(b);for(i=0;i<n;i++){I=I+4*y[(2*i)+1];}for(i=1;i<n;i++){I=I+2*y[2*i];}I=(h/6)*I;printf("%lf\n",I);}double f(double x){double f;f=1.0/(1.0+(x*x));return(f);}/*梯形*//*辛普森*/。
东南大学计算方法与实习上机实验一
东南大学计算方法与实习实验报告学院：电子科学与工程学院学号：06A*****姓名：***指导老师：***实习题14、设S N=Σ(1)编制按从大到小的顺序计算S N的程序；(2)编制按从小到大的顺序计算S N的程序；(3)按两种顺序分别计算S1000，S10000，S30000，并指出有效位数。
解析：从大到小时，将S N分解成S N-1=S N-,在计算时根据想要得到的值取合适的最大的值作为首项；同理从小到大时，将S N=S N-1+ ，则取S2=1/3。
则所得式子即为该算法的通项公式。
(1)从大到小算法的C++程序如下：/*从大到小的算法*/#include<iostream>#include<iomanip>#include<cmath>using namespace std;const int max=34000; //根据第(3)问的问题，我选择了最大数为34000作为初值void main(){int num;char jus;double cor,sub;A: cout<<"请输入你想计算的值"<<'\t';cin>>num;double smax=1.0/2.0*(3.0/2.0-1.0/max-1.0/(max+1)),temps;double S[max];// cout<<"s["<<max<<"]="<<setprecision(20)<<smax<<'\n';for(int n=max;n>num;){temps=smax;S[n]=temps;n--;smax=smax-1.0/((n+1)*(n+1)-1.0);}cor=1.0/2.0*(3.0/2.0-1.0/num-1.0/(num+1.0)); //利用已知精确值公式计算精确值sub=fabs(cor-smax); //double型取误差的绝对值cout<<"用递推公式算出来的s["<<n<<"]="<<setprecision(20)<<smax<<'\n';cout<<"实际精确值为"<<setprecision(20)<<cor<<'\n';cout<<"则误差为"<<setprecision(20)<<sub<<'\n';cout<<"是否继续计算S[N],是请输入Y，否则输入N！"<<endl;cin>>jus;if ((int)jus==89||(int)jus==121) goto A;}(2)从小到大算法的C++程序如下：/*从小到大的算法*/#include<iostream>#include<iomanip>#include<cmath>using namespace std;void main(){int max;A: cout<<"请输入你想计算的数，注意不要小于2"<<'\t';cin>>max;double s2=1.0/3.0,temps,cor,sub;char jus;double S[100000];for(int j=2;j<max;){temps=s2;S[j]=temps;j++;s2+=1.0/(j*j-1.0);}cor=1.0/2.0*(3.0/2.0-1.0/j-1.0/(j+1.0)); //利用已知精确值公式计算精确值sub=fabs(cor-s2); //double型取误差的绝对值cout<<"用递推公式算出来的s["<<j<<"]="<<setprecision(20)<<s2<<'\n';cout<<"实际精确值为"<<setprecision(20)<<cor<<'\n';cout<<"则误差为"<<setprecision(20)<<sub<<'\n';cout<<"是否继续计算S[N],是请输入Y，否则输入N！"<<endl;cin>>jus;if ((int)jus==89||(int)jus==121) goto A;}(3)(注：因为程序中setprecision(20)表示输出数值小数位数20，则程序运行时所得到的有效数字在17位左右)ii.选择从小到大的顺序计算S1000、S10000、S30000的值需要计算的项S1000S10000S30000计算值0.74900049950049996 0.74966672220370571 0.74996666722220728实际精确值0.74900049950049952 0.74990000499950005 0.74996666722220373误差 4.4408920985006262*10-16 5.6621374255882984*10-15 3.5527136788005009*10-15有效数字17 17 17附上部分程序运行图：iii.实验分析通过C++程序进行计算验证采用从大到小或者从小到大的递推公式算法得到的数值基本稳定且误差不大。
东南大学计算方法上机报告实验报告完整版
实习题11. 用两种不同的顺序计算644834.11000012≈∑=-n n，试分析其误差的变化解：从n=1开始累加，n 逐步增大，直到n=10000；从n=10000开始累加，n 逐步减小，直至1。
算法1的C 语言程序如下： #include<stdio.h> #include<math.h> void main() { float n=0.0; int i; for(i=1;i<=10000;i++) { n=n+1.0/(i*i); } printf("%-100f",n); printf("\n"); float m=0.0; int j; for(j=10000;j>=1;j--) { m=m+1.0/(j*j); } printf("%-7f",m); printf("\n"); }运行后结果如下：结论： 4.设∑=-=Nj N j S 2211，已知其精确值为)11123(21+--N N 。
1）编制按从大到小的顺序计算N S 的程序； 2）编制按从小到大的顺序计算N S 的程序；3）按2种顺序分别计算30000100001000,,S S S ，并指出有效位数。
解：1）从大到小的C语言算法如下：#include<stdio.h>#include<math.h>#include<iostream>using namespace std;void main(){float n=0.0;int i;int N;cout<<"Please input N"<<endl;cin>>N;for(i=N;i>1;i--){n=n+1.0/(i*i-1);N=N-1;}printf("%-100f",n);printf("\n");}执行后结果为：N=2时，运行结果为：N=3时，运行结果为：N=100时，运行结果为：N=4000时，运行结果为：2）从小到大的C语言算法如下：#include<stdio.h>#include<math.h>#include<iostream>using namespace std;void main(){float n=0.0;int i;int N;cout<<"Please input N"<<endl;cin>>N;for(i=2;i<=N;i++){n=n+1.0/(i*i-1);}printf("%-100f",n);printf("\n");}执行后结果为：N=2时，运行结果为：N=3时，运行结果为：N=100时，运行结果为：N=4000时，运行结果为：结论：通过比较可知：N 的值较小时两种算法的运算结果相差不大，但随着N 的逐渐增大，两种算法的运行结果相差越来越大。
应用计算方法
结果分析：比较发现，经过两种改进迭代法，求重根时迭代速度明显加快。
3-4 试验目的体验 Steffensen’s method 加速技巧 试验内容：先用 Newton 法求解方程 x-tanx=0 再用 Steffensen 法求解，比较迭代步数。精确到 10-4。 迭代公式： P(k+1)=P(k)-(P(k)-tan(P(k)))/(1-(sec(P（k）))^2) 初值 P0=1 Newton 法： function [ p,k ]=fnewton( p0,max,tol ) for k=1:max p=p0-(p0-tan(p0))/(1-(sec(p0))^2); if abs(p-p0)<tol break; end p0=p; end disp(p); disp(k) % fnewton( 1,100,10^(-4) ) Steffensen 法： function [ p1,k ]=steffensen( p0,max,tol ) n=1; p(1)=p0; while n<=max for k=1:2 p(k+1)=p(k)-(p(k)-tan(p(k)))/(1-(sec(p(k)))^2); end p1=p(1)-(p(2)-p(1))^2/(p(3)-2*p(2)+p(1)); f0=p1-(p1-tan(p1))/(1-(sec(p1))^2); if abs(f0)<tol break; end n=n+1; p(1)=p1; end disp(p1); disp(n) % steffensen( 1,100,10^(-4) ) >> steffensen( 1,100,10^(-4) ) -1.3387e-07 3
3-2 试验目的:考察 Newton 法求单根的收敛速度 应用 Newton 迭代法求 3-1 中的方程，并与 3-1 中的迭代法相比较，考察收敛速度。精确到时 10-4 迭代公式： P(k+1)= P(k)-(2P(k)-eP(k)+3)/(2- eP(k)) 初值 P0=1 和-1。 Newton 法: function [ p,k ] = fnewton( p0,max,tol ) for k=1:max p=p0-(2*p0-exp(p0)+3)/(2-exp(p0)); if abs(p-p0)<tol break; end
计算方法大作业1 克服Runge现象
x3
x2
x
1
S1 ( x)
-0.34685
0.2086
0.073964
0.038462
S2 (x)
S (xi 0 ) S x(i 0 )
S
'
(xi
0) S
xi' (
0 )i
S
'
'
x(i
0)S
xi' ' (
0)
1 ,n2, . . . , 1
（1）
这里共有了 3n-3 个条件，再加上条件（2）中的 n+1 个插值条件，共有 4n-2 个条件，
因此还需要 2 个方程才能确定 S (x) .通常可在区间[a, b]的端点 a x0,b xn 上各加一个边
dn1
1
2
Mn
dn
（6）
2 1
2
2
2
1 M1 d1
M2
d2
n 1
2
n
1
M
n
1
dn1
n
n 2 M n dn
由式（1）内点拼接条件，可得
i M i1 2M i i M i1 d j i 1, 2,..., n 1
（3） （4）
其中
i
hi 1 hi1
, hi
i
hi hi 1
《线性代数》上机作业一(格式)
课程
《线性代数》
上机
内容
成绩
姓名
雷祥
专业
班级
建筑环境与设备工程
学号
1008020229
教学班
指导
教师
上机
日期
一、上机目的
一、突出教学的பைடு நூலகம்点，减少教学的难点；
二、掌握常用计算方法和处理问题的方法.
二、上机内容
1、向量与矩阵运算。
三、上机作业
1.试分别生成4阶的单位阵和5阶均匀分布的随机矩阵。
2．生成列向量x=[1, 3, 5, 7, 9,…, 15]。
3.试分别生成4阶矩阵A和4阶矩阵B；求A+B，A-B，A*B，A\B，5*B，A^2，A.*B，A.^2及矩阵A的行列式，特征值，特征向量，秩和行最简形.
四、上机心得体会
考试可以带个电脑去就好咯。。。
数值计算方法上机实验报告
数值计算方法上机实验报告实验目的：复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。
利用计算机语言独立编写、调试数值计算方法程序，培养学生利用计算机和所学理论知识分析解决实际问题的能力。
上机练习任务：利用计算机基本C 语言编写并调试一系列数值方法计算通用程序，并能正确计算给定题目，掌握调试技能。
掌握文件使用编程技能，如文件的各类操作，数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。
一、各算法的算法原理及计算机程序框图1. 列主元高斯消去法算法原理：高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘一个方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上对上三角方程组求解。
列选住院是当高斯消元到第k 步时，从k 列的kk a 以下（包括kk a ）的各元素中选出绝对值最大的，然后通过行交换将其交换到kk a 的位置上。
交换系数矩阵中的两行（包括常数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结果。
●源程序：#define N 200#include "stdio.h"#include "math.h"FILE *fp1,*fp2;void LZ(){int n,i,j,k=0,l;double d,t,t1;static double x[N],a[N][N];fp1=fopen("a1.txt","r");fp2=fopen("b1.txt","w");fscanf(fp1,"%d",&n);for(i=0;i<n;++i)for(j=0;j<=n;++j){fscanf(fp1,"%lf",&a[i][j]);}{d=a[k][k];l=k;i=k+1;do{if(fabs(a[i][k])>fabs(d)) /*选主元*/{d=a[i][k];l=i;}i++;}while(i<n);if(d==0){printf("\n输入矩阵有误！\n");}else{ /*换行*/if(l!=k){for(j=k;j<=n;j++){t=a[l][j];a[l][j]=a[k][j];a[k][j]=t;}}}for(j=k+1;j<=n;j++) /*正消*/ a[k][j]/=a[k][k];for(i=k+1;i<n;i++)for(j=k+1;j<=n;j++)a[i][j]-=a[i][k]*a[k][j];k++;}while(k<n);if(k!=0){for(i=n-1;i>=0;i--) /*回代*/ {t1=0;for(j=i+1;j<n;j++)t1+=a[i][j]*x[j];x[i]=a[i][n]-t1;}for(i=0;i<n;i++)fprintf(fp2,"\n 方程组的根为x[%d]=%lf",i+1,x[i]); fclose(fp1); fclose(fp2); }main() { LZ(); }● 具体算例及求解结果：用列选主元法求解下列线性方程组⎪⎩⎪⎨⎧=++=++=-+28x x 23x 2232832321321321x x x x x x 输入3 输出结果：方程组的根为x[1]=6.0000001 2 -3 8 方程组的根为x[2]=4.000000 2 1 3 22 方程组的根为x[3]=2.000000 3 2 1 28● 输入变量、输出变量说明：输入变量：ij a 系数矩阵元素，i b 常向量元素 输出变量：12,,n b b b 解向量元素2. 杜里特尔分解法解线性方程● 算法原理：求解线性方程组Ax b =时，当对A 进行杜里特尔分解，则等价于求解LUx b =，这时可归结为利用递推计算相继求解两个三角形（系数矩阵为三角矩阵）方程组，用顺代，由Ly b =求出y ，再利用回带，由Ux y =求出x 。
数值计算方法上机题
习题二问题：1．编制通用子程序对n+1个节点x i及y i=f(x i) (i=1,…n) （1）n次拉格朗日插值计算公式；（2）n次牛顿向前插值计算公式；（3）n次牛顿向后插值计算公式；（一）程序流程图（1）拉格朗日插值程序流程图（2）牛顿向前插值程序流程图（3）牛顿向后插值程序流程图。
（二）源程序见主程序清单问题：2．计算（1）已知f(x)=lnx,，[a，b]=[1，2]，取h=0.1，x i=1+ih，i=0，1， (10)用通用程序（1），（3）计算ln1.54及ln1.98的近似值；（一）程序清单/* program of question 2.1, page 61 */#include "stdio.h"#include "math.h"main(){ int i,flag=0;double z1,z2,x[11],y[11],t,s1,s2,z[10],c[11][11],log(double),ntb(),L(); for(i=0;i<=10;i++){x[i]=1+0.1*i;y[i]=log(x[i]);}printf("data x:\n");for(i=0;i<=10;i++){flag++;printf("%11.6f",x[i]);if(flag%4==0)printf("\n");}printf("\ndata y:\n");flag=0;for(i=0;i<=10;i++){flag++;printf("%11.6f",y[i]);if(flag%4==0)printf("\n");}printf("\nThe true value:\n");printf(" ln1.54=%f ln1.98=%f\n",log(1.54),log(1.98));z1=L(x,y,10,1.54);z2=L(x,y,10,1.98);t=(1.54-x[10])/0.1;s1=ntb(y,10,t,z,c);s2=ntb(y,10,t,z,c);t=(1.98-x[10])/0.1;printf("The approximate value:\n");printf(" L(1.54)=%f L(1.98)=%f\n",z1,z2);printf(" NTB(1.54)=%f NTB(1.98)=%f\n",s1,s2);}double L(double x[],double y[],int n,double t){ int i,k; double z=0.0,s;if(n==1)z=y[0];for(k=0;k<=n;k++){ s=1.0;for(i=0;i<=n;i++) if(i!=k)s=s*(t-x[i])/(x[k]-x[i]);z=z+s*y[k]; }return z;}double ntb(double y[],int n,double t,double z[],double c[][11]) { int i,j,sn=n;double s;z[0]=t;for(i=1;i<=n-1;i++) z[i]=z[i-1]*(t+i)/(i+1);for(i=0;i<=n;i++) c[i][0]=y[sn--];for(j=1;j<=n;j++)for(i=0;i<=n-j;i++) c[i][j]=c[i][j-1]-c[i+1][j-1];s=y[n];for(i=0;i<=n-1;i++) s=s+z[i]*c[0][i+1];return s;}(二)运行结果data x:1.000000 1.100000 1.200000 1.3000001.400000 1.500000 1.600000 1.7000001.800000 1.9000002.000000data y:0.000000 0.095310 0.182322 0.2623640.336472 0.405465 0.470004 0.5306280.587787 0.641854 0.693147The true value:ln1.54=0.431782 ln1.98=0.683097The approximate value:L(1.54)=0.431782 L(1.98)=0.683097NTB(1.54)=0.431782 NTB(1.98)=0.683097问题：（2）f(x)=1/(1+25x2), |x|≤1取等距节点n=5和n=10，用通用程序（1），（2）依次计算x=-0.95+ih（i=0，1， (19)h=0.1）处f(x)的近似值，并将其结果与其真实值相比较。
数值计算方法上机实习题答案.doc
1．设I n 1 x ndx ，0 5 x（ 1）由递推公式 I n 5I n 11，从 I 0的几个近似值出发，计算I 20；n解：易得： I 0 ln6-ln5=0.1823, 程序为：I=0.182;for n=1:20I=(-5)*I+1/n;endI输出结果为： I 20= -3.0666e+010（ 2）粗糙估计 I 20，用 I n 1 1I n 1 1 ，计算 I 0；5 5n0.0079 1 x 20 1 x 200.0095因为dx I 20dx 6 5所以取 I 20 1(0.0079 0.0095) 0.0087 2程序为： I=0.0087;for n=1:20I=(-1/5)*I+1/(5*n);endII 0= 0.0083（ 3）分析结果的可靠性及产生此现象的原因(重点分析原因 )。
首先分析两种递推式的误差；设第一递推式中开始时的误差为E0 I 0 I 0，递推过程的舍入误差不计。
并记 E n I n I n，则有 E n 5E n 1 ( 5) n E0。
因为 E20( 5) 20 E0 I 20，所此递推式不可靠。
而在第二种递推式中E0 1E1 (1)n E n，误差在缩小，5 5所以此递推式是可靠的。
出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，即算法是否数值稳定。
2．求方程e x10x 2 0 的近似根，要求x k 1x k 5 10 4，并比较计算量。
（1）在 [0， 1]上用二分法；程序： a=0;b=1.0;while abs(b-a)>5*1e-4c=(b+a)/2;if exp(c)+10*c-2>0b=c;else a=c;endendc结果： c =0.0903（ 2）取初值x0 0，并用迭代 x k 1 2 e x ；10程序： x=0;a=1;while abs(x-a)>5*1e-4a=x;x=(2-exp(x))/10;endx结果： x =0.0905（3）加速迭代的结果；程序： x=0;a=0;b=1;while abs(b-a)>5*1e-4a=x;y=exp(x)+10*x-2;z=exp(y)+10*y-2;x=x-(y-x)^2/(z-2*y+x);b=x;endx结果： x =0.0995（ 4）取初值x00 ，并用牛顿迭代法；程序： x=0;a=0;b=1;while abs(b-a)>5*1e-4a=x;x=x-(exp(x)+10*x-2)/(exp(x)+10); b=x;end x 结果： x =0.0905（ 5） 分析绝对误差。
数值分析上机作业1-1解析
数值计算方法上机题目11、实验1. 病态问题实验目的：算法有“优”与“劣”之分，问题也有“好”和“坏”之别。
所谓坏问题就是问题本身的解对数据变化的比较敏感，反之属于好问题。
希望读者通过本实验对此有一个初步的体会。
数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。
病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。
问题提出：考虑一个高次的代数多项式∏=-=---=201)()20)...(2)(1()(k k x x x x x p （E1-1）显然该多项式的全部根为l ，2，…，20，共计20个，且每个根都是单重的（也称为简单的）。
现考虑该多项式方程的一个扰动0)(19=+xx p ε (E1-2)其中ε是一个非常小的数。
这相当于是对（E1-1）中19x 的系数作一个小的扰动。
我们希望比较（E1-1）和（E1-2）根的差别，从而分析方程（E1-1）的解对扰动的敏感性。
实验内容：为了实现方便，我们先介绍两个 Matlab 函数：“roots ”和“poly ”，输入函数u ＝roots （a ）其中若变量a 存储1+n 维的向量，则该函数的输出u 为一个n 维的向量。
设a 的元素依次为121,...,,+n a a a ，则输出u 的各分量是多项式方程0...1121=++++-n n n n a x a x a x a的全部根，而函数b=poly(v)的输出b 是一个n ＋1维变量，它是以n 维变量v 的各分量为根的多项式的系数。
可见“roots ”和“Poly ”是两个互逆的运算函数.ve=zeros(1,21); ve(2)=ess;roots(poly(1:20))+ve)上述简单的Matlab 程序便得到（E1-2）的全部根，程序中的“ess ”即是（E1-2）中的ε。
实验要求：（1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。
计算方法实验上机报告(完整版)
简单迭代法#include<stdio.h>#include<math.h>#define x0 3.0#define MAXREPT 1000#define EPS 1E-6#define G(x) pow(12*x+sin(x)-1,1.0/3)void main(){int i;double x_k=x0,x_k1=x0;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k1);x_k1=G(x_k);if (fabs(x_k1-x_k)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k1,i);return;}x_k=x_k1;}printf("AFTER %d repeate,no solved.\n",MAXREPT);}结果牛顿迭代法一#include<stdio.h>#include<math.h>#define x0 3.0#define MAXREPT 1000#define EPS 1E-6#define G(x) x-(pow(x,3)-sin(x)-12*x+1)/(3*pow(x,2)-cos(x)-12) void main(){int i;double x_k=x0,x_k1=x0;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k1);x_k1=G(x_k);if (fabs(x_k1-x_k)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k1,i);return;}x_k=x_k1;}printf("AFTER %d repeate,no solved.\n",MAXREPT);}结果埃特金加速法#include<stdio.h>#include<math.h>#define x0 3.0#define MAXREPT 1000#define EPS 1E-6#define G(x) (pow(x,3)-sin(x)+1)/12void main(){int i;double x1=x0,x2=x0;double z,y;printf("k\tx1\tx2\txk\n");for(i=0;i<MAXREPT;i++){if(i==0)printf("%d\t\t\t%g\n",i,x2);elseprintf("%d\t%g\t%g\t%g\n",i,y,z,x2);y=G(x1);z=G(y);x2=z-((z-y)*(z-y))/(z-2*y+x1);if (fabs(x2-x1)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x2,i);return;}x1=x2;}printf("AFTER %d repeate,no solved.\n",MAXREPT);} 结果牛顿迭代法二#include<stdio.h>#include<math.h>#define x0 1.5#define MAXREPT 1000#define EPS 1E-6#define G(x) x-(pow(x,3)+pow(x,2)-3*x-3)/(3*pow(x,2)+2*x-3) void main(){int i;double x_k=x0,x_k1=x0;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k1);x_k1=G(x_k);if (fabs(x_k1-x_k)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k1,i);return;}x_k=x_k1;}printf("AFTER %d repeate,no solved.\n",MAXREPT);}结果弦截法#include<stdio.h>#include<math.h>#define x0 0#define x1 1.5#define MAXREPT 1000#define EPS 1E-6#define G(x) pow(x,3)+pow(x,2)-3*x-3void main(){int i;double x_k=x0,x_k1=x1,x_k2=0;double y,z;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k2);y=G(x_k);z=G(x_k1);x_k2=x_k1-(z*(x_k1-x_k))/(z-y);if (fabs(x_k2-x_k1)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k2,i);return;}x_k=x_k1;x_k1=x_k2;}printf("AFTER %d repeate,no solved.\n",MAXREPT); } 结果高斯顺序消元法#include<stdio.h>#include<math.h>#define N 4static double aa[N][N+1]={{2,4,0,1,1},{3,8,2,2,3},{1,3,3,0,6},{2,5,2,2,3}}; int gauss(double a[][N+2],double x[]);void putout(double a[][N+2]);void main(){int i,j,det;double a[N+1][N+2],x[N+1];for(i=1;i<=N;i++)for(j=1;j<=N+1;j++)a[i][j]=aa[i-1][j-1];det=gauss(a,x);if(det!=0)for(i=1;i<=N;i++)printf(" x[%d]=%g",i,x[i]);printf("\n");}int gauss(double a[][N+2],double x[]){int i,j,k;double c;putout(a);for(k=1;k<=N-1;k++){ if(fabs(a[k][k])<1e-17){printf("\n pivot element is 0.fail!\n");return 0;}for(i=k+1;i<=N;i++){c=a[i][k]/a[k][k];for(j=k;j<=N+1;j++){a[i][j]=a[i][j]-c*a[k][j];}}putout(a);}if(fabs(a[N][N])<1e-17){printf("\n pivot element is 0.fail!\n");return 0;}for(k=N;k>=1;k--){x[k]=a[k][N+1];for(j=k+1;j<=N;j++){x[k]=x[k]-a[k][j]*x[j];}x[k]=x[k]/a[k][k];}return(1);}void putout(double a[][N+2]){for(int i=1;i<=N;i++){for(int j=1;j<=N+1;j++)printf("%-15g",a[i][j]);printf("\n");}printf("\n");}结果雅克比迭代法#include<stdio.h>#include<math.h>#define N 5#define EPS 0.5e-4static double aa[N][N]={{4,-1,0,-1,0},{-1,4,-1,0,-1},{0,-1,4,-1,0},{-1,0,-1,4,-1},{0,-1,0,-1,4}}; static double bb[N]={2,1,2,1,2};void main(){int i,j,k,NO;double a[N+1][N+1],b[N+1],x[N+1],y[N+1];double d,sum,max;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];}printf("\n 请输入最大迭代次数（尽量取大值！）NO:");scanf("%d",&NO);printf("\n");for(i=1;i<=N;i++)x[i]=0;k=0;printf(" k",' ');for(i=1;i<=N;i++)printf("%8cx[%d]",' ',i);printf("\n 0");for(i=1;i<=N;i++)printf("%12.8g",x[i]);printf("\n");do{for(i=1;i<=N;i++){sum=0.0;for(j=1;j<=N;j++)if(j!=i) sum=sum+a[i][j]*x[j];y[i]=(-sum+b[i])/a[i][i];}max=0.0;for(i=0;i<=N;i++){d=fabs(y[i]-x[i]);if(max<d) max=d;x[i]=y[i];}printf("%6d",k+1);for(i=1;i<=N;i++)printf("%12.8g",x[i]);printf("\n");k++;}while((max>=EPS)&&(k<NO));printf("\nk=%d\n",k);if(k>=NO) printf("\nfail!\n");elsefor(i=1;i<=N;i++)printf("x[%d]=%g\t",i,x[i]);}结果拉格朗日插值多项式#include<stdio.h>#include<math.h>#define N 4doublex[N]={0.56160,0.56280,0.56401,0.56521},y[N]={0.82741,0.82659,0.82577,0.82495}; void main(){double x=0.5635;double L(double xx);double lagBasis(int k,double xx);void output();output();printf("\n近似值L(%g)=%g\n",x,L(x));}double lagBasis(int k,double xx){double lb=1;int i;for(i=0;i<N;i++)if(i!=k) lb*=(xx-x[i])/(x[k]-x[i]);return lb;}double L(double xx){double s=0;int i;for(i=0;i<=N;i++)s+=lagBasis(i,xx)*y[i];return s;}void output(){int i;printf("\n各节点信息:\nxi:");for(i=0;i<N;i++)printf("\t%g",x[i]);printf("\nyi:");for(i=0;i<N;i++)printf("\t%g",y[i]);}结果牛顿插值多项式#include <math.h>#include <stdio.h>int a;#define M 4double x[M+1]={0.4,0.55,0.65,0.8,0.9},y[M+1]={0.41075,0.57815,0.69675,0.88811,1.02652}; void main(){double x;printf("输入x=");scanf("%lf",&x);printf("次数：");scanf("%d",&a);double N(double xx,int a);void output();output();printf("\n%d次牛顿插值多项式N(%g)=%g\n",a,x,N(x,a));}double N(double xx,int a){double s=y[0],d=1;int i,j;double df[M+1][M+1];for(i=0;i<=M;i++)df[i][0]=y[i];for(j=1;j<=a;j++)for(i=j;i<=a;i++)df[i][j]=(df[i][j-1]-df[i-1][j-1])/(x[i]-x[i-j]);printf("\nx\tf(x)\t");for(j=1;j<=a;j++) printf("%5d阶",j);for(i=0;i<=a;i++){printf("\n%g\t%g",x[i],y[i]);for(j=1;j<=i;j++)printf("\t%7.5g",df[i][j]);}for(i=1;i<=a;i++){d*=(xx-x[i-1]);s+=df[i][i]*d;}return s;}void output(){int i;printf("\n各节点信息:\nxi:");for(i=0;i<=M;i++)printf("\t%7g",x[i]);printf("\nyi:");for(i=0;i<=M;i++)printf("\t%7g",y[i]);}结果复合梯形公式#include<stdio.h>#include<math.h>#define f(x) 1/(x*x+1)#define Pi 3.1415926void main(){double a=0,b=1;double T,h,x;int n,i;printf("please input n:");scanf("%d",&n);h=(b-a)/n;x=a;T=0;for(i=1;i<n;i++){x+=h;T+=f(x);}T=(f(a)+2*T+f(b))*h/2;printf("T(%d)=%g\n",n,T);printf("The exact value is %g\n",Pi/4);}复合辛普森公式#include<stdio.h>#include<math.h>#define f(x) 1/(1+x*x)#define Pi 3.1415926void main(){double a=0,b=1;double S,h,x;int n,i;printf("please input Even n:");scanf("%d",&n);h=(b-a)/n;x=a; S=0;for(i=1;i<n;i++){x+=h;if(i%2==0) S+=2*f(x);else S+=4*f(x);}S=(f(a)+S+f(b))*h/3;printf("S(%d)=%g\n",n,S);printf("The exact value is %g\n",Pi/4);}龙贝格公式加速#include<stdio.h>#include<math.h>#define f(x) sin(x)/(1+x)#define M 3void main(){double a=0,b=1;double Long(double a,double b);printf("近似值I=%g\n",Long(a,b));}double Long(double a,double b){int n=1,i=1,j=1;double T[M+1][M+1],h,x,sum;h=b-a;T[0][0]=(f(a)+f(b))/2;for(j=1;j<=3;j++){x=a;h/=2;n*=2;sum=0;for(i=1;i<=n;i+=2){x=a+i*h;sum+=f(x);}T[j][0]=T[j-1][0]/2+h*sum;}for(i=1;i<=M;i++)for(j=1;j<=i;j++){T[i][j]=(pow(4,j)*T[i][j-1]-T[i-1][j-1])/(pow(4,j)-1);}printf("k\tT0\tT1\tT2\tT3\n");for(i=0;i<=M;i++){printf("%d",i);for(j=0;j<=i;j++)printf(" %g",T[i][j]);printf("\n");}return T[M][M];}。
计算方法作业
应用计算方法学院：自动化学院班级：硕1104班姓名：陈兆辉计算方法2班教师：张晓丹目录上机作业1 (3)上机作业2 (6)上机作业3 (10)上机作业4 (12)上机作业5 (17)上机作业6 (18)上机作业7 (19)上机作业11.分别用不动带你迭代法与Newton 法求解方程2x- ex +3=0的正根和负根解：（1）用不动点迭代法求方程的正根1）迭代公式：()231ln k P k P ++=，初值p0=1.02）程序：k=1;Tol=0.00001;p0=1.0;N=100;while k<=Np=log(2*p0+3);if abs(p-p0)<Tolbreak;endk=k+1;p0=p;enddisp(p);disp(k)3)显示结果：P= 1.9239K=11（2） 用不动点迭代法求解方程的负根1）迭代公式： 1(3)/2k p k p e +=-，初值p0=-1.02）程序：k=1;Tol=0.00001;p0=-1.0;N=100;while k<=Np=(exp(p0)-3)/2;if abs(p-p0)<Tolbreak;endk=k+1;p0=p;enddisp(p);disp(k)3）显示结果：P= -1.3734k=7(3)用牛顿法求方程的正根1）迭代公式：()1232k k p k k k p p e p p e +-+=--，初值p0=1.02）程序：k=1;Tol=0.00001;p0=1.0;N=100;while k<=Np=p0-(2*p0-exp(p0)+3)/(2-exp(p0));if abs(p-p0)<Tolbreak;endk=k+1;p0=p;enddisp(p);disp(k)3）显示结果：P=1.9239K=8（4）用牛顿法求解方程的负根1）迭代公式：()1232k k p k k k p p e p p e +-+=--，初值p0=-1.02）程序：k=1;Tol=0.00001;p0=-1.0;N=100;while k<=Np=p0-(2*p0-exp(p0)+3)/(2-exp(p0));if abs(p-p0)<Tolbreak;endk=k+1;p0=p;enddisp(p);disp(k)3）显示结果：P= -1.3734K= 44）结果分析：不动点迭代法是求解非线性方程的主要方法，其中牛顿法应用最广泛，它求解方程的单根时具有二阶收敛性，但是它要求较好的初始近似值，牛顿法比普通的迭代法收敛速度快，能较少的迭代达到理想的值。
计算方法上机作业
计算⽅法上机作业请注意! Pay attention!计算⽅法作业要求1.上机内容：只要求完成《计算⽅法》实验报告那⼀部分的内容。
2.试验报告要求：（1）要求写出解题的程序代码和结果，以及必要的⽂字说明。
（2）报告要求打印，封⾯上写清班级、学号和姓名。
（3）提交报告时间：2013年6⽉13号（16周周四）。
⽬录第⼀部分基本操作 (1)§1.Matlab的使⽤ (1)1-1.直接输⼊命令 (1)1-2.⽤M⽂件开发程序 (1)§2.M⽂件程序的主要语句和主要函数 (2)2-1.Matlab的数字特征 (2)2-2.主要语句 (3)2-3.常⽤函数 (4)2-4.⼏个常⽤的命令 (5)§3.矩阵的有关计算 (5)3-1.矩阵的输⼊ (5)3-2.矩阵/向量的运算 (6)3-3.矩阵的范数 (6)3-4.向量的范数 (7)3-5.矩阵的条件数 (7)3-6.矩阵的特征值和特征向量 (8)§4.Matlab绘图 (9)4-1.绘图的基本命令 (9)4-2.图形的交互编辑 (11)第⼆部分数值计算 (12)§1.⽅程求根 (12)1-1.⽜顿迭代法 (12)1-2.图解法确定迭代的初始点 (13)§2.线性⽅程组 (13)2-1.迭代法的收敛性 (13)2-2.线性⽅程组的病态问题 (14)2-3.求解线性⽅程组 (15)§3.插值和拟合 (16)/doc/b4c2d2dfc850ad02de8041fc.html grange插值 (16)3-2.代数多项式插值 (17)3-3.插值误差 (17)3-4.分段线性插值 (18)3-5.数据的曲线拟合 (18)§4.数值积分 (20)4-1.复合梯形求积公式 (20)4-2.复合Simpson求积公式 (20)§5.常微分⽅程的数值解法 (21)5-1.Euler⽅法 (21)5-2.改进的Euler⽅法 (23)5-3.四阶龙格-库塔⽅法 (24)习题 (27)⼀、⽅程求根 (27)⼆、线性⽅程组 (27)三、插值与拟合 (28)四、数值积分 (29)五、常微分⽅程 (30)《计算⽅法》实验报告 (31)⼀、⽅程求根 (31)⼆、线性⽅程组 (31)三、插值与拟合 (32)四、数值积分 (32)五、常微分⽅程 (33)第⼀部分基本操作§1.1.MatlabMatlab 的使⽤Matlab 的使⽤⽅法有两种：（1）在Matlab 的命令窗⼝（Matlab Command Windows ）中直接输⼊命令，即可得到结果；（2）在Matlab 的编辑窗⼝（Matlab Editor ）内编写M ⽂件，然后在命令窗⼝执⾏该⽂件，得到所需的结果。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
计算方法第四次上机报告
2.用欧拉方法解初值
y’=10x(1-y) 0<=x<=1
Y(0)=0
取步长h=0.1,保留5位有效数字，并与准确解相比较
分析：该题目考察欧拉方法解初值问题
程序如下：
function Heun(a,b,y0,n)
h=(b-a)/n;x=a:h:b;
y=y0*ones(1,n+1);
for j=2:n+1
yp=y(j-1)+h*f(x(j-1),y(j-1));
yc=y(j-1)+h*f(x(j),yp);
y(j)=1/2*(yp+yc);
end
for k=1:n+1
fprintf('x[%d]=%f\ty[%d]=%f\n',k-1,x(k),k-1,y(k)); end
function z=f(xx,yy)
z=10*xx*(1-yy);
运行结果：
>> Heun(0,1,0,10)
x[0]=0.000000 y[0]=0.000000
x[1]=0.100000 y[1]=0.050000
x[2]=0.200000 y[2]=0.183000
x[3]=0.300000 y[3]=0.362740
x[4]=0.400000 y[4]=0.547545
x[5]=0.500000 y[5]=0.705905
x[6]=0.600000 y[6]=0.823543
x[7]=0.700000 y[7]=0.901184
x[8]=0.800000 y[8]=0.947627
x[9]=0.900000 y[9]=0.973290
x[10]=1.000000 y[10]=0.986645
>>
分析：
该结果与准确结果相比比较接近，但是有一定的误差。
6.用四阶龙格—库塔公式解第三题中的初值问题，取步长h=0.2，保留五位有效数字。
题目目的分析：
该题考查四阶龙格-库塔方法和改进欧拉方法求解精确度问题。
程序：
改进欧拉法：
function Heun(a,b,y0,n)
h=(b-a)/n;x=a:h:b;
y=y0*ones(1,n+1);
for j=2:n+1
yp=y(j-1)+h*f(x(j-1),y(j-1));
yc=y(j-1)+h*f(x(j),yp);
y(j)=1/2*(yp+yc);
end
for k=1:n+1
fprintf('x[%d]=%f\ty[%d]=%f\n',k-1,x(k),k-1,y(k));
end
四阶龙格-库塔方法：
function odeRK4(a,b,y0,n)
h=(b-a)/n;x=a:h:b;
y=y0*ones(1,n+1);
for j=2:n+1
k1=f(x(j-1),y(j-1));
k2=f(x(j-1)+h/2,y(j-1)+h/2*k1);
k3=f(x(j-1)+h/2,y(j-1)+h/2*k2);
k4=f(x(j-1)+h,y(j-1)+h*k3);
y(j)=y(j-1)+h/6*(k1+k4)+h/3*(k2+k3);
end
for k=1:n+1
fprintf('x[%d]=%f\ty[%d]=%f\n',k-1,x(k),k-1,y(k)); end
function z=f(xx,yy)
z=-yy;
结果：
>> Heun(0,1.0,1,10)
x[0]=0.000000 y[0]=1.000000
x[1]=0.100000 y[1]=0.905000
x[2]=0.200000 y[2]=0.819025
x[3]=0.300000 y[3]=0.741218
x[4]=0.400000 y[4]=0.670802
x[5]=0.500000 y[5]=0.607076
x[6]=0.600000 y[6]=0.549404
x[7]=0.700000 y[7]=0.497210
x[8]=0.800000 y[8]=0.449975
x[9]=0.900000 y[9]=0.407228
x[10]=1.000000 y[10]=0.368541
>>
>> odeRK4(0,1.0,1,5)
x[0]=0.000000 y[0]=1.000000
x[1]=0.200000 y[1]=0.818733
x[2]=0.400000 y[2]=0.670324
x[3]=0.600000 y[3]=0.548817
x[4]=0.800000 y[4]=0.449335
x[5]=1.000000 y[5]=0.367885
>>
结果分析：
由程序结果得改进欧拉法和四阶龙格-库塔法运行结果上在小数第二位后才开始有所不同，但结果相差不大，精确度都是挺高的
12.将二阶方程
化为一阶方程组，取h=0.1,用四阶龙格-库塔法求有y（0.2）的近似值，保留5位有效数字。
解：
目的分析：
该题目是要将二阶方程化为一解方程组，并用四阶龙格-库塔法求近似值。
程序如下：
公式部分：
function z=f1(xx,yy)
z=5*xx-6*yy;
四阶龙格-库塔法：
function odeRK4_1(a,b,y0,n)
h=(b-a)/n;x=a:h:b;
y=y0*ones(1,n+1);
z=-ones(1,n+1);
for j=2:n+1
k1=z(j-1);
l1=f1(z(j-1),y(j-1));
k2=z(j-1)+h/2*l1;
l2=f1(z(j-1)+h/2*l1,y(j-1)+h/2*k1);
k3=z(j-1)+h/2*l2;
l3=f1(z(j-1)+h*l2/2,y(j-1)+h/2*k2);
k4=z(j-1)+h*l3;
l4=f1(z(j-1)+h*l3,y(j-1)+h*k3);
y(j)=y(j-1)+h*(k1+2*k2+2*k3+k4)/6;
z(j)=z(j-1)+h*(l1+2*l2+2*l3+l4)/6;
end
for k=1:n+1
fprintf('x[%d]=%f\ty[%d]=%f\n',k-1,x(k),k-1,y(k));
end
end
结果图：
结果分析：
由程序结果得四阶龙格-库塔法所得y(0.2)=0.501088
所给初值问题精确解为,y(0.2)=0.5009424.两结果在数值上相差不大，精确度高。