东南大学数值分析上机解剖

数值分析上机实验报告目录1.chapter1舍入误差及有效数 (1)2.chapter2Newton迭代法 (3)3.chapter3线性代数方程组数值解法-列主元Gauss消去法 (7)4.chapter3线性代数方程组数值解法-逐次超松弛迭代法 (8)5.chapter4多项式插值与函数最佳逼近 (10)1.chapter1舍入误差及有效数1.1题目设S N =∑1j 2−1N j=2，其精确值为)11123(21+--N N 。

（1）编制按从大到小的顺序11131121222-+⋯⋯+-+-=N S N ，计算S N 的通用程序。

（2）编制按从小到大的顺序1211)1(111222-+⋯⋯+--+-=N N S N ，计算S N 的通用程序。

（3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。

（编制程序时用单精度）（4）通过本次上机题，你明白了什么？ 1.2编写相应的matlab 程序 clear;N=input('please input N:'); AValue=((3/2-1/N-1/(N+1))/2); sn1=single(0); sn2=single(0); for i=2:Nsn1=sn1+1/(i*i-1); %从大到小相加的通用程序% endep1=abs(sn1-AValue); for j=N:-1:2sn2=sn2+1/(j*j-1); %从小到大相加的通用程序% endep2=abs(sn2-AValue);fprintf('精确值为:%f\n',AValue);fprintf('从大到小的顺序累加得sn=%f\n',sn1); fprintf('从大到小相加的误差ep1=%f\n',ep1); fprintf('从小到大的顺序累加得sn=%f\n',sn2); fprintf('从小到大相加的误差ep2=%f\n',ep2); disp('================================='); 1.3matlab 运行程序结果 >> chaper1please input N:100 精确值为:0.740050从大到小的顺序累加得sn=0.740049 从大到小相加的误差ep1=0.000001 从小到大的顺序累加得sn=0.740050 从小到大相加的误差ep2=0.000000 >> chaper1please input N:10000 精确值为:0.749900从大到小的顺序累加得sn=0.749852 从大到小相加的误差ep1=0.000048 从小到大的顺序累加得sn=0.749900 从小到大相加的误差ep2=0.000000please input N:1000000精确值为:0.749999从大到小的顺序累加得sn=0.749852 从大到小相加的误差ep1=0.000147 从小到大的顺序累加得sn=0.749999 从小到大相加的误差ep2=0.0000001.4结果分析以及感悟按照从大到小顺序相加的有效位数为：5,4,3。

实习题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 的逐渐增大，两种算法的运行结果相差越来越大。

第七章 偏微分方程数值解法——Crank-Nicolson 格式****（学号） *****（姓名）上机题目要求见教材P346,10题。

一、算法原理本文研究下列定解问题(抛物型方程)22(,) (0,0)(,0)() (0)(0,)(), (1,)() (0)u ua f x t x l t T t x u x x x l u t t u t t t T ϕαβ⎧∂∂-=<<≤≤⎪∂∂⎪=≤≤⎨⎪==<≤⎪⎩(1)的有限差分法，其中a 为正常数，,,,f ϕαβ为已知函数，且满足边界条件和初始条件。

关于式(1)的求解，采用离散化方法，剖分网格，构造差分格式。

其中，网格剖分是将区域{}0,0D x l t T =≤≤≤≤用两簇平行直线(0)(0)i k x x ih i M t t k k N τ==≤≤⎧⎨==≤≤⎩ 分割成矩形网格，其中,l Th M Nτ==分别为空间步长和时间步长。

将式(1)中的偏导数使用不同的差商代替，将得到不同的差分格式，如古典显格式、古典隐格式、Crank-Nicolson 格式等。

其中，Crank-Nicolson 格式具有更高的收敛阶数，应用更广泛，故本文采用Crank-Nicolson 格式求解抛物型方程。

Crank-Nicolson 格式推导：在节点(,)2i k x t τ+处考虑式(1)，有22(,)(,)(,)222i k i k i k u u x t a x t f x t t x τττ∂∂+-+=+∂∂ (2)对偏导数(,)2i k u x t t τ∂+∂用中心差分展开 []2311+131(,)(,)(,)(,) ()224k k i k i k i k i i k i k u ux t u x t u x t x t t t t ττηητ++∂∂+=--<<∂∂ (3) 将22(,)2i k u x t x τ∂+∂在节点(,)i k x t 和1(,)i k x t +表示为222+122224+1221(,)=(,)+(,)22 (,) ()8i k i k i k k k i i k i k u u ux t x t x t x x x ux t t x tττηη⎡⎤∂∂∂+⎢⎥∂∂∂⎣⎦∂-<<∂∂ (4)对以上两个偏导数用二阶差分展开[]2112224i+141(,)(,)2(,)(,) (,) ()12i k i k i k i k kk i k i i u x t u x t u x t u x t x hh u t x x x ξξ+-∂=-+∂∂-<<∂ (5)[]211111122241i+141(,)(,)2(,)(,) (,) ()12i k i k i k i k k ki k i i u x t u x t u x t u x t x hh u t x x xξξ++++-++∂=-+∂∂-<<∂ (6)将式(4)(5)(6)分别代入式(3)，略去高阶小量，用k i u 代替(,)i k u x t 并化简得()()()2111111112122,22k k k k k k k k i i i i i i i i i k a u u u u u u u u f x t h ττ+++++-+-⎛⎫⎡⎤---++-+=+ ⎪⎣⎦⎝⎭ (7) 令2/r a h τ=，将式(7)联合式(1)初始条件和边界条件，用矩阵的形式表示为：11112212211111221122221122221122 k k k k k kM M k k M M r r r r u u r r r r r ru u r r r r u u r r u u r r r r +++--+--⎡⎤⎡⎤+---⎢⎥⎢⎥⎢⎥⎢⎥⎡⎤⎡⎤⎢⎥⎢⎥-+--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+--⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥⎣⎦⎣⎦()()()()()()112211,22,2 ,2,22k k k k M k M k k k r f x t t t f x t f x t r f x t t t ττααττττττββ+--+⎡⎤⎛⎫+++ ⎪⎢⎥⎝⎭⎢⎥⎢⎥⎛⎫+ ⎪⎢⎥⎝⎭⎢⎥+⎢⎥⎢⎥⎛⎫⎢⎥+ ⎪⎢⎥⎝⎭⎢⎥⎛⎫⎢⎥+++ ⎪⎢⎥⎝⎭⎣⎦(8) Crank-Nicolson 格式的截断误差为22()R O h τ=+，具有较高的精度。

实用标准文案文档大全上机作业题报告2015.1.9 USER1.Chapter 11.1题目设S N =∑1j 2−1N j=2，其精确值为)11123(21+--N N 。

（1）编制按从大到小的顺序11131121222-+⋯⋯+-+-=N S N ，计算S N 的通用程序。

（2）编制按从小到大的顺序1211)1(111222-+⋯⋯+--+-=N N S N ，计算S N 的通用程序。

（3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。

（编制程序时用单精度） （4）通过本次上机题，你明白了什么？1.2程序1.3运行结果1.4结果分析按从大到小的顺序，有效位数分别为：6，4，3。

按从小到大的顺序，有效位数分别为：5，6，6。

可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。

当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。

因此，采取从小到大的顺序累加得到的结果更加精确。

2.Chapter 22.1题目（1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。

（2）给定方程03)(3=-=x xx f ,易知其有三个根3,0,3321=*=*-=*x x x○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收敛于根x2*。

试确定尽可能大的δ。

○2试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。

（3）通过本上机题，你明白了什么？2.2程序2.3运行结果（1）寻找最大的δ值。

算法为：将初值x0在从0开始不断累加搜索精度eps，带入Newton迭代公式，直到求得的根不再收敛于0为止，此时的x0值即为最大的sigma值。

运行Find.m，得到在不同的搜索精度下的最大sigma值。

数值分析上机报告XX：学号：2013年12月22日第四章38.（上机题）3次样条插值函数（1）编制求第一型3次样条插值函数的通用程序;端点条件为'0y =0.8，'10y =0.2。

用所编制程序求车门的3次样条插值函数S(x),并打印出S(i+0.5)(i=0,1,…9)。

解：（1）#include <iostream.h> #include <math.h>floatx1[100],f1[100],f2[99],f3[98],m[100],a[100][101],x,d[100]; float c[99],e[99],h[99],u[99],w[99],y_0,y_n,arr ,s; int i,j,k,n,q;void selectprint(float y) {if ((y>0)&&(y!=1)) cout<<"+"<<y; else if (y==1) cout<<"+"; else if (y<0) cout<<y; }void printY(float y){ if (y!=0) cout<<y; }float calculation(float x){ for (j=1;j<=n;j++) if (x<=x1[j]) {s=(float)(f1[j-1]+c[j-1]*(x-x1[j-1])+m[j-1]/2.0*(x-x1[j-1])*(x-x1[j-1])+e[j-1]*(x-x1[j-1])*(x-x1[j-1])*(x-x1[j-1])); break; }return s; }void main() {do{cout<<"请输入n值:";cin>>n;if ((n>100)||(n<1)) cout<<"请重新输入整数(1..100):"<<endl;}while ((n>100)||(n<1));cout<<"请输入Xi(i=0,1,...,"<<n<<"):";for (i=0;i<=n;i++) cin>>x1[i];cout<<endl;cout<<"请输入Yi(i=0,1,...,"<<n<<"n):";for (i=0;i<=n;i++) cin>>f1[i];cout<<endl;cout<<"请输入第一型边界条件Y'0，Y'n:";cin>>y_0>>y_n;cout<<endl;for (i=0;i<n;i++) h[i]=x1[i+1]-x1[i];for (i=1;i<n;i++) u[i]=(float) (h[i-1]/(h[i-1]+h[i]));for (i=1;i<n;i++) w[i]=(float) (1.0-u[i]);for (i=0;i<n;i++) f2[i]=(f1[i+1]-f1[i])/h[i]; //一阶差商for (i=0;i<n-1;i++) f3[i]=(f2[i+1]-f2[i])/(x1[i+2]-x1[i]); //二阶差商for (i=1;i<n;i++) d[i]=6*f3[i-1]; //求出所有的d[i]d[0]=6*(f2[0]-y_0)/h[0];d[n]=6*(y_n-f2[n-1])/h[n-1];for (i=0;i<=n;i++)for (j=0;j<=n;j++)if (i==j) a[i][j]=2;else a[i][j]=0;a[0][1]=1;a[n][n-1]=1;for (i=1;i<n;i++){a[i][i-1]=u[i];a[i][i+1]=w[i];}for (i=0;i<=n;i++) a[i][n+1]=d[i];for (i=1;i<=n;i++) //用追赶法解方程,得M[i]{arr=a[i][i-1];for (j=0;j<=n+1;j++)a[i][j]=a[i][j]-a[i-1][j]*arr/a[i-1][i-1];}m[n]=a[n][n+1]/a[n][n];for (i=n-1;i>=0;i--) m[i]=(a[i][n+1]-a[i][i+1]*m[i+1])/a[i][i];for (i=0;i<n;i++) //计算S（x）的表达式c[i]=(float) (f2[i]-(1.0/3.0*m[i]+1.0/6.0*m[i+1])*h[i]);for (i=0;i<n;i++)e[i]=(m[i+1]-m[i])/(6*h[i]);for (i=0;i<n;i++){cout<<"X属于区间["<<x1[i]<<","<<x1[i+1]<<"]时"<<endl<<endl;cout<<"S(x)=";printY(f1[i]);if (c[i]!=0){selectprint(c[i]);cout<<"(x";printY(-x1[i]);cout<<")";}if (m[i]!=0){selectprint((float)(m[i]/2.0));for (q=0;q<2;q++){cout<<"(x";printY(-x1[i]);cout<<")";}}if (e[i]!=0){selectprint(e[i]);for (q=0;q<3;q++){cout<<"(x";printY(-x1[i]);cout<<")";}}cout<<endl<<endl;}cout<<"针对本题，计算S(i+0.5)(i=0,1,...,9):"<<endl;for (i=0;i<10;i++){if ((i+0.5<=x1[n])&&(i+0.5>=x1[0])){calculation((float)(i+0.5));cout<<"S("<<i+0.5<<")="<<s<<endl;}else cout<<i+0.5<<"超出定义域"<<endl;};cout<<endl;}（2）编制的程序求车门的3次样条插值函数S(x)：x属于区间[0,1]时；S(x)=2.51+0.8(x)-0.0014861(x)(x)-0.00851395(x)(x)(x)x属于区间[1,2]时；S(x)=3.3+0.771486(x-1)-0.027028(x-1)(x-1)-0.00445799(x-1)(x-1)(x-1) x属于区间[2,3]时；S(x)=4.04+0.704056(x-2)-0.0404019(x-2)(x-2)-0.0036543(x-2)(x-2)(x-2)x属于区间[3,4]时；S(x)=4.7+0.612289(x-3)-0.0513648(x-3)(x-3)-0.0409245(x-3)(x-3)(x-3) x属于区间[4,5]时；S(x)=5.22+0.386786(x-4)-0.174138(x-4)(x-4)+0.107352(x-4)(x-4)(x-4) x属于区间[5,6]时；S(x)=5.54+0.360567(x-5)+0.147919(x-5)(x-5)-0.268485(x-5)(x-5)(x-5) x属于区间[6,7]时；S(x)=5.78-0.149051(x-6)-0.657537(x-6)(x-6)+0.426588(x-6)(x-6)(x-6) x属于区间[7,8]时；S(x)=5.4-0.184361(x-7)+0.622227(x-7)(x-7)-0.267865(x-7)(x-7)(x-7) x属于区间[8,9]时；S(x)=5.57+0.256496(x-8)-0.181369(x-8)(x-8)+0.0548728(x-8)(x-8)(x-8) x属于区间[9,10]时；S(x)=5.7+0.058376(x-9)-0.0167508(x-9)(x-9)+0.0583752(x-9)(x-9)(x-9)S(0.5)=2.90856 S(1.5)=3.67843 S (2.5)=4.38147 S(3.5)=4.98819 S(4.5)=5.38328 S(5.5)=5.7237S(6.5)=5.59441 S(7.5)=5.42989 S(8.5)=5.65976 S(9.5)=5.7323第六章21．（上机题）常微分方程初值问题数值解（1）编制RK4方法的通用程序；（2）编制AB4方法的通用程序（由RK4提供初值）；（3）编制AB4-AM4预测校正方法的通用程序（由RK4提供初值）；（4）编制带改进的AB4-AM4预测校正方法的通用程序（由RK4提供初值）;（5）对于初值问题h=,应用（1）~（4）中的四种方法进行计算，并将计算结果和精确解取步长0.13=+作比较；y x x()3/(1)（6）通过本上机题，你能得到哪些结论？解：#include<iostream.h>#include<fstream.h>#include<stdlib.h>#include<math.h>ofstream outfile("data.txt");//此处定义函数f(x,y)的表达式//用户可以自己设定所需要求得函数表达式double f1(double x,double y){double f1;f1=(-1)*x*x*y*y;return f1;}//此处定义求函数精确解的函数表达式double f2(double x){double f2;f2=3/(1+x*x*x);return f2;}//此处为精确求函数解的通用程序void accurate(double a,double b,double h){double x[100],accurate[100];x[0]=a;int i=0;outfile<<"输出函数准确值的程序结果:\n";do{x[i]=x[0]+i*h;accurate[i]=f2(x[i]);outfile<<"accurate["<<i<<"]="<<accurate[i]<<'\n';i++;}while(i<(b-a)/h+1);}//此处为经典Runge-Kutta公式的通用程序void RK4(double a,double b,double h,double c) {int i=0;double k1,k2,k3,k4;double x[100],y[100];y[0]=c;x[0]=a;outfile<<"输出经典Runge-Kutta公式的程序结果:\n"; do{x[i]=x[0]+i*h;k1=f1(x[i],y[i]);k2=f1((x[i]+h/2),(y[i]+h*k1/2));k3=f1((x[i]+h/2),(y[i]+h*k2/2));k4=f1((x[i]+h),(y[i]+h*k3));y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;outfile<<"y"<<"["<<i<<"]="<<y[i]<<'\n';i++;}while(i<(b-a)/h+1);}//此处为4阶Adams显式方法的通用程序void AB4(double a,double b,double h,double c) {double x[100],y[100],y1[100];double k1,k2,k3,k4;y[0]=c;x[0]=a;outfile<<"输出4阶Adams显式方法的程序结果:\n";for(int i=0;i<=2;i++){x[i]=x[0]+i*h;k1=f1(x[i],y[i]);k2=f1((x[i]+h/2),(y[i]+h*k1/2));k3=f1((x[i]+h/2),(y[i]+h*k2/2));k4=f1((x[i]+h),(y[i]+h*k3));y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;}int j=3;y1[0]=y[0];y1[1]=y[1];y1[2]=y[2];y1[3]=y[3];do{x[j]=x[0]+j*h;y1[j+1]=y1[j]+(55*f1(x[j],y1[j])-59*f1(x[j-1],y1[j-1])+37*f1(x[j-2],y1[j-2])-9*f1(x[j-3],y1[ j-3]))*h/24;outfile<<"y1"<<"["<<j<<"]="<<y1[j]<<'\n';j++;}while(j<(b-a)/h+1);}//主函数void main(void){double a,b,h,c;cout<<"输入上下区间、步长和初始值:\n";cin>>a>>b>>h>>c;accurate(a,b,h);RK4(a,b,h,c);AB4(a,b,h,c);}float si(int u,int v){float sum=0; int q;for(q=0;q<k;q++)sum+=a[u][q]*a[q][v];sum=a[u][v]-sum;return sum;}void exchange(int g){int t; float temp;for(t=0;t<n;t++){temp=a[k][t];a[k][t]=a[g][t];a[g][t]=temp;}}void analyze(){int t;float si(int u,int v);for(t=k;t<n;t++)a[k][t]=si(k,t);for(t=(k+1);t<m;t++)a[t][k]=(float)(si(t,k)/a[k][k]);}void ret(){int t,z;float sum;x[m-1]=(float)a[m-1][m]/a[m-1][m-1];for(t=(m-2);t>-1;t--){sum=0;for(z=(t+1);z<m;z++)sum+=a[t][z]*x[z];x[t]=(float)(a[t][m]-sum)/a[t][t];}}（5）由经典Runge-Kutta公式得出的结果列在下面的表格中，以及精确值y(x i)和精确值和数值解的误差：由AB4方法得出的结果为：Y1[0]=3 y1[1]=2.997 y1[2]=2.97619 y1[3]=2.92113 y1[4]=2.81839 y1[5]=2.66467 y1[6]=2.4652 y1[7]=2.23308 y1[8]=1.98495 y1[9]=1.73704 y1[10]=1.50219 y1[11]=1.28876 y1[12]=1.10072 y1[13]=0.93871 y1[14]=0.801135y1[15]=0.685335（6）本次上机作业让我知道了在遇到复杂问题中，无法给出解析式的情况下，要学会灵活使用各种微分数值解法，并且能计算出不同方法的精度大小。

数值分析上机题第一章（17题）(1)从2依次累加到N的程序function sn = sum1( n )sn=0;sn=single(sn);for i=2:nai=1/(i^2-1);sn=sn+ai;endend（2）从N依次累加到2的程序function sn = sum2( n )sn=0;sn=single(sn);for i=n:-1:2ai=1/(i^2-1);sn=sn+ai;endend（3）编制求精确值的求和函数sum0function sn = sum0( n )sn=0;sn=single(sn);sn=1/2*(3/2-1/n-1/(n+1));end按第一种顺序得到的值及有效位数如下：N=100时sn0=sum0(100);sn=sum1(100)n=fix(-log10(2*abs(sn-sn0)))得到：sn =0.7400495 n =7N=10e4时sn0=sum0(10e4);sn=sum1(10e4)n=fix(-log10(2*abs(sn-sn0)))得到：sn =0.7498521 n =3N=10e6时sn0=sum0(10e6);sn=sum1(10e6)n=fix(-log10(2*abs(sn-sn0)))得到：sn =0.7498521 n =3按第二种顺序得到的值及有效位数如下：N=100时sn0=sum0(100);sn=sum2(100)n=fix(-log10(2*abs(sn-sn0)))得到：sn =0.7400495 n =7N=10e4时sn0=sum0(10e4);sn=sum2(10e4)n=fix(-log10(2*abs(sn-sn0)))得到：sn =0.7499900 n =7N=10e6时sn0=sum0(10e6);sn=sum2(10e6)n=fix(-log10(2*abs(sn-sn0)))得到：sn =0.7499999 n =7（4）通过这道上机题，我明白了，应用计算机进行求和运算时，求和的顺序不同对结果的精度是有影响的。

第六章 常微分方程数值解法——RK 4法、AB 4法******（学号） *****（姓名）上机题目要求见教材P307,23题。

一、算法原理题目要求采用RK 4法和AB 4法求解最简单的常微分方程初值问题(,),()y f x y a x by a η'=≤≤⎧⎨=⎩ (1)为求解式(1)，采用离散化方法，就是寻求解)(x y 在区间],[b a 上的一系列点<<<<<n x x x x 321上的近似值 ,,,,21n y y y 。

记1(1,2,)i i i h x x i -=-=表示相邻两个节点的间距，称为步长。

求微分方程数值解的主要问题：(1) 如何将微分方程(,)y f x y '=离散化,并建立求其数值解的递推公式; (2) 递推公式的局部截断误差、数值数n y 与精确解)(n x y 的误差估计； (3) 递推公式的稳定性与收敛性. a) Runge-Kutta 方法基本思想：通过在1[,]i i x x +多预报几个点求斜率，并将其加权平均作为k *的近似值，以此构造更高精度的计算公式。

如果每步计算四次函数 的值，完全类似的，可以导出局部截断误差为)(5h O 的四阶Runge-Kutta 公式(RK 4)：1123412132431(22),6(,),(,),221(,),22(,).n n n n n n n n n n y y k k k k k f x y h h k f x y k h k f x h y k k f x h y hk +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ (2)b) Adams 显式公式Runge-Kutta 方法是单步法，计算1+n y 时，只用到n y ， 而已知信息1-n y 、2-n y 等没有被直接利用。

