当前位置:文档之家› 数值分析第一次作业

数值分析第一次作业

问题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)首先我们要知道牛顿插值多项式的系数,即均差表中得部分均差。

解:由matlab 计算得: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 代表端点的一阶导。

相关主题