计算方法第七章上机报告
计算方法上机实验报告
《计算方法》上机实验报告班级：XＸＸXXX小组成员:XXＸXXＸＸXXXＸXXXXXＸＸXXXＸXXXXXＸ任课教师：XXX二〇一八年五月二十五日前言通过进行多次得上机实验,我们结合课本上得内容以及老师对我们得指导,能够较为熟练地掌握Ｎewton迭代法、Ｊacｏｂi迭代法、Gausｓ－Seｉｄｅl 迭代法、Ｎewｔon 插值法、Ｌａgｒanｇe 插值法与Ｇaｕｓs 求积公式等六种算法得原理与使用方法，并参考课本例题进行了MATＬAB程序得编写。
以下为本次上机实验报告,按照实验内容共分为六部分.实验一：一、实验名称及题目:Newｔon迭代法例2、7(P3８）:应用Newｔon迭代法求在附近得数值解，并使其满足、二、解题思路：设就是得根,选取作为初始近似值，过点做曲线得切线,得方程为,求出与轴交点得横坐标,称为得一次近似值,过点做曲线得切线,求该切线与轴得横坐标称为得二次近似值，重复以上过程，得得近似值序列，把称为得次近似值，这种求解方法就就是牛顿迭代法。
三、Matlab程序代码：fｕncｔion newtｏn_itｅrａｔioｎ(x０,tol）syms z %定义自变量fｏrmaｔlｏnｇ％定义精度ｆ=ｚ*z＊z-z－1;f1=diff(f）；％求导y=suｂs(ｆ，ｚ,x0);y1=sｕbs(f１,ｚ，x0)；%向函数中代值x1=x0－y/y1; k=1；wｈilｅａｂｓ(ｘ1—x0)〉=ｔolx０=x1;y=subｓ（f,z,ｘ0)；y1=subｓ(ｆ1,ｚ,ｘ0);x1＝x0-ｙ/y1；k=k+1;ｅndｘ=douｂle（x１)Ｋ四、运行结果：实验二:一、实验名称及题目：Jacｏbｉ迭代法例３、7(P74）:试利用Jacobi迭代公式求解方程组要求数值解为方程组得精确解、二、解题思路：首先将方程组中得系数矩阵分解成三部分,即:，为对角阵,为下三角矩阵,为上三角矩阵。
之后确定迭代格式,,（, 即迭代次数),称为迭代矩阵。
计算方法上机实验报告-MATLAB
《计算方法》实验报告指导教师：学院：班级：团队成员：一、题目例2.7应用Newton 迭代法求方程210x x --=在1x =附近的数值解k x ，并使其满足8110k k x x ---<原理：在方程()0f x =解的隔离区间[],a b 上选取合适的迭代初值0x ，过曲线()y f x =的点()()00x f x ，引切线()()()1000:'l y f x f x x x =+-其与x 轴相交于点：()()0100 'f x x x f x =-，进一步,过曲线()y f x =的点()()11x f x , 引切线()()()2111: 'l y f x f x x x =+-其与x 轴相交于点：()()1211 'f x x x f x =-如此循环往复，可得一列逼近方程()0f x =精确解*x 的点01k x x x ,,,,,其一般表达式为：()()111 'k k k k f x x x f x ---=-该公式所表述的求解方法称为Newton 迭代法或切线法。
程序：function y=f(x)%定义原函数y=x^3-x-1;endfunction y1=f1(x0)%求导函数在x0点的值syms x;t=diff(f(x),x);y1=subs(t,x,x0);endfunction newton_iteration(x0,tol)%输入初始迭代点x0及精度tol x1=x0-f(x0)/f1(x0);k=1;%调用f函数和f1函数while abs(x1-x0)>=tolx0=x1;x1=x0-f(x0)/f1(x0);k=k+1;endfprintf('满足精度要求的数值为x(%d)=%1.16g\n',k,x1);fprintf('迭代次数为k=%d\n',k);end结果：二、题目例3.7试利用Jacobi 迭代公式求解方程组1234451111101112115181111034x x x x ----⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪= ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪---⎝⎭⎝⎭⎝⎭ 要求数值解()k X 满足()4210k --≤X X ，其中T (1,2,3,4)=X 为方程组的精确解。
(完整word版)计算方法上机报告-备份
计算方法上机作业报告姓名：李小盼班级：计算方法B3班学号：6230489477专业：热能工程2016年《计算方法B》上机题目一.计算机语言要求使用语言不做限制，可以使用C、C++、FORTRAN、VC、VB、C#、Matlab、PHP、JavaScript等语言。
二.实习报告内容上机题目完成后，必须交一份上机报告。
上机报告中应对每道题目包括：（1）上机题目内容;（2）详细说明实现的思想、算法依据、算法实现的结构;（3）详细完整的源程序，并附相关的注释说明;（4）给出必要的计算结果（若数据量大，则只需列出主要的数据内容），并对结果进行分析;（5）对上机中出现的问题进行分析总结；三.实习报告要求1.提供一份完整的上机报告的电子文档；然后再提供一份与电子文档对应的纸质上机报告。
2.电子文档中提供上机过程中的所有源程序、输出数据，以及可以运行的文件。
如果程序的运行环境特殊，请附上运行程序所需要的软件环境。
3.上机报告严禁抄袭，如发现有抄袭现象，所有涉及抄袭的上机报告将被以作弊处理，并按零分处理,不再另行通知。
4.上机报告电子版一、二、三班发送到邮箱:xjtujsff@,上机作业纸面作业请送到:理科楼338。
上机实习题目1.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。
已探测到一组等分点位置的深度数据(单位:米)如下表所示：(1)请用合适的曲线拟合所测数据点；(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；解：1.1 实现思想选用曲线拟合数据点时，一方面要满足插值条件，即保证所得曲线经过以上所有数据点，另一方面也要保证曲线具有足够的“光滑性”，故这里采用三次样条插值法构造插值函数。
为构成封闭方程组，边界条件选取自然三次样条，以获得三次样条插值的关键参数M0和M n。
1.2 算法的依据以及结构三次样条插值具有较好的光滑性，可以用于对数据点的光滑拟合。
计算方法上机实验
1.拉格朗日插值多项式，用于离散数据的拟合#include <stdio.h>#include <conio.h>#include <alloc.h>float lagrange(float *x,float *y,float xx,int n) /*拉格朗日插值算法*/{ int i,j;float *a,yy=0.0; /*a作为临时变量，记录拉格朗日插值多项式*/a=(float *)malloc(n*sizeof(float));for(i=0;i<=n-1;i++){ a[i]=y[i];for(j=0;j<=n-1;j++)if(j!=i) a[i]*=(xx-x[j])/(x[i]-x[j]);yy+=a[i];}free(a);return yy;}main(){ int i,n;float x[20],y[20],xx,yy;printf("Input n:");scanf("%d",&n);if(n>=20) {printf("Error!The value of n must in (0,20)."); getch();return 1;} if(n<=0) {printf("Error! The value of n must in (0,20)."); getch(); return 1;} for(i=0;i<=n-1;i++){ printf("x[%d]:",i);scanf("%f",&x[i]);}printf("\n");for(i=0;i<=n-1;i++){ printf("y[%d]:",i);scanf("%f",&y[i]);}printf("\n");printf("Input xx:");scanf("%f",&xx);yy=lagrange(x,y,xx,n);printf("x=%f,y=%f\n",xx,yy);getch();}2.牛顿插值多项式，用于离散数据的拟合#include <stdio.h>#include <conio.h>#include <alloc.h>void difference(float *x,float *y,int n){ float *f;int k,i;f=(float *)malloc(n*sizeof(float));for(k=1;k<=n;k++){ f[0]=y[k];for(i=0;i<k;i++)f[i+1]=(f[i]-y[i])/(x[k]-x[i]);y[k]=f[k];}return;}main(){ int i,n;float x[20],y[20],xx,yy;printf("Input n:");scanf("%d",&n);if(n>=20) {printf("Error! The value of n must in (0,20)."); getch(); return 1;} if(n<=0) {printf("Error! The value of n must in (0,20).");getch(); return 1;} for(i=0;i<=n-1;i++){ printf("x[%d]:",i);scanf("%f",&x[i]);}printf("\n");for(i=0;i<=n-1;i++){ printf("y[%d]:",i);scanf("%f",&y[i]);}printf("\n");difference(x,(float *)y,n);printf("Input xx:");scanf("%f",&xx);yy=y[20];for(i=n-1;i>=0;i--) yy=yy*(xx-x[i])+y[i];printf("NewtonInter(%f)=%f",xx,yy);getch();}3.高斯列主元消去法,求解其次线性方程组第一种#include<stdio.h>#include <math.h>#define N 20int main(){ int n,i,j,k;int mi,tmp,mx;float a[N][N],b[N],x[N];printf("\nInput n:");scanf("%d",&n);if(n>N){ printf("The input n should in(0,N)!\n");getch();return 1;}if(n<=0){ printf("The input n should in(0,N)!\n");getch();return 1;}printf("Now input a(i,j),i,j=0...%d:\n",n-1); for(i=0;i<n;i++){ for(j=0;j<n;j++)scanf("%f",&a[i][j]);}printf("Now input b(i),i,j=0...%d:\n",n-1);for(i=0;i<n;i++)scanf("%f",&b[i]);for(i=0;i<n-2;i++){ for(j=i+1,mi=i,mx=fabs(a[i][j]);j<n-1;j++) if(fabs(a[j][i])>mx){ mi=j;mx=fabs(a[j][i]);}if(i<mi){ tmp=b[i];b[i]=b[mi];b[mi]=tmp;for(j=i;j<n;j++){ tmp=a[i][j];a[i][j]=a[mi][j];a[mi][j]=tmp;}}for(j=i+1;j<n;j++){ tmp=-a[j][i]/a[i][i];b[j]+=b[i]*tmp;for(k=i;k<n;k++)a[j][k]+=a[i][k]*tmp;}}x[n-1]=b[n-1]/a[n-1][n-1];for(i=n-2;i>=0;i--){ x[i]=b[i];for(j=i+1;j<n;j++)x[i]-=a[i][j]*x[j];x[i]/=a[i][i];}for(i=0;i<n;i++)printf("Answer:\n x[%d]=%f\n",i,x[i]); getch();return 0;}第二种#include<math.h>#include<stdio.h>#define NUMBER 20#define Esc 0x1b#define Enter 0x0dfloat A[NUMBER][NUMBER+1] ,ark;int flag,n;exchange(int r,int k);float max(int k);message();main(){float x[NUMBER];int r,k,i,j;char celect;clrscr();printf("\n\nUse Gauss.");printf("\n\n1.Jie please press Enter."); printf("\n\n2.Exit press Esc.");celect=getch();if(celect==Esc)exit(0);printf("\n\n input n=");scanf("%d",&n);printf(" \n\nInput matrix A and B:"); for(i=1;i<=n;i++){printf("\n\nInput a%d1--a%d%d and b%d:",i,i,n,i);for(j=1;j<=n+1;j++) scanf("%f",&A[i][j]);}for(k=1;k<=n-1;k++){ark=max(k);if(ark==0){printf("\n\nIt's wrong!");message();}else if(flag!=k)exchange(flag,k);for(i=k+1;i<=n;i++)for(j=k+1;j<=n+1;j++)A[i][j]=A[i][j]-A[k][j]*A[i][k]/A[k][k];}x[n]=A[n][n+1]/A[n][n];for( k=n-1;k>=1;k--){float me=0;for(j=k+1;j<=n;j++){me=me+A[k][j]*x[j];}x[k]=(A[k][n+1]-me)/A[k][k];}for(i=1;i<=n;i++){printf(" \n\nx%d=%f",i,x[i]);}message();}exchange(int r,int k){int i;for(i=1;i<=n+1;i++)A[0][i]=A[r][i];for(i=1;i<=n+1;i++)A[r][i]=A[k][i];for(i=1;i<=n+1;i++)A[k][i]=A[0][i];}float max(int k){int i;float temp=0;for(i=k;i<=n;i++)if(fabs(A[i][k])>temp){temp=fabs(A[i][k]);flag=i;}return temp;}message(){printf("\n\n Go on Enter ,Exit press Esc!");switch(getch()){case Enter: main();case Esc: exit(0);default:{printf("\n\nInput error!");message();} }}4.龙贝格求积公式，求解定积分#include<stdio.h>#include<math.h>#define f(x) (sin(x)/x)#define N 20#define MAX 20#define a 2#define b 4#define e 0.00001float LBG(float p,float q,int n){ int i;float sum=0,h=(q-p)/n;for (i=1;i<n;i++)sum+=f(p+i*h);sum+=(f(p)+f(q))/2;return(h*sum);}void main(){ int i;int n=N,m=0;float T[MAX+1][2];T[0][1]=LBG(a,b,n);n*=2;for(m=1;m<MAX;m++){ for(i=0;i<m;i++)T[i][0]=T[i][1];T[0][1]=LBG(a,b,n);n*=2;for(i=1;i<=m;i++)T[i][1]=T[i-1][1]+(T[i-1][1]-T[i-1][0])/(pow(2,2*m)-1);if((T[m-1][1]<T[m][1]+e)&&(T[m-1][1]>T[m][1]-e)){ printf("Answer=%f\n",T[m][1]); getch();return ;}}}5.牛顿迭代公式，求方程的近似解#include<stdio.h>#include<math.h>#include<conio.h>#define N 100#define PS 1e-5#define TA 1e-5float Newton(float (*f)(float),float(*f1)(float),float x0 ) { float x1,d=0;int k=0;do{ x1= x0-f(x0)/f1(x0);if((k++>N)||(fabs(f1(x1))<PS)){ printf("\nFailed!");getch();exit();}d=(fabs(x1)<1?x1-x0:(x1-x0)/x1);x0=x1;printf("x(%d)=%f\n",k,x0);}while((fabs(d))>PS&&fabs(f(x1))>TA) ;return x1;}float f(float x){ return x*x*x+x*x-3*x-3; }float f1(float x){ return 3.0*x*x+2*x-3; }void main(){ float f(float);float f1(float);float x0,y0;printf("Input x0: ");scanf("%f",&x0);printf("x(0)=%f\n",x0);y0=Newton(f,f1,x0);printf("\nThe root is x=%f\n",y0); getch();}。
数值计算方法上机实验报告
数值计算方法上机实验报告
一、实验目的
本次实验的主要目的是熟悉和掌握数值计算方法，学习梯度下降法的
原理和实际应用，熟悉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库的使用方法。
计算方法上上机实习报告
计算方法上上机实习报告在本次计算方法的上机实习中，我深入体验了数值计算的魅力和挑战，通过实际操作和实践，对计算方法有了更深刻的理解和认识。
实习的目的在于将课堂上学到的理论知识运用到实际的计算中，熟悉各种数值算法的实现过程，提高编程能力和解决实际问题的能力。
我们使用了具体编程语言和软件名称进行编程和计算。
在实习过程中，我首先接触到的是数值逼近的相关内容。
通过多项式插值和曲线拟合的练习，我明白了如何用简单的函数去近似复杂的曲线。
例如，拉格朗日插值法和牛顿插值法让我能够根据给定的离散数据点构建出一个连续的函数，从而对未知点进行预测。
在实际操作中，我需要仔细处理数据的输入和输出，以及算法中的细节，如边界条件和误差控制。
数值积分是另一个重要的部分。
通过梯形公式和辛普森公式，我学会了如何对给定的函数进行数值积分。
在编程实现时，要合理地选择积分区间和步长，以达到所需的精度。
同时，我也了解到了数值积分方法的误差来源和误差估计方法，这对于评估计算结果的可靠性非常重要。
线性方程组的求解是计算方法中的核心内容之一。
我分别使用了高斯消元法和迭代法（如雅克比迭代法和高斯赛德尔迭代法）来求解线性方程组。
在实际编程中，我深刻体会到了算法的效率和稳定性的重要性。
对于大规模的线性方程组，选择合适的算法可以大大提高计算速度和精度。
在非线性方程求根方面，我运用了二分法、牛顿法和割线法等方法。
这些方法各有特点，二分法简单但收敛速度较慢，牛顿法收敛速度快但需要计算导数。
在实际应用中，需要根据方程的特点和求解的要求选择合适的方法。
在实习中，我也遇到了不少问题和挑战。
首先是编程中的错误，如语法错误、逻辑错误等，这需要我耐心地调试和修改代码。
其次，对于一些复杂的算法，理解其原理和实现细节并不容易，需要反复查阅资料和思考。
还有就是数值计算中的误差问题，有时候由于误差的积累，导致计算结果与预期相差较大，需要通过调整算法参数或者采用更精确的算法来解决。
现代数值计算上机实验报告资料
太原科技大学现代数值计算方法上机报告院系：华科学院专业年级：计算机科学与技术学生姓名：张栩嘉学号：201522030129指导教师：2017年5月12日数值计算方法上机实习题1． 设⎰+=105dx xx I nn ， （1） 由递推公式nI I n n 151+-=-，从0I 的几个近似值出发，计算20I ； （2） 粗糙估计20I ，用nI I n n 51511+-=-，计算0I ； （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。
%上机习题1 %（1） I=0.182; for n=1:20I=(-5)*I+1/n; endfprintf('I20 的值是 %e\n',I); %(2)I=0.0087; for n=1:20I=(-1/5)*I+1/(5*n); endfprintf('I0 的值是 %f\n',I);3） 现象产生的原因：假设S n 的真值为S n ∗，误差为εn ，即εn =S n ∗−S n 对于真值，我们也有关系式S n ∗+5S n−1∗=1n 综合两个递推公式有εn =−5εn−1 这就意味着哪怕开始只有一点点误差，只要n 足够大，按照这种每计算一步误差增长5倍的方式，所得结果总是不可信的。
因此整个算法是不稳定的。
对于第二种方法 εn =−15εn−1 误差会以每计算一步缩小到1/5的方式进行，所以以这种方法计算的结果是可靠的，整个算法是稳定的。
2. 求方程0210=-+x e x的近似根，要求41105-+⨯<-k k x x ，并比较计算量。
（1） 在[0，1]上用二分法； （2） 取初值00=x ，并用迭代1021x k e x -=+；（3） 加速迭代的结果；（4） 取初值00=x ，并用牛顿迭代法； （5） 分析绝对误差。
%（1）在[0，1]上用二分法；a=0;b=1.0;ci=0;while abs(b-a)>5*1e-4c=(b+a)/2;ci=ci+1;if exp(c)+10*c-2>0b=c;else a=c;endendfprintf('二分法求得 c 的值是 %f\t',c);fprintf('迭代次数是 %d\n',ci);%（2）不动点迭代法，x=0;a=1;ci=0;while abs(x-a)>5*1e-4a=x;x=(2-exp(x))/10;ci=ci+1;endfprintf('不动点迭代求得 x 的值是 %f\t',x); fprintf('迭代次数是 %d\n',ci);%（3）艾特肯加速迭代x=0;a=0;b=1;ci=0;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;ci=ci+1;endfprintf('艾特肯加速迭代求得 x 的值是 %f\t',x); fprintf('迭代次数是 %d\n',ci);%（4）用牛顿迭代法x=0;a=0;b=1;ci=0;while abs(b-a)>5*1e-4a=x;x=x-(exp(x)+10*x-2)/(exp(x)+10);b=x;ci=ci+1;endfprintf('牛顿迭代法求得 x 的值是 %f\t',x); fprintf('迭代次数是 %d\n',ci); %（5）分析绝对误差solve('exp(x)+10*x-2=0'); fprintf('x 的真值是%.8f',x);运行结果：二分法求得 c 的值是 0.090332 迭代次数是 11 绝对为误差为0.000193 不动点迭代求得 x 的值是 0.090513 迭代次数是 4 绝对为误差为0.000012 艾特肯加速迭代求得 x 的值是 0.099488 迭代次数是 3 绝对为误差为0.008963 牛顿迭代法求得 x 的值是 0.090525 迭代次数是 2 绝对为误差为0.000000 x 的真值是0.090525>>分析可知加速迭代绝对误差最大，二分法次之，牛顿法和不动点法迭代效果最好。
东南大学计算方法上机报告实验报告完整版
实习题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 的逐渐增大，两种算法的运行结果相差越来越大。
计算方法与实习上机实验报告
计算方法与实习上机实验报告一、引言本文旨在介绍和展示我们在“计算方法”课程中的实习上机实验环节所完成的一些关键任务和所取得的成果。
该实验课程的目标是让我们更深入地理解和应用各种计算方法，并在实际操作中提高我们的编程和问题解决能力。
二、实验内容与目标实验的主要内容是利用各种计算方法解决实际数学问题。
我们被要求使用编程语言（如Python或Java）来实现和解决这些问题。
这些问题包括使用牛顿法求解平方根，使用蒙特卡洛方法计算圆周率，以及使用最优化方法求解函数的最小值等。
实验的目标不仅是让我们掌握计算方法的基本理论，更是要让我们能够在实际操作中运用这些方法。
我们需要在实习过程中，通过与同伴们合作，共同解决问题，提高我们的团队合作能力和问题解决能力。
三、实验过程与问题解决策略在实验过程中，我们遇到了许多问题，如编程错误、理解困难和时间压力等。
我们通过相互讨论、查阅资料和寻求教师帮助等方式，成功地解决了这些问题。
例如，在实现牛顿法求解平方根时，我们一开始对导数的计算和理解出现了一些错误。
但我们通过查阅相关资料和讨论，最终理解了导数的正确计算方法，并成功地实现了牛顿法。
四、实验结果与结论通过这次实习上机实验，我们不仅深入理解了计算方法的基本理论，还在实际操作中提高了我们的编程和问题解决能力。
我们的成果包括编写出了能有效求解平方根、计算圆周率和求解函数最小值的程序。
这次实习上机实验非常成功。
我们的团队不仅在理论学习和实践操作上取得了显著的进步，还在团队合作和问题解决方面积累了宝贵的经验。
这次实验使我们对计算方法有了更深的理解和认识，也提高了我们的编程技能和解决问题的能力。
五、反思与展望回顾这次实验，我们意识到在实验过程中，我们需要更好地管理我们的时间和压力。
在解决问题时，我们需要更有效地利用我们的知识和资源。
在未来，我们希望能够更加熟练地运用计算方法，并能够更有效地解决问题。
我们也希望能够将所学的计算方法应用到更广泛的领域中，如数据分析、科学研究和工业生产等。
计算方法第七章上机报告
实验报告名称求解常微分方程班级：020991 学号：02099037 姓名：杜凡成绩：1实验目的1）熟悉求解常微分方程初值问题的有关方法和理论，主要是改进欧拉法、四阶龙格-库塔法与阿当姆斯方法。
2）会变质上述方法的计算程序，包括求解常微分方程组的计算程序。
3）通过对各种求解方法的计算实习，体会各种解法的功能、优缺点及适用场合，会选取适当的求解方法。
2 实验内容实习题7.1用改进欧拉法与四阶龙格-库塔公式求解所给微分方程初值问题；7.2 用四阶龙格-库塔公式解下列微分方程初值问题；7.3用阿当姆斯方法解微分方程初值问题；3实验步骤7.11）根据改进欧拉法的算法编写改进欧拉法求微分方程的函数// 实验环境的配置，例如添加什么函数，库，头文件等，以及你的思路都可以写。
3 程序设计// 程序流程图、代码。
以下均用matlab编写1）改进欧拉法function Heun2(f,a,b,y0,n)h=(b-a)/n;x=a:h:b;%ytrue=f1(-1*x);y=y0*ones(1,n+1);for j=2:n+1yp=y(j-1)+h*f(x(j-1),y(j-1));yc=y(j-1)+h*f(x(j),yp);y(j)=(yp+yc)/2;endfor i=1:n+1fprintf('x[%d]=%f\t y[%d]=%f\n',i-1,x(i),i-1,y(i));%fprintf('x[%d]=%f\t y[%d]=%f\t ytrue[%d]=%f\n',i-1,x(i),i-1,y(i),i-1,%ytrue(i));end4实验结果及分析// 程序运行的结果，可以添加截图以说明问题。
7.11）改进欧拉法2）四阶龙格库塔公式解方程组3）阿当姆斯方法解方程//实验结果分析，包括误差分析和结论。
2）实验结果分析改进欧拉公式的局部截断误差O（h^3），h=0.1,则绝对误差e<1.0*10^2.四阶龙格库塔方法的局部截断误差为O(h^5),h=0.1,则绝对误差e<1.0*10^4.阿当姆斯方法的局部截断误差为O(h^5),h=0.1,则绝对误差e<1.0*10^4.5总结// 通过本实验掌握的内容，以及在实验中遇到的问题及解决方法。
哈工大计算方法上机实验报告
实验报告一题目： 非线性方程求解摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。
本实验采用两种常见的求解方法二分法和Newton 法及改进的Newton 法。
前言：（目的和意义）掌握二分法与Newton 法的基本原理和应用。
数学原理：对于一个非线性方程的数值解法很多。
在此介绍两种最常见的方法：二分法和Newton 法。
对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b ]上连续，f(a)f(b)<0，且f(x)在[a,b ]内仅有一个实根x *，取区间中点c ，若，则c 恰为其根，否则根据f(a)f(c)<0是否成立判断根在区间[a,c ]和[c,b ]中的哪一个，从而得出新区间，仍称为[a,b ]。
重复运行计算，直至满足精度为止。
这就是二分法的计算思想。
Newton 法通常预先要给出一个猜测初值x 0，然后根据其迭代公式)()('1k k k k x f x f x x -=+ 产生逼近解x *的迭代数列{x k }，这就是Newton 法的思想。
当x 0接近x *时收敛很快，但是当x 0选择不好时，可能会发散，因此初值的选取很重要。
另外，若将该迭代公式改进为)()('1k k k k x f x f r x x -=+ 其中r 为要求的方程的根的重数，这就是改进的Newton 法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton 法快的多。
程序设计：本实验采用Matlab 的M 文件编写。
其中待求解的方程写成function 的方式，如下 function y=f(x);y=-x*x-sin(x);写成如上形式即可，下面给出主程序。
二分法源程序：clear%%%给定求解区间b=1.5;a=0;%%%误差R=1;k=0;%迭代次数初值while (R>5e-6) ;c=(a+b)/2;if f12(a)*f12(c)>0;a=c;elseb=c;endR=b-a;%求出误差k=k+1;endx=c%给出解Newton法及改进的Newton法源程序：clear%%%% 输入函数f=input('请输入需要求解函数>>','s')%%%求解f(x)的导数df=diff(f);%%%改进常数或重根数miu=2;%%%初始值x0x0=input('input initial value x0>>');k=0;%迭代次数max=100;%最大迭代次数R=eval(subs(f,'x0','x'));%求解f(x0)，以确定初值x0时否就是解while (abs(R)>1e-8)x1=x0-miu*eval(subs(f,'x0','x'))/eval(subs(df,'x0','x'));R=x1-x0;x0=x1;k=k+1;if (eval(subs(f,'x0','x'))<1e-10);breakendif k>max ;%如果迭代次数大于给定值，认为迭代不收敛，重新输入初值 ss=input('maybe result is error,choose a new x0,y/n?>>','s');if strcmp(ss,'y')x0=input('input initial value x0>>');k=0;elsebreakendendendk;%给出迭代次数x=x0；%给出解结果分析和讨论：1. 用二分法计算方程02sin 2=-x x 在[1，2]内的根。
实用数值计算方法上机实验报告讲解
实用数值计算方法上机实验报告学院：化学工程学院姓名：陶明专业：工业催化学号： 21113011681. 问题来源某公司饲养实验用的动物以供出售,已知这些动物的生长对饲料中3种营养成分（蛋白质,矿物质和维生素）特别敏感,每个动物每周至少需要蛋白质60g,矿物质3g,维生素8mg,该公司能买到5种不同的饲料,每种饲料1kg 所含各种营养成分和成本如表1所示,如果每个小动物每周食用饲料不超过52kg,求既满足动物生长需要,又能使总成本最低的饲料配方。
数学模型 设需要饲料A1,A2,A3,A4,A5分别为x1,x2,x3,x4,x5（单位kg ）12345min 0.20.70.40.30.5S x x x x x =++++1234512345123451234512350.3x +2x +x +0.6x +1.8x 600.1x +0.05x +0.02x +0.2x +0.05x 3.0.05x +0.1x +0.02x +0.2x +0.08x 8x +x +x +x +x 52,,,4,0s t x x x x x ≥⎧⎪≥⎪⎪≥⎨⎪≤⎪⎪≥⎩在LINGO 的MODEL 窗口内输入如下模型：Min =0.2*x1+0.7*x2+0.4*x3+0.3*x4+0.5*x5; 0.3*x1+2*x2+x3+0.6*x4+1.8*x5>60;0.1*x1+0.05*x2+0.02*x3+0.2*x4+0.05*x5>3; 0.05*x1+0.1*x2+0.02*x3+0.2*x4+0.08*x5>8; x1+x2+x3+x4+x5<52; end求解输出结果如下：Global optimal solution found.Objective value: 22.40000Infeasibilities: 0.000000Total solver iterations: 3Variable Value Reduced Cost X1 0.000000 0.7000000 X2 12.00000 0.000000 X3 0.000000 0.6166667 X4 30.00000 0.000000 X5 10.00000 0.000000 Row Slack or Surplus Dual Price1 22.40000 -1.0000002 0.000000 -0.58333333 4.100000 0.0000004 0.000000 -4.1666675 0.000000 0.8833333结果分析：因此每周每个动物的配料为饲料A2,A4,A5分别为12kg,30kg,10kg,可使得成本达到最低,最低成本为22.4元。
计算机上机报告
计算方法上机实验报告上课时间: 2014-2015学年秋学期, 6～14周一. 拉格朗日插值------------------------------------------------------1二. 牛顿插值------------------------------------------------------------3三. 改进欧拉法---------------------------------------------------------5四. 四阶龙格-库塔-----------------------------------------------------7五. 牛顿迭代------------------------------------------------------------9 六．复化Simpson公式------------------------------------------------11七. Romberg算法------------------------------------------------------14八. Seidel迭代法------------------------------------------------------17九. Gauss列主元消去法----------------------------------------------20一.拉格朗日插值1.程序代码#include<iostream.h>void Lagrange(){int i=0;double a[10],b[10],L,L1,L2,L3,L4,x;cout<<"x=";for(i=0;i<4;i++){cin>>a[i];}cout<<"y=";for(i=0;i<4;i++){cin>>b[i];}cout<<"x=";cin>>x;L1=(x-a[1])*(x-a[2])*(x-a[3])*b[0]/(a[0]-a[1])/(a[0]-a[2])/(a[0]-a[3]);L2=(x-a[0])*(x-a[2])*(x-a[3])*b[1]/(a[1]-a[0])/(a[1]-a[2])/(a[1]-a[3]);L3=(x-a[0])*(x-a[1])*(x-a[3])*b[2]/(a[2]-a[0])/(a[2]-a[1])/(a[2]-a[3]);L4=(x-a[0])*(x-a[1])*(x-a[2])*b[3]/(a[3]-a[0])/(a[3]-a[1])/(a[3]-a[2]);L=L1+L2+L3+L4;cout<<"L="<<L;}void main(){Lagrange();cout<<endl;}2.例子3.运行结果二. 牛顿插值1.程序代码#include <iostream.h>#include<conio.h>void main(){int n,i,j;double A[50][50],*x,*y;cout<<"请输入插值节点数: ";cin>>n;x=new double[n];y=new double[n];cout<<"请输入这"<<n<<"个插值节点（xi,yi）: "<<endl;for(i=0;i<=n-1;i++)cin>>x[i]>>y[i];double K=1,xx,N=0,P;for(i=0;i<=n-1;i++){A[i][0]=x[i];A[i][1]=y[i];}for(j=2;j<=n;j++)for(i=1;i<=n-1;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-1;i++)cout<<"输出第"<<i+1<<"阶差商为: "<<A[i][i+1]<<endl;cout<<"请输入预求值x=";cin>>xx;for(i=0;i<n-1;i++){K*=xx-x[i];N+=A[i+1][i+2]*K;P=A[0][1]+N;}cout<<"Newton插值结果为: y="<<P<<endl;getch();}2.例子3.运行结果三. 改进欧拉法1.程序代码#include<iostream.h>#include<conio.h>double fun(double x,double y){return(-0.9*y/(1+2*x));}void main(){double a,b,*y,h,*x,yp,yc;int n,k;cout<<"常微分方程为y'=-0.9*y/(1+2*x)"<<endl;cout<<"其中0<=x<=1"<<endl;cout<<"初值为y(0)=1"<<endl;cout<<"请输入计算区间(a,b): ";cin>>a>>b;cout<<"请输入步长h: ";cin>>h;cout<<"请输入计算次数: ";cin>>n;y=new double[n];x=new double[n];cout<<"请输入初值y(0)=";cin>>y[0];x[0]=a;for(k=0;k<=n;k++){yp=y[k]+h*fun(x[k],y[k]);yc=y[k]+h*fun(x[k]+h,yp);y[k+1]=0.5*(yp+yc);x[k+1]=x[k]+h;}cout<<"迭代结果为: "<<endl;for(k=0;k<=n;k++)cout<<"x"<<k<<"="<<x[k]<<'\t'<<"y"<<k<<"="<<y[k]<<endl;getch();}2.例子3.运行结果四. 四阶龙格-库塔1.程序代码#include<iostream.h>#include<conio.h>double fun(double x,double y){return(x-y);}void main(){double a,b,*y,h,x,k1,k2,k3,k4;int n,k;cout<<"常微分方程为y'=x-y"<<endl;cout<<"其中0<=x<=1"<<endl;cout<<"初值为y(0)=0"<<endl;cout<<"请输入计算区间(a,b): ";cin>>a>>b;cout<<"请输入步长h: ";cin>>h;cout<<"请输入计算次数: ";cin>>n;y=new double[n];cout<<"请输入初值y(0): ";cin>>y[0];x=a;cout<<"Runge-Kutta迭代法结果为: "<<endl;cout<<"x0="<<x<<'\t'<<"y0="<<y[0]<<endl;for(k=0;k<=n-1;k++){k1=fun(x,y[k]);k2=fun(x+h/2,y[k]+k1*h/2);k3=fun(x+h/2,y[k]+k2*h/2);k4=fun(x+h,y[k]+k3*h);y[k+1]=y[k]+(h/6)*(k1+2*(k2+k3)+k4);cout<<"x"<<k+1<<"="<<x+h<<'\t'<<"y"<<k+1<<"="<<y[k+1]<<endl;x=x+h;}getch();}2.例子3.运行结果五. 牛顿迭代法1.程序代码(C++代码)#include<iostream>#include<cmath>using namespace std;double newtondiedai(double a,double b,double c,double d,double x); int main(){double a,b,c,d;double x=1.5;cout<<"请依次输入方程四个系数: ";cin>>a>>b>>c>>d;x=newtondiedai(a,b,c,d,x);cout<<x<<endl;return 0;}double newtondiedai(double a,double b,double c,double d,double x) {while(abs(a*x*x*x+b*x*x+c*x+d)>0.000001){x=x-(a*x*x*x+b*x*x+c*x+d)/(3*a*x*x+2*b*x+c);}return x;}2.例子3.运行结果六. 复化Simpson公式1.程序代码（C++代码）#include<iostream.h>#include<math.h>double function1(double x)//被积函数{double s;s=x/(4+x*x);return s;}double function2(double x)//被积函数{double s;s=sqrt(x);return s;}double ReiterationOfSimpson(double a,double b,double n,double f(double x))//复化Simpson公式{double h,fa,fb,xk,xj;h=(b-a)/n;fa=f(a);fb=f(b);double s1=0.0;double s2=0.0;for(int k=1;k<n;k++){xk=a+k*h;s1=s1+f(xk);}for(int j=0;j<n;j++){xj=a+(j+0.5)*h;s2=s2+f(xj);}double sn;//和sn=h/6*(fa+fb+2*s1+4*s2);//复化Simpson公式return sn;}main(){double a,b,Result,n;cout<<"请输入积分下限:"<<endl;cin>>a;cout<<"请输入积分上限:"<<endl;cin>>b;cout<<"请输入分割区间数n:"<<endl;cin>>n;cout<<"复化Simpson 公式计算结果:";Result=ReiterationOfSimpson(a, b, n,function1); cout<<Result<<endl;}2.例子 （620(n 3)4x I dx x ==+⎰）3.运行结果七. Romberg算法1.程序代码（C++代码）#include<iostream>#include<cmath>using namespace std;#define f(x) (4/(1+x*x))#define epsilon 0.0001#define MAXREPT 10double Romberg(double aa,double bb) { int m,n;double h,x;double s,q;double ep;double *y =new double[MAXREPT]; double p;h=bb-aa;y[0]=h*(f(aa)+f(bb))/2.0;m=1;n=1;ep=epsilon+1.0;while((ep>=epsilon)&&(m<MAXREPT)) { p =0.0;for(int i=0;i<n;i++){ x=aa+(i+0.5)*h;p=p+f(x);}p=(y[0] + h*p)/2.0;s=1.0;for(int k=1;k<=m;k++){ s=4.0*s;q=(s*p-y[k-1])/(s-1.0); y[k-1]=p;p=q;}p=fabs(q-y[m-1]);m=m+1;y[m-1]=q;n=n+n;h=h/2.0;}return (q);}int main(){double a,b;cout<<"Romberg积分,请输入积分范围a,b:"<<endl; cin>>a>>b;cout<<"积分结果:"<<Romberg(a,b)<<endl;system("pause");return 0;}2.例子（1241I dxx=+⎰）3.运行结果八. Seidel迭代法1.程序代码（C++代码）# include <math.h># include <stdio.h># define max 100# define EPS 1e-6float a[3][3]={{10,-1,-2},{-1,10,-2},{-1,-1,5}}; float b[3]={7.2,8.3,4.2};float x[3]={0,0,0};float y[3];float S(int m){int n;float S=0;float y;for(n=0;n<3;n++){if(m==n){}else{S+=a[m][n]*x[n];}}y=(b[m]-S)/a[m][m];return y;}void main(){int i;int F,T=1,k=1;do{F=0;if(T){T=0;}else{k++;for(i=0;i<3;i++){x[i]=y[i];}}for(i=0;i<3;i++){y[i]=S(i);}for(i=0;i<3;i++){if(fabs(x[i]-y[i])>EPS){F=1;}}printf("%d\n",k);}while(((F==1)&&k!=max));printf("迭代次数:%d\n",k);for(i=0;i<3;i++){printf("x[%d]=%1.6f\n",i,y[i]); }}2.例子（1012237.21102238.31253 4.2x x xx x xx x x--=⎧⎪-+-=⎨⎪--+=⎩）3.运行结果九. Gauss列主元消去法1.程序代码（C++代码）#include <stdio.h>#include<conio.h>#include <math.h>#define max_dimension 20int n;static float a[max_dimension][max_dimension]; static float b[max_dimension];static float x[max_dimension];void main(){int i;int j;int d;int row;float temp;float known_items;float l[max_dimension][max_dimension];printf("请输入阶数:");scanf("%d",&n);printf("\n");printf("请输入系数矩阵的值: ");printf("\n");for(i=0; i<n; i++){ printf("输入第%d行的值:",i+1);for (j=0; j<n; j++){scanf("%f",&a[i][j]);}printf("\n");}printf("请输入常数项的值: ");for(i=0; i<n; i++)scanf("%f",&b[i]);printf("\n");for(d=0; d<n-1; d++){row=d;for(i=d+1; i<n; i++){if(fabs(a[i][d])>fabs(a[row][d]))row=i;}if(row!=d){for(j=d; j<n; j++){temp=a[row][j];a[row][j]=a[d][j];a[d][j]=temp;}temp=b[row];b[row]=b[d];b[d]=temp;}for(i=d+1; i<n; i++){l[i][d]=-a[i][d]/a[d][d];for (j=d; j<n; j++){a[i][j]=a[i][j]+a[d][j]*l[i][d];}b[i]=b[i]+b[d]*l[i][d];}}printf("\n");for (i=n-1; i>-1; i--){known_items=0;for(j=1; j<n-i; j++){known_items=known_items+a[i][i+j]*x[i+j]; }x[i]=(b[i]-known_items)/a[i][i];}printf("方程组的根为:\n\n");for(i=0; i<n; i++)printf("%.5f ",x[i]);printf("\n\n");getch();}2.例子3.运行结果。
西安交通大学 计算方法上机报告
计算方法上机报告姓名：学号：班级：机械硕4002上课班级：02班说明：本次上机实验使用的编程语言是Matlab 语言，编译环境为MATLAB 7.11.0，运行平台为Windows 7。
1. 对以下和式计算：∑∞⎪⎭⎫ ⎝⎛+-+-+-+=0681581482184161n n n n S n，要求：① 若只需保留11个有效数字，该如何进行计算； ② 若要保留30个有效数字，则又将如何进行计算；（1） 算法思想1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为： 142111416818485861681n n n a n n n n n ε⎛⎫=---<< ⎪+++++⎝⎭； 2、为了保证计算结果的准确性，写程序时，从后向前计算； 3、使用Matlab 时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数)（2）算法结构1. ;0=s⎪⎭⎫⎝⎛+-+-+-+=681581482184161n n n n t n； 2. for 0,1,2,,n i =⋅⋅⋅ if 10m t -≤end;3. for ,1,2,,0n i i i =--⋅⋅⋅;s s t =+（3）Matlab源程序clear; %清除工作空间变量clc; %清除命令窗口命令m=input('请输入有效数字的位数m='); %输入有效数字的位数s=0;for n=0:50t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));if t<=10^(-m) %判断通项与精度的关系break;endend;fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值for i=n-1:-1:0t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));s=s+t; %求和运算ends=vpa(s,m) %控制s的精度（4）结果与分析当保留11位有效数字时，需要将n值加到n=7，s =3.1415926536；当保留30位有效数字时，需要将n值加到n=22,s =3.14159265358979323846264338328。
数值计算方法上机实验报告
数值计算方法上机实验报告实验目的：复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。
利用计算机语言独立编写、调试数值计算方法程序，培养学生利用计算机和所学理论知识分析解决实际问题的能力。
上机练习任务：利用计算机基本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 。
华科计算方法上机实验报告
《计算方法》实验报告专业及班级 : 姓名：学号：日期：2013年6月12日一，方程求根1.用牛顿迭代法求解下列方程的正根（1）log(1+x)-x^2=0 x=7.468818e-001求解代码：牛顿迭代法的函数M文件：function x=Newt_n(f_name,x0)x=x0;xb=x+1;k=0;n=50;del_x=0.01;while abs(x-xb)>0.000001k=k+1;xb=x;if k>=n break;end;y=feval(f_name,x);y_driv=(feval(f_name,x+del_x)-y)/del_x;x=xb-y/y_driv;fprintf('k=%3.0f,x=%12.5e,y=%12.5e,yd=%12.5e\n',k,x,y,y_driv); end;fprintf('\n Final answer=%12.6e\n',x);待求方程的函数M文件：function y=eqn_1(x)y=log(x+1)-x^2;运行结果：Newt_n('eqn_1',1)k= 1,x=7.96954e-001,y=-3.06853e-001,yd=-1.51125e+000k= 2,x=7.50200e-001,y=-4.90424e-002,yd=-1.04895e+000k= 3,x=7.46936e-001,y=-3.07003e-003,yd=-9.40663e-001k= 4,x=7.46882e-001,y=-5.03365e-005,yd=-9.33074e-001k= 5,x=7.46882e-001,y=-6.30906e-007,yd=-9.32949e-001Final answer=7.468818e-001（2）e^x-5x=0 x=6.052671e-001待解方程：function y=eqn_1(x)y=exp(x)-5*x^2;运行结果：Newt_n('eqn_2',2)k= 1,x=1.00102e+000,y=-1.26109e+001,yd=-1.26239e+001 k= 2,x=6.88531e-001,y=-2.28918e+000,yd=-7.32552e+000 k= 3,x=6.11607e-001,y=-3.79584e-001,yd=-4.93453e+000 k= 4,x=6.05365e-001,y=-2.69226e-002,yd=-4.31343e+000 k= 5,x=6.05268e-001,y=-4.13312e-004,yd=-4.26254e+000 k= 6,x=6.05267e-001,y=-3.99547e-006,yd=-4.26175e+000 Final answer=6.052671e-001ans =0.6053（3）x^3+2x-1=0 x=4.533977e-001待解方程：function y=eqn_1(x)y=x^3+2*x-1;运行结果：Newt_n('eqn_3',3)k= 1,x=1.89997e+000,y=3.20000e+001,yd=2.90901e+001k= 2,x=1.15047e+000,y=9.65861e+000,yd=1.28868e+001k= 3,x=6.80277e-001,y=2.82368e+000,yd=6.00536e+000k= 4,x=4.82154e-001,y=6.75369e-001,yd=3.40884e+000k= 5,x=4.53984e-001,y=7.63945e-002,yd=2.71198e+000k= 6,x=4.53401e-001,y=1.53570e-003,yd=2.63202e+000k= 7,x=4.53398e-001,y=8.46834e-006,yd=2.63042e+000k= 8,x=4.53398e-001,y=4.41262e-008,yd=2.63041e+000 Final answer=4.533977e-001ans =0.4534（4）(x+2)^0.5-x=0 x=2.0000待解方程：function y=eqn_1(x)y=(x+2)^0.5-x;运行结果：Newt_n('eqn_4',1)k= 1,x=2.02879e+000,y=7.32051e-001,yd=-7.11565e-001 k= 2,x=2.00002e+000,y=-2.16052e-002,yd=-7.51049e-001 k= 3,x=2.00000e+000,y=-1.72788e-005,yd=-7.50157e-001 k= 4,x=2.00000e+000,y=-3.60276e-009,yd=-7.50156e-001 Final answer=2.000000e+000ans =2.00002.先用图解法确定初始点，然后再求方程 x sin x −1 = 0 x ∈[0,10]的所有根。
计算方法上机报告
计算方法上机作业计算方法上机作业1.1 共轭梯度法1.1.1 方法描述共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次泛函极小值问题，即：**=min ⎧Φ⎪⎨ΦΦ⎪⎩求，使得（）（） (0-1)其中，泛函1(,)(,)2K F ΦΦΦ-Φ（）=(0-2)共轭梯度法在形式上迭代法的特征(k 1)(k)()k k d α+Φ=Φ+只要确定了初始向量(0)Φ，然后确定确定最佳步长k α和搜索方向()k d 就可以得出迭代格式。
给定初始向量(0)Φ后，由于负梯度方向是下降最快的方向，所以第一次迭代的搜索方向(0)(0)(0)(0)d r F K =-∇Φ-Φ（）==，得(1)(0)(0)0d αΦ=Φ+；第二次迭代时，选取的搜索方向不再是 (1)r ，而是(1)(1)(0)0d r d β=+使得(1)d 与(0)d 是关于矩阵K 的共轭向量，即(1)(0)(,)0d Kd= ，由此可以求得参数(1)(0)0(0)T (0)T r Kd d Kdβ=-，然后沿方向(1)d 进行搜索得(2)(1)(1)1d αΦ=Φ+。
最佳步长k α的计算是在给定迭代点()k Φ以及搜索方向()k d 后，要求选取非负实数使得(k)()k k d αΦ+（）达到最小，即(k)()min ()=k f d ααΦ+（）对α求导，令'(k)()()()0k T k f d d αα=∇Φ+=（），得：()()()()()()()()0k k T k k k T k K b Kd d r Kd d ααΦ-+=-+=所以()()()()k T k k k T k r d d Kdα=由此确定了每一步的最佳步长。
实际上,k k αβ的计算还可以化简，在此不再推导，化简后的表达式在下列算法中。
1.1.2 程序流程程序分为两个部分，主程序读取或产生矩阵和右端向量A,b 以及初始值x0和误差限和最大迭代次数tol 、Imax ，然后调用方法部分进行计算得出x ，方法部分放在iter_solver 模块中，调用函数接口为call CG_method(A,b,x,x0,tol,Imax)。
计算方法上机报告
计算方法上机报告硕5028李兰鑫31150350051. 计算以下和式：0142118184858616nn S n n n n ∞=⎛⎫=--- ⎪++++⎝⎭∑，要求： (1)若保留11个有效数字，给出计算结果，并评价计算的算法； (2)若要保留30个有效数字，则又将如何进行计算。
1.1算法思想由题可得S<(1/16^n)*(4/(8*n+1))保留11个有效数字，即误差ε<0.5*10^（1-11） 保留30个有效数字，即误差ε<0.5*10^（1-30） 1.2源程序%有效数字为11的情况digita=11;（第二问中改为30即可） for i=1:1:1000S0=(1/(16^(i-1)))*(4/(8*i-7)); S1=(1/16^i)*(4/(8*i+1)); if ((S0-S1)<0.5*10^(1-digita)) n=i; break ; end end S=0;digits(digita);for i=0:1:nS=S+vpa(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); end1.3运行结果2.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。
已探测(1)请用合适的曲线拟合所测数据点；(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；2.1算法思想由题可知应使用多项式插值，考虑到本例中的拟合数据点较多，使用多项式插值时会因龙格现象造成误差较大，故决定采用三次样条插值。
2.2源程序%------三次样条插值多项式------clc;clear;%-----插值函数和插值点-------n=20;%插值点个数11i=0:n;%子元素从0到nx=10*i;f=[9.01 8.96 7.96 7.97 8.02 9.05 10.13 11.18 12.26 13.28 13.32 12.61 11.29 10.229.15 7.90 7.95 8.86 9.81 10.80 10.93];%--输出插值点---fprintf('当n=%d时\n',n)fprintf('插值点x坐标为：\n')disp(x);fprintf('插值点y坐标f（x）为：\n')disp(f);%----计算差商表，并将有用的差商进行选择和保存-----y=f;%存放差商d=zeros(1,n+1);%存放三弯矩方程组的右端项for i=1:n%计算第i阶差商for j=n+1:-1:i+1y(j)=(y(j)-y(j-1))/(x(j)-x(j-i));endif i==2%将第2阶差商赋值给dd(2:n)=6.*y(3:n+1);endif i==3%三次样条插值的第三类边界条件d(1)=-12*(x(2)-x(1))*y(i+1);d(n+1)=12*(x(n+1)-x(n))*y(n+1); endend%-------用追赶法求三弯矩方程组------a=ones(1,n+1);%存储下三角的非零元素b=2*ones(1,n+1);%存储对角元素c=ones(1,n+1);%存储上三角的非零元素for i=2:na(i)=(x(i)-x(i-1))/(x(i+1)-x(i-1));c(i)=1-a(i);enda(n+1)=-2;c(1)=-2;%-------计算矩阵L和U-----------%初始化L和Ul=zeros(1,n+1);u=zeros(1,n+1);r=zeros(1,n+1);u(1)=b(1);r(1)=c(1);for k=2:n+1l(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*r(k-1);r(k)=c(k);end%--------使用追赶法进行求解----------M=zeros(1,n+1);%方程组的解M(1)=d(1);for k=2:n+1M(k)=d(k)-l(k)*M(k-1);endM(n+1)=M(n+1)/u(n+1);for k=n:-1:1M(k)=(M(k)-r(k)*M(k+1))/u(k);endfprintf('三弯矩方程组的解M为：\n')disp(M);%------计算f（x）、N(x)和S(x)-------Sy1=zeros(1,200);x1=zeros(1,200);for i=1:200x1(i)=i-1;endfor i=1:200for j=1:nif x1(i)>=x(j)&&x1(i)<=x(j+1)temp=x(j+1)-x(j);Sy1(i)=M(j).*(x(j+1)-x1(i)).^3/6/temp+M(j+1).*(-x(j)+x1(i)).^3/6/temp+(f(j)-M(j)*te mp^2/6).*(x(j+1)-x1(i))/temp+(f(j+1)-M(j+1)*temp^2/6).*(-x(j)+x1(i))/temp;endendendfprintf('由三次样条插值方法计算的各年的产量为：\n')disp(Sy1)plot(x,f,'*')%画出由插值点构成的曲线hold onplot(x1,Sy1,'r')%用红线画出三次样条插值拟合的曲线hold onlegend('三次样条插值');%------求解电缆长度------Cordinate=[0:20/(length(Sy1)-1):20;Sy1];diff_Cordinate=[diff(Cordinate(1,:));diff(Cordinate(2,:))];disp('电缆长度为：');anwser=sum(sqrt(diff_Cordinate(1,:).*diff_Cordinate(1,:)+diff(Cordinate(2,:)).*diff( Cordinate(2,:))))2.3运行结果由图可知，三次样条插值拟合数据点比较好，没有龙格现象。
数值计算方法上机实验报告
（2）计算机程序框图：（见下页）
（3）输入变量、输出变量说明：
输入变量： 系数矩阵元素， 常向量元素
称为改进欧拉公式。
（2）计算机程序框图：（见下页）
（3）输入变量、输出变量说明：
输入变量： 处置点， 区间长度， 计算次数
输出变量： 初值问题的数值解法结果
（4）具体算例及求解结果：
例：求解初值问题（课本P242例7-2）
求解结果：
0.1
1.095909
1.095909
0.6
1.485956
1.485955
输出变量： 解向量元素
（4）具体算例及求解结果：
例：用列选主元法求解下列线性方程组（课本P65例3-3）
求解结果：
3、 分解法求解线性方程组
（1）算法原理：
求解线性方程组 时，当对 进行 分解，则等价于求解 ，这时可归结为利用递推计算相继求解两个三角形（系数矩阵为三角矩阵）方程组，用顺代，由
求出 ，再利用回带，由 求出 。
例：根据给定的函数 的实例数据表，试用最小二乘法求二次拟合多项式。（课本P186习题3）
求解结果：
6、变步长梯形求积分
（1）算法原理：
设将积分区间 分成 等份，即有 个子区间，分点 ，其中步长
对于子区间 ，利用体型求其积分近似值
对于子区间 有
对于子区间 再取其中点
作新节点，此时区间数增加了一倍为 ，
0.2
1.184097
计算方法实验上机报告(完整版)
简单迭代法#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];}。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验报告名称求解常微分方程
班级：020991 学号：02099037 姓名：杜凡成绩：
1实验目的
1）熟悉求解常微分方程初值问题的有关方法和理论，主要是改进欧拉法、四阶龙格-库塔法与阿当姆斯方法。
2）会变质上述方法的计算程序，包括求解常微分方程组的计算程序。
3）通过对各种求解方法的计算实习，体会各种解法的功能、优缺点及适用场合，会选取适当的求解方法。
2 实验内容
实习题7.1用改进欧拉法与四阶龙格-库塔公式求解所给微分方程初值问题；
7.2 用四阶龙格-库塔公式解下列微分方程初值问题；
7.3用阿当姆斯方法解微分方程初值问题；
3实验步骤
7.1
1）根据改进欧拉法的算法编写改进欧拉法求微分方程的函数
// 实验环境的配置，例如添加什么函数，库，头文件等，以及你的思路都可以写。
3 程序设计
// 程序流程图、代码。
以下均用matlab编写
1）改进欧拉法
function Heun2(f,a,b,y0,n)
h=(b-a)/n;
x=a:h:b;
%ytrue=f1(-1*x);
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)=(yp+yc)/2;
end
for i=1:n+1
fprintf('x[%d]=%f\t y[%d]=%f\n',i-1,x(i),i-1,y(i));
%fprintf('x[%d]=%f\t y[%d]=%f\t ytrue[%d]=%f\n',i-1,x(i),i-1,y(i),i-1,%ytrue(i));
end
4实验结果及分析
// 程序运行的结果，可以添加截图以说明问题。
7.1
1）改进欧拉法
2）四阶龙格库塔公式解方程组3）阿当姆斯方法解方程
//实验结果分析，包括误差分析和结论。
2）实验结果分析
改进欧拉公式的局部截断误差O（h^3），h=0.1,则绝对误差e<1.0*10^2.
四阶龙格库塔方法的局部截断误差为O(h^5),h=0.1,则绝对误差e<1.0*10^4.
阿当姆斯方法的局部截断误差为O(h^5),h=0.1,则绝对误差e<1.0*10^4.
5总结
// 通过本实验掌握的内容，以及在实验中遇到的问题及解决方法。
通过本实验，掌握了求解常微分方程的几种方法，包括改进欧拉法，四阶龙格库塔方法与阿当姆斯方法，并且学会了编制上述方法的程序，了解了各种解法的适用场合。
6参考资料
// 学习相关理论、编写程序及为了完成实验查阅的书籍和文献
// 英文参考文献格式
// 期刊
// [序号] 主要责任者. 文献题名[J]. 刊名, 年, 卷(期): 起止页码.
// 专著、论文集、学位论文、报告
// [序号] 主要责任者. 文献题名[文献类型标识]. 出版地: 出版者, 出版年. 起止页码.。