(完整版)数值分析第一次作业
数值分析
数值分析第一次作业信计2 20121314044 王峥虹一、实验内容：1、已知函数在下列各点的值为：38.064.081.092.098.0|0.18.06.04.02.0|y x -------------------试用4次牛顿插值多项式)(4x P 及三次样条函数)(x S （自然边界条件）对数据进行插值，用图给出(){}10,11,1,008.02.0,=+=i i x y x i i i ，，，)(4x P 及)(x S 。
分析：先求4次插值多项式:根据差分形式的牛顿差值公式：))...(](,...,,[...))(](,,[)](,[)()(1010102100100---++--+-+=n n n x x x x x x x f x x x x x x x f x x x x f x f x Px=[0.2,0.4,0.6,0.8,1.0];y=[0.98,0.92,0.81,0.64,0.38];n=length(y);z=zeros(n,n);for i=1:nz(i,1)=y(i);endfor k=2:nfor l=k:nz(l,k)=(z(l,k-1)-z(l-1,k-1))/(x(l)-x(l-k+1));endendz结果：4次牛顿插值多项式为：)6.0)(4.0)(2.0(2083.0)4.0)(2.0(625.0)2.0(3.098.04---------=x x x x x x P )8.0)(6.0)(4.0)(2.0(5208.0-----x x x x再求三次样条插值函数：由上面及已知的：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡---=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡075.65.475.30200005.025.00005.025.00005.025.00000243210M M M M M 程序如下：A=[2,0,0,0,0;0.5,2,0.5,0,0;0,0.5,2,0.5,0;0,0,0.5,2,0.5;0,0,0,0,2];B=[0,-3.75,-4.5,-6.75,0]';M=inv(A)*B结果：则由表达式：j j j j j j j j j j j j j j j j h x x h M y h x x h M y h x x M h x x M x S -⎪⎪⎭⎫ ⎝⎛-+-⎪⎪⎭⎫ ⎝⎛-+-+-=+++++666)(6)()(2111231311,...,1,0-=n j得，三次样条插值多项式为：⎪⎪⎩⎪⎪⎨⎧∈-+-+----∈-+-+----∈-+-+----∈-+-+---=]0.1,8.0[),8.0(9.1)0.1(3036.3)8.0(0)0.1(5893.2]8.0,6.0[),6.0(3036.3)8.0(0857.4)6.0(5893.2)8.0(8929.0]6.0,4.0[),4.0(0857.4)6.0(6536.4)4.0(8929.0)6.0(3393.1]4.0,2.0[),2.0(6536.4)4.0(9.4)2.0(3393.1)4.0(0)(3333333x x x x x x x x x x x x x x x x x x x x x S 绘制4次插值多项式及三次样条插值多项式的图像：代码：x=[0.2,0.4,0.6,0.8,1.0];y=[0.98,0.92,0.81,0.64,0.38];plot(x,y)hold onfor i=1:1:5y(i)=0.98-0.3*(x(i)-0.2)-0.625*(x(i)-0.2)*(x(i)-0.4)-0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0.5 2083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)endk=[0 1 10 11];x0=0.2+0.08*k;y0=zeros(4);for i=1:1:4y0(i)=0.98-0.3*(x(i)-0.2)-0.625*(x(i)-0.2)*(x(i)-0.4)-0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0. 52083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)endplot(x0,y0,'o',x0,y0)hold ony1=spline(x,y,x0)plot(x0,y1,'o')hold ons=csape(x,y,'variational')fnplt(x,'r')hold ongtext('原图像')gtext('三次样条自然边界')gtext('4次牛顿插值')一、实验内容：2、在区间[]11，-上分别取20,10=n 用两组等距节点对龙格函数22511)(x x f +=作多项式插值及三次样条插值，对每个n 值，分别画出插值函数及)(x f 的图形。
数值分析大作业一
数值分析大作业一一、算法设计方案1、求λ1和λ501的值：思路：采用幂法求出按模最大特征值λmax，该值必为λ1或λ501，若λmax小于0，则λmax=λ1；否则λmax=λ501。
再经过原点平移，使用幂法迭代出矩阵A-λmax I的特征值，此时求出的按模最大特征值即为λ1和λ501的另一个值。
2、求λs的值：采用反幂法求出按模最小的特征值λmin即为λs，其中的方程组采用LU分解法进行求解。
3、求与μk最接近的特征值：对矩阵A采用带原点平移的反幂法求解最小特征值，其中平移量为：μk。
4、A的条件数cond(A)=| λmax/λmin|；5、A的行列式的值：先将A进行LU分解，再求U矩阵对角元素的乘积即为A 行列式的值。
二、源程序#include<iostream>#include<iomanip>#include<math.h>#define N 501#define E 1.0e-12 //定义精度常量#define r 2#define s 2using namespace std;double a[N];double cc[5][N];void init();double mifa();double fmifa();int max(int aa,int bb);int min(int aa,int bb);int max_3(int aa,int bb,int cc);void LU();void main(){double a1,a2,d1,d501=0,ds,det=1,miu[39],lamta,cond;int i,k;init();/*************求λ1和λ501********************/a1=mifa();if(a1<0)d1=a1; //若小于0则表示λ1的值elsed501=a1; //若大于0则表示λ501的值for(i=0;i<N;i++)a[i]=a[i]-a1;a2=mifa()+a1;if(a2<0)d1=a2; //若小于0则表示λ1的值elsed501=a2; //若大于0则表示λ501的值cout<<"λ1="<<setiosflags(ios::scientific)<<setprecision(12)<<d1<<"\t";cout<<"λ501="<<setiosflags(ios::scientific)<<setprecision(12)<<d501<<endl;/**************求λs*****************/init();ds=fmifa();cout<<"λs="<<setiosflags(ios::scientific)<<setprecision(12)<<ds<<endl;/**************求与μk最接近的特征值λik**************/cout<<"与μk最接近的特征值λik:"<<endl;for(k=0;k<39;k++){miu[k]=d1+(k+1)*(d501-d1)/40;init();for(i=0;i<N;i++)a[i]=a[i]-miu[k];lamta=fmifa()+miu[k];cout<<"λi"<<k+1<<"\t\t"<<setiosflags(ios::scientific)<<setprecision(12)<<lamta<<en dl;}/**************求A的条件数**************/cout<<"矩阵A的条件式";cond=abs(max(abs(d1),abs(d501))/ds);cout<<"cond="<<setiosflags(ios::scientific)<<setprecision(12)<<cond<<endl;/**************求A的行列式**************/cout<<"矩阵A的行列式";init();LU();for(i=0;i<N;i++){det*=cc[2][i];}cout<<"det="<<setiosflags(ios::scientific)<<setprecision(12)<<det<<endl;system("pause");}/**************初始化函数，给a[N]赋值*************/void init(){int i;for(i=1;i<=501;i++)a[i-1]=(1.64-0.024*i)*sin((double)(0.2*i))-0.64*exp((double)(0.1/i)); }/**************幂法求最大绝对特征值**************/double mifa(){int i,k=0;double u[N],y[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++) //控制最大迭代次数为2000{/***求y(k-1)***/double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;}/****求新的uk****/u[0]=a[0]*y[0]+b*y[1]+c*y[2];u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3]; //前两列和最后两列单独拿出来求中D间的循环求for(i=2;i<N-2;i++){u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];}u[N-2]=c*y[N-4]+b*y[N-3]+a[N-2]*y[N-2]+b*y[N-1];u[N-1]=c*y[N-3]+b*y[N-2]+a[N-1]*y[N-1];/***求beta***/double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}//cout<<"Beta"<<k<<"="<<Beta<<"\t"; 输出每次迭代的beta /***求误差***/error=abs(Beta-Beta_)/abs(Beta);if(error<=E) //若迭代误差在精度水平内则可以停止迭代{return Beta;} //控制显示位数Beta_=Beta; //第个eta的值都要保存下来，为了与后个值进行误差计算 }if(k==2000){cout<<"error"<<endl;return 0;} //若在最大迭代次数范围内都不能满足精度要求说明不收敛}/**************反幂法求最小绝对特¬征值**************/double fmifa(){int i,k,t;double u[N],y[N]={0},yy[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++){double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;yy[i]=y[i]; //用重新赋值，避免求解方程组的时候改变y的值}/****LU分解法解方程组Au=y，求新的***/LU();for(i=2;i<=N;i++){double temp_b=0;for(t=max(1,i-r);t<=i-1;t++)temp_b+=cc[i-t+s][t-1]*yy[t-1];yy[i-1]=yy[i-1]-temp_b;}u[N-1]=yy[N-1]/cc[s][N-1];for(i=N-1;i>=1;i--){double temp_u=0;for(t=i+1;t<=min(i+s,N);t++)temp_u+=cc[i-t+s][t-1]*u[t-1];u[i-1]=(yy[i-1]-temp_u)/cc[s][i-1];}double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}error=abs(Beta-Beta_)/abs(Beta);if(error<=E){return (1/Beta);}Beta_=Beta;}if(k==2000){cout<<"error"<<endl;return 0;} }/**************求两数最大值的子程序**************/int max(int aa,int bb){return(aa>bb?aa:bb);}/**************求两数最小值的子程序**************/int min(int aa,int bb){return(aa<bb?aa:bb);}/**************求三数最大值的子程序**************/int max_3(int aa,int bb,int cc){ int tt;if(aa>bb)tt=aa;else tt=bb;if(tt<cc) tt=cc;return(tt);}/**************LU分解**************/void LU(){int i,j,k,t;double b=0.16,c=-0.064;/**赋值压缩后矩阵cc[5][501]**/for(i=2;i<N;i++)cc[0][i]=c;for(i=1;i<N;i++)cc[1][i]=b;for(i=0;i<N;i++)cc[2][i]=a[i];for(i=0;i<N-1;i++)cc[3][i]=b;for(i=0;i<N-2;i++)cc[4][i]=c;for(k=1;k<=N;k++){for(j=k;j<=min(k+s,N);j++){double temp=0;for(t=max_3(1,k-r,j-s);t<=k-1;t++)temp+=cc[k-t+s][t-1]*cc[t-j+s][j-1];cc[k-j+s][j-1]=cc[k-j+s][j-1]-temp;}//if(k<500){for(i=k+1;i<=min(k+r,N);i++){double temp2=0;for(t=max_3(1,i-r,k-s);t<=k-1;t++)temp2+=cc[i-t+s][t-1]*cc[t-k+s][k-1];cc[i-k+s][k-1]=(cc[i-k+s][k-1]-temp2)/cc[s][k-1];}}}}三、程序结果。
数值分析课第一次作业答案answer1
计算机习题： 1. 作多项式 p，以 −1，0，1 为零点，首项系数为 2，并计算 p(3)。 4
答案：p = poly ([−1, 0, 1])，s = polyval(p, 3)。 2. 已知函数在下列各点的值为 xi 0.2 0.4 0.6 0.8 1.0
2
a 6 6e+154 0 1 1
b 10 10e+154 1 -1e+5 -4
c -4 -4e+154 1 1 3.999999
-1e+155 -7e+155 1e+155 答案：第二种方法更准确，因为第一种方法是一个累加的过程。 matlab 的 x = a : h : b 和 x = a + (0 : n) ∗ h 是第二种方法实现的。 代码： format long e a = 0; b = 8; n = 9; h = (b-a)/n; x(1) = a; y(1) = a; for j = 1:n, x(j+1) = x(j) + h; y(j+1) = y(1) + j*h; end [x',y',(a:h:b)',a+(0:n)’*h] 第二章 插值法 1. 当 x = 1, −1, 2 时，f (x) = 0, −3, 4，求 f (x) 的二次插值多项式。 （计算两遍，分别用拉格朗日插值和牛顿插值）
5
f (xi ) 0.98 0.92 0.81 0.64 0.38 求 4 次牛顿插值多项式 P4 (x) 并画图。 答案： 代码： x=0.2:0.2:1.0; y=[0.98,0.92,0.81,0.64,0.38]; n = length(y); if length(x)~=n, error('x and y are not compatible'); end D = zeros(n,n); D(:,1)=y(:); for j=2:n for i=j:n D(i,j) = (D(i,j-1)-D(i-1,j-1))/(x(i)-x(i-j+1)); end end p=D(1,1)*[zeros(1,n-1),1]; for k=2:n p=p+D(k,k)*[zeros(1,n-k),poly(x(1:k-1))]; end x=0.2:0.01:1.0; z=polyval(p,x); plot(x,z) 比较：p = polyf it(x, y, 4)。
JZX高等数值分析第一次实验作业
相对残差 6.1302e-16 8.7797e-09 8.0295e-09 8.5677e-09 9.1433e-09
a、 m=1 (左为相对残差，右为取对数情况)
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1
0 1
1.5
0
-5
-10
-15
-20
-25
-30
-35
-40
1000
4.761361
2.9675e-11
5.7069e-04
3.5336e+11
1500
15.646571
1.0778e-11
6.8236e-04
2.3672e+10
3000
132.198513
6.3164e-13
1.1700e-04
2.8110e+11
a、n=1000 时步数与相对残差关系图：（上为相对残差，下为取对数情况下结果）
（2）当 A 最大特征值远大于第二个特征值，最小特征值远小于第二个最小特征值时收敛
性情况。
思路：构造题目要求的矩阵 A。首先随机生成 n 阶矩阵 B，B 不满秩，构造对角阵 A1（最
大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值），则由此构
造出对称正定矩阵 A： b1=B’*B; A=b1’*A1*b1。同样设定精确解 Xj 为元素全部为 1 的 n
5、 构造对称不定的矩阵，验证 Lanczos 方法的近似中断，观察收敛曲线中的峰点个数和特
征值的分布关系；观察当出现峰点时，MINRES 方法的收敛性态怎样。
解：思路：类似前两题，首先构造出一个 n 阶对角阵 D，其对角线上有 m 个负值，再对随
数值分析第一次作业
《数值分析》计算作业院系：航空科学与工程学院学号： SY1005512姓名：王天龙日期： 2010年10月31日计算实习说明书目的：训练运用计算机进行科学与工程计算的能力。
要求：1．独立进行算法设计、程序设计和上机运算，并得出正确的结果。
2．编制程序时全部采用双精度，要求按题目的要求设计输出，并执行打印。
3．只能根据题目给出的信息并且只允许一次计算得出全部结果。
题目：第一题 设有501×501的矩阵123499500501a b c b a b cc b a b c A c b a b c c b a b c ba ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦其中0.1(1.640.024)s i n (0.2)0.64 (125i i a i i e i =--= ，，，；0.16b =；0.064c =-。
矩阵A 的特征值12501()ii λ= ，，，满足 125011501||min ||S i i λλλλλ≤≤<<<= ，试求：1.1λ，501λ和S λ的值。
2.A 的与数5011140k kλλμλ-=+最接近的特征值(1239)ik k λ= ，，，。
3.A 的（谱范数）条件数2()cond A 和行列式det A 。
说明：1．在所有的算法中，凡是要给出精度水平ε的，都取1210ε-=。
2．选择算法时，应使A 的所有零元素都不存储。
3．打印以下内容： （1）算法的设计方案。
（2）全部源程序（要求注明主程序和每个子程序的功能）。
（3）特征值1λ，501λ，S λ和(1239)ik k λ= ，，，以及2()cond A ，det A 的值。
4.采用e 型输出所有计算结果，并至少显示12位有效数字。
一、程序算法的设计算法设计方案如下：二、全部源程序编程软件：Fortran：三、计算结果1.特征值1λ，501λ，S λ1-.107001135582E+02λ=，501 .972463398616E+01λ=，-.555823879237E-02S λ=2. (1239)ik k λ= ，，，如下表所示（ZK 代表ik λ）3.A 的条件数2()cond A 和行列式det A 的值2() .192509100000E+04cond A =，det .277059968428+119A =五、讨论这里选取的初始向量为X(i)=1,X={x1,x2,x3,,,,,,,x501},当初始向量与特征向量较近时，收敛较快，若初始向量与特征向量正交，则求解可能失真。
清华大学高等数值分析 第一次实验作业
10
-10
0
100
200
300
400
500
600
700
800
900
迭代次数
图9
m=100时，Lanczos法求解Ax=b的收敛曲线
高等数值分析实验作业一
10
4
Lanzcos 算法的收敛曲线 (阶数 n=1002）
10
2
10
0
||rk||/||b||
10
-2
10
-4
10
-6
10
-8
10
-10
0
200
迭代次数
图12 m=10时，Minres法求解Ax=b的收敛曲线
10
2
Minres 算法的收敛曲线 (阶数 n=1002）
10
0
10
-2
||rk||/||b||
10
-4
10
-6
10
-8
10
-10
0
100
200
300
400
500
600
700
迭代次数
图13
10
2
m=50时，Minres法求解Ax=b的收敛曲线
10
0
Lanzcos 算法的收敛曲线 (阶数 n=1002)
m=10 m=50 m=100 m=400 m=800
10
-2
10
-4
||rk||/||b||
10
-6
10
-8
10
-10
10
-12
0
2
4
6
8
10
12
14
16
数值分析第一次作业解答
数值分析第一次作业解答1：(a) —个问题的病态性如何，与求解它的算法有关系。
x ;(b) 无论问题是否病态，好的算法都会得到它好的近似解。
x ；(C)计算中使用更高的精度，可以改善问题的病态性。
X ；(d) 用一个稳定的算法计算一个良态问题，一定会得到他好的近似解。
V;(e) 浮点数在整个数轴上是均匀分布。
x ;(f) 浮点数的加法满足结合律。
x(g) 浮点数的加法满足交换律。
X ;(h) 浮点数构成有效集合。
V;(i) 用一个收敛的算法计算一个良态问题，一定得到它好的近似解2: 程序t=0.1;n=1:10;e=n/10-n*te = 1.0e-015 *[ 0 0 -0.0555 0 0-0.1110 -0.1110 0 0 0] 由舍人误差造成n=3，6，7 时的结果不为零。
4:两种等价的一元二次方程求解公式-b - Pb2 - 4acx =2a2cx 二-b b2 - 4ac对a=1, b=-100000000, c=1，应采用哪种算法?A二[1,-100000000,1];roots(A);可得：X1 = 100000000;x2=0a=1;b=-100000000;c=1;x1仁(-b-sqrt(b*b-4*a*c))/(2*a)x12=(-b+sqrt(b*b-4*a*c))/(2*a)x2仁2*c/(-b-sqrt(b*b-4*a*c))x22=2*c/(-b+sqrt(b*b-4*a*c))由第一种算法：X1 = 100000000;x2=7.45058 X10由第二种算法：X1 = 13417728;x2=-1.0 X108原因：太小的数作分母。
5:程序:fun cti on y=tt(x)s=0;t=x;n=1;while s+t~=s;s=s+t;t=-x A2/(( n+1)*( n+2))*tn=n+2;endntt(2n 1)eps)(a)t小于计算机的计算精度。
数值分析作业题(1)
第一章 误差与算法1. 误差分为有__模型误差___, _观测误差___, __方法误差____, ___舍入误差____, Taylor 展开式近似表达函数产生的误差是_方法误差 .2. 插值余项是插值多项式的 方法误差。
0.2499作为1/4的近似值, 有几位有效数字？00.24990.249910,0m =⨯=即，031|0.2499|0.00010.5100.510,34m n n ---=<⨯=⨯=即22 3.1428751...,7=作为圆周率的近似值，误差和误差限分别是多少，有几位有效数字？2133.142875 3.14159260.00126450.5100.510---=<⨯=⨯有3位有效数字.* 有效数字与相对误差的关系3. 利用递推公式计算积分110,1,2,...,9n x n I x e dx n -==⎰错误!未找到引用源。
, 建立稳定的数值算法。
该算法是不稳定的。
因为:11()()...(1)!()n n n I n I n I εεε-=-==-111n n I I n n -=-， 10110I =4. 衡量算法优劣的指标有__时间复杂度,__空间复杂度_.时间复杂度是指: , 两个n 阶矩阵相乘的乘法次数是 , 则称两个n 阶矩阵相乘这一问题的时间复杂度为 .二 代数插值1.根据下表数据建立不超过二次的Lagrange 和Newton 插值多项式, 并写出误差估计式, 以及验证插值多项式的唯一性。
x 0 1 4f(x) 1 9 3Lagrange:设0120120,1,4;()1()9()3x x x f x f x f x ======则，， 对应 的标准基函数 为:1200102()()(1)(x 4)1()(1)(x 4)()()(01)(04)4x x x x x l x x x x x x ----===------ 1()...l x =2()...l x =因此, 所求插值多项式为:220()()()....i i i P x f x l x ===∑ (3)2()()(0)(1)(x 4)3!f R x x x ξ=--- Newton:构造出插商表:xi f(xi ) 一 二 三0 11 9 84 3 -2 -5/2所以, 所求插值多项式为:2001001201()()[,]()[,,]()()518(0)(0)(1)2...P x f x f x x x x f x x x x x x x x x x =+-+--=+----=插值余项: 2()[0,1,4,](0)(1)(x 4)R x f x x x =---2. 已知函数f(0)=1,f(1)=3,f(2)=7,则f[0,1]=___2________, f[0,1,2]=____1______)('],[000x f x x f =3.过0,1两节点构造三次Hermite 插值多项式, 使得满足插值条件: f(0)=1. .’(0)=... f(1.=2. .’(1)=1设0101010,1,()1()2'()0,'()1x x f x f x f x f x ======则，， 写出插商表:xi f(xi) 一 二 三0 10 1 01 a 1 11 a 1 0 a-1因此, 所求插值多项式为:插值余项:222()[0,0,1,1,](1)R x f x x x =-4.求f(x)=sinx 在[a,b]区间上的分段线性插值多项式, 并写出误差估计式。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
问题1：20.给定数据如下表：试求三次样条插值S(x)，并满足条件 (1)S`(0.25)=1.0000，S`(0.53)=0.6868； （2）S ’’(0.25)=S ’’(0.53)=0。
分析：本问题是已知五个点，由这五个点求一三次样条插值函数。
边界条件有两种，（1）是已知一阶倒数，（2）是已知自然边界条件。
对于第一种边界(已知边界的一阶倒数值)，可写出下面的矩阵方程。
⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡432104321034322110d M M M M M 200020000020022d d d d λμμλμλμλ其中μj =j1-j 1-j h h h +，λi=j1-j j h h h +，dj=6f[x j-1,x j ,x j+1]， μn =1，λ0=1对于第一种边界条件d 0=0h 6（f[x 0,x 1]-f 0`），d n =1-n h 6（f`n-f `[x n-1,x n ]） 解：由matlab 计算得：由此得矩阵形式的线性方程组为：⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡ 2.1150-2.4286-3.2667-4.3143-5.5200-M M M M M 25714.00001204286.000004000.026000.0006429.023571.0001243210解得 M 0=-2.0286;M 1=-1.4627;M 2= -1.0333; M 3= -0.8058; M 4=-0.6546S(x)=⎪⎪⎩⎪⎪⎨⎧∈-+-+-∈-+-+-∈-+-+-∈-+-+-]53.0,45.0[x 5.40x 9.1087x 35.03956.8.450-x 1.3637-x .5301.67881- ]45.0,39.0[x 9.30x 11.188x 54.010.418793.0-x 2.2384-x .450(2.87040-]39.0,30.0[x 03.0x 6.9544x 9.30 6.107503.0-x 1.9136-x .3902.708779-]30.0,25.0[x 5.20x 10.9662x 0.3010.01695.20-x 4.8758-x .3006.76209-33333333），（）（）（）（），（）（）（）），（）（）（）（），（）（）（）（Matlab 程序代码如下：function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。
x=[0.25 0.30 0.39 0.45 0.53];y=[0.5 0.5477 0.6245 0.6708 0.7280];n=5,s=1.0,t=0.6868for j=1:1:n-1h(j)=x(j+1)-x(j);endfor j=2:1:n-1r(j)=h(j)/(h(j)+h(j-1));endfor j=1:1:n-1u(j)=1-r(j);endfor j=1:1:n-1f(j)=(y(j+1)-y(j))/h(j);endfor j=2:1:n-1d(j)=6*(f(j)-f(j-1))/(h(j-1)+h(j));endd(1)=6*(f(1)-s)/h(1)d(n)=6*(t-f(n-1))/h(n-1)a=zeros(n,n);for j=1:1:na(j,j)=2;endr(1)=1;u(n)=1;for j=1:1:n-1a(j+1,j)=u(j+1);a(j,j+1)=r(j);endb=inv(a)m=b*d'p=zeros(n-1,4); %p矩阵为S(x)函数的系数矩阵for j=1:1:n-1p(j,1)=m(j)/(6*h(j));p(j,2)=m(j+1)/(6*h(j));p(j,3)=(y(j)-m(j)*(h(j)^2/6))/h(j);p(j,4)=(y(j+1)-m(j+1)*(h(j)^2/6))/h(j);end对于第二中边界，已知边界二阶倒数，可写出下面的矩阵：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡432104321034322110d M M M M M 2000200002002002d d d d λμμλμλμλ 其中μj =j1-j 1-j h h h +，λi=j1-j j h h h +，dj=6f[x j-1,x j ,x j+1]，μn =λ0=0 d 0=d n =0解：由matlab 计算得：由此得矩阵形式的线性方程组为：⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡0 2.4286-3.2667-4.3143-0M M M M M 25714.00000204286.000004000.026000.0006429.023571.000243210解得M 0=0 ;M 1= -1.8795;M 2= -0.8636; M 3= -1.0292; M 4=0S(x)=⎪⎪⎩⎪⎪⎨⎧∈-+-+--∈-+-+---∈-+-+---∈-+-+--]53.0,45.0[x )45.0(1000.953.03987.845.00-x .5302.1442-]45.0,39.0[x 9.30x 1903.11x 54.0417.1093.0-x 8590.2x .4503990.2]39.0,30.0[x 03.0x 9518.6x 9.301137.603.0-x 5993.1x .3904806.3]30.0,25.0[x )25.0(9697.1030.0 010.25-x 2652.6x 0.30 033333333，）（）（）（），（）（）（）（），（）（）（）（，）（）（）（x x x x x matlab 程序代码如下：function tgsanci(n,s,t) %n 代表元素数， x=[0.25 0.30 0.39 0.45 0.53];y=[0.5 0.5477 0.6245 0.6708 0.7280];n=5for j=1:1:n-1h(j)=x(j+1)-x(j);endfor j=2:1:n-1r(j)=h(j)/(h(j)+h(j-1));endfor j=1:1:n-1u(j)=1-r(j);endfor j=1:1:n-1f(j)=(y(j+1)-y(j))/h(j);endfor j=2:1:n-1d(j)=6*(f(j)-f(j-1))/(h(j-1)+h(j));endd(1)=0d(n)=0a=zeros(n,n);for j=1:1:na(j,j)=2;endr(1)=0;u(n)=0;for j=1:1:n-1a(j+1,j)=u(j+1);a(j,j+1)=r(j);endb=inv(a)k=am=b*d'p=zeros(n-1,4); %p矩阵为S(x)函数的系数矩阵for j=1:1:n-1p(j,1)=m(j)/(6*h(j));p(j,2)=m(j+1)/(6*h(j));p(j,3)=(y(j)-m(j)*(h(j)^2/6))/h(j);p(j,4)=(y(j+1)-m(j+1)*(h(j)^2/6))/h(j);endp由下面的一段程序，画出自然边界与第一边界的图形：程序代码如下：x=[0.25 0.30 0.39 0.45 0.53];y=[0.5 0.5477 0.6245 0.6708 0.7280];s=csape(x,y,'variational')fnplt(s,'r')hold onxlabel('x')ylabel('y')gtext('三次样条自然边界')plot(x,y,'o',x,y,':g')hold ons.coefsr=csape(x,y,'complete',[1.0 0.6868])fnplt(r,'b')gtext('三次样条第一边界')hold on r.coefs图中的三条线几乎重合。
xy图20-1 在一小段区间内的图形xy图20-2 在整个给出的区间的图形问题2：1已知函数在下列各点的值为试用4次牛顿插值多项式P 4（x ）及三次样条函数S （x ）（自然边界条件）对数据进行插值。
用图给出{（x i ，y i ），x i =0.2+0.08i ，i=0，1, 11, 10}，P 4（x ）及S （x ）。
分析：（1）要用4次牛顿插值多项式处理数据， P n =f(x 0)+f[x 0,x 1](x-x 0)+ f[x 0,x 1,x 2](x-x 0) (x-x 1)+···+ f[x 0,x 1，···x n ](x-x 0) ···(x-x n-1) 首先我们要知道牛顿插值多项式的系数，即均差表中得部分均差。
解：由所以有四次插值牛顿多项式为：P 4（x ）=0.98-0.3(x-0.2)-0.62500 (x-0.2)(x-0.4) -0.20833 (x-0.2)(x-0.4)(x-0.6)-0.52083(x-0.2)(x-0.4)(x-0.6)（x-0.8）计算牛顿插值中多项式系数的程序如下： function varargout=newtonliu(varargin) clear,clcx=[0.2 0.4 0.6 0.8 1.0];fx=[0.98 0.92 0.81 0.64 0.38]; newtonchzh(x,fx);function newtonchzh(x,fx) %由此函数可得差分表 n=length(x);fprintf('*****************差分表*****************************\n'); FF=ones(n,n); FF(:,1)=fx'; for i=2:nfor j=i:nFF(j,i)=(FF(j,i-1)-FF(j-1,i-1))/(x(j)-x(j-i+1)); end endfor i=1:nfprintf('%4.2f',x(i)); for j=1:ifprintf('%10.5f',FF(i,j)); endfprintf('\n'); end（2）用三次样条插值函数由上题分析知,要求各点的M 值： 有matlab 编程计算求出个三次样条函数：⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡06.7500-4.5000-3.7500-0M M M M M 2.5000000020.50000000.50002.500000.50002.5000000 243210 解得：M 0=0 ;M 1= -1.6071;M 2= -1.0714; M 3= -3.1071; M 4=0⎪⎪⎩⎪⎪⎨⎧∈-+-+--∈-+-+-∈-+-+--∈-+-+---]0.1,8.0[x )8.0(9.10.1 3036.38.00-x 0.12.5893-]8.0,6.0[x 6.0x 3036.3x 8.00857.4.60-x 5893.2-x .800.8929-]6.0,4.0[x 4.0x 7508.4x 6.04.6536 .40-x 8929.0x .601.3393- ]4.0,2.0[x )2.0(6536.44.0900.42.0 1.33934.00 3333333，）（）（）（），（）（）（）（），（）（）（）（，）（）（）（x x x x x x x 三次样条插值函数计算的程序如下：function tgsanci(n,s,t) %n 代表元素数，s,t 代表端点的一阶导。