可以设想如果充分利用已知信息1-n y ，2-n y ，…来计算1+n y ，那么不但有可能提高精度，而且大大减少了计算量，这就是构造所谓线性多步法的基本思想。

第六章 常微分方程数值解法——RK 4法、AB 4法******（学号） *****（姓名）上机题目要求见教材P307,23题。

一、算法原理题目要求采用RK 4法和AB 4法求解最简单的常微分方程初值问题(,),()y f x y a x by a η'=≤≤⎧⎨=⎩ (1)为求解式(1)，采用离散化方法，就是寻求解)(x y 在区间],[b a 上的一系列点<<<<<n x x x x 321上的近似值 ,,,,21n y y y 。

记1(1,2,)i i i h x x i -=-=表示相邻两个节点的间距，称为步长。

求微分方程数值解的主要问题：(1) 如何将微分方程(,)y f x y '=离散化,并建立求其数值解的递推公式; (2) 递推公式的局部截断误差、数值数n y 与精确解)(n x y 的误差估计； (3) 递推公式的稳定性与收敛性. a) Runge-Kutta 方法基本思想：通过在1[,]i i x x +多预报几个点求斜率，并将其加权平均作为k *的近似值，以此构造更高精度的计算公式。

如果每步计算四次函数 的值，完全类似的，可以导出局部截断误差为)(5h O 的四阶Runge-Kutta 公式(RK 4)：1123412132431(22),6(,),(,),221(,),22(,).n n n n n n n n n n y y k k k k k f x y h h k f x y k h k f x h y k k f x h y hk +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ (2)b) Adams 显式公式Runge-Kutta 方法是单步法，计算1+n y 时，只用到n y ， 而已知信息1-n y 、2-n y 等没有被直接利用。

可以设想如果充分利用已知信息1-n y ，2-n y ，…来计算1+n y ，那么不但有可能提高精度，而且大大减少了计算量，这就是构造所谓线性多步法的基本思想。