数值分析实验题(-华科)

数值分析实验作业

专业： ： 学号：

实验2.1 多项式插值的振荡现象

[问题提出]：

考虑在一个固定的区间上用插值逼近一个函数，显然Lagrange插值中使用的节点越多，插值多项式的次数就越高，我们自然关心插值多项的次数增加时，Ln(x)是否也更加靠近逼近的函数，Runge给出的例子是极著名并富有启发性的，设区间[-1,1]上函数

f(x)[实验内容]：

考虑区间[-1，1]的一个等距离划分，分点为

1

125x2xi1则拉格朗日插值多项式为

2i, i0,1,2,...,n nLn(x)1l(x) 2ii0125xin其中，li(x)，i=0,1,2,…,n是n次Lagrange插值函数。

[实验要求]：

〔1〕选择不断增大的分点数目n=2，3，…画出原函数f(x)及插值多项式函数Ln(x)在[-1，1]上的图像，比较并分析实验结果。

〔2〕选择其他的函数，例如定义在区间[-5，5]上的函数，

h(x)x,g(x)arctanx 41x重复上述的实验看其结果如何。

解：以下的f(x)、h(x)、g(x)的为插值点用“*”表示，朗格朗日拟合曲线用连续曲线表示。通过三个函数的拉格朗日拟合可以看到，随着插值点的增加，产生Rung现象。 〔1〕 f(x)

多项式求值的振荡现象 n=210.90.80.70.60.50.40.30.20.10 -1-0.8-0.6-0.4-0.20x0.20.40.60.81f(x)lagrange(x) 10.90.80.70.60.50.40.30.20.10 -1-0.8-0.6-0.4-0.20x0.20.40.60.81f(x)lagrange(x)多项式求值的振荡现象 n=3 y

y

第 1 页 共 13 页

多项式求值的振荡现象 n=41 f(x)0.8lagrange(x)0.60.4y0.20-0.2-0.4-1 -0.8-0.6-0.4-0.200.20.40.60.81x多项式求值的振荡现象 n=61.2 f(x)1lagrange(x)0.80.6y0.40.20-0.2-1 -0.8-0.6-0.4-0.200.20.40.60.81x多项式求值的振荡现象 n=81 f(x)0.8lagrange(x)0.60.40.2y0-0.2-0.4-0.6-0.8-1-1 -0.8-0.6-0.4-0.200.20.40.60.81x多项式求值的振荡现象 n=101.6 f(x)1.4lagrange(x)1.210.8y0.60.40.20-0.2-0.4-1 -0.8-0.6-0.4-0.200.20.40.60.81x多项式求值的振荡现象 n=51.2 f(x)1lagrange(x)0.80.6y0.40.20-0.2-1 -0.8-0.6-0.4-0.200.20.40.60.81x

多项式求值的振荡现象 n=71 f(x)0.9lagrange(x)0.80.70.6y0.50.40.30.20.10-1 -0.8-0.6-0.4-0.200.20.40.60.81x

多项式求值的振荡现象 n=91 f(x)0.8lagrange(x)0.60.4y0.20-0.2-0.4-1 -0.8-0.6-0.4-0.200.20.40.60.81x

多项式求值的振荡现象 n=111 f(x)0.9lagrange(x)0.80.70.6y0.50.40.30.20.10-1 -0.8-0.6-0.4-0.200.20.40.60.81x

第 2 页 共 13 页

(2) h(x)

多项式求值的振荡现象 n=20.6 h(x)0.4lagrange(x)0.20y-0.2-0.4-0.6-0.8-5 -4-3-2-1012345x多项式求值的振荡现象 n=40.6 h(x)0.4lagrange(x)0.20y-0.2-0.4-0.6-0.8-5 -4-3-2-1012345x多项式求值的振荡现象 n=60.6 h(x)0.4lagrange(x)0.20y-0.2-0.4-0.6-0.8-5 -4-3-2-1012345x多项式求值的振荡现象 n=80.6 h(x)0.4lagrange(x)0.20y-0.2-0.4-0.6-0.8-5 -4-3-2-1012345x多项式求值的振荡现象 n=30.6 h(x)0.4lagrange(x)0.20y-0.2-0.4-0.6-0.8-5 -4-3-2-1012345x

多项式求值的振荡现象 n=50.8 h(x)0.6lagrange(x)0.40.2y0-0.2-0.4-0.6-0.8-5 -4-3-2-1012345x

多项式求值的振荡现象 n=71.5 h(x)lagrange(x)10.5y0-0.5-1-1.5-5 -4-3-2-1012345x

多项式求值的振荡现象 n=93 h(x)lagrange(x)21y0-1-2-3-5 -4-3-2-1012345x

第 3 页 共 13 页

多项式求值的振荡现象 n=101 h(x)0.8lagrange(x)0.60.40.2y0-0.2-0.4-0.6-0.8-1-5 -4-3-2-1012345x

(3) g(x)

多项式求值的振荡现象 n=21.5 g(x)lagrange(x)10.5y0-0.5-1-1.5-5 -4-3-2-1012345x多项式求值的振荡现象 n=41.5 g(x)lagrange(x)10.5y0-0.5-1-1.5-5 -4-3-2-1012345x多项式求值的振荡现象 n=61.5 g(x)lagrange(x)10.5y0-0.5-1-1.5-5 -4-3-2-1012345x多项式求值的振荡现象 n=115 h(x)4lagrange(x)321y0-1-2-3-4-5-5 -4-3-2-1012345x

多项式求值的振荡现象 n=32 g(x)1.5lagrange(x)10.5y0-0.5-1-1.5-2-5 -4-3-2-1012345x

多项式求值的振荡现象 n=51.5 g(x)lagrange(x)10.5y0-0.5-1-1.5-5 -4-3-2-1012345x

多项式求值的振荡现象 n=72 g(x)1.5lagrange(x)10.5y0-0.5-1-1.5-2-5 -4-3-2-1012345x

第 4 页 共 13 页

多项式求值的振荡现象 n=82 g(x)1.5lagrange(x)10.5y0-0.5-1-1.5-2-5 -4-3-2-1012345x多项式求值的振荡现象 n=101.5 g(x)lagrange(x)10.5y0-0.5-1-1.5-5 -4-3-2-1012345x多项式求值的振荡现象 n=91.5 g(x)lagrange(x)10.5y0-0.5-1-1.5-5 -4-3-2-1012345x

多项式求值的振荡现象 n=112.5 g(x)2lagrange(x)1.510.5y0-0.5-1-1.5-2-2.5-5 -4-3-2-1012345x

第 5 页 共 13 页

实验3.1 最小二乘法拟合

编制以函数{x}k0为基的多项式最小二乘拟合程序，并用于对表中的数据作三次多项式最小二乘拟合。 knxi yi -1.0 -4.447 -0.5 -0.452 *n0.0 0.551 0.5 0.048 1.0 -0.447 21.5 0.549 2.0 4.552 取权数i1,求拟合曲线xk0*kk中的参数{k}，平方误差，并作离散数据{xi,yi}的拟合函数y*(x)的图形。 解：三次多项式的拟合曲线为： y(x)a0a1xa2x2a3x3 此题中权函数(x)1，即W=(1,1,1,1,1,1,1) 利用法方程AAa=AY 求解这个方程组，就可以得到系数a。

解之得：00.54912,13.9683105,22.9977,31.9991 故拟合的函数为:y0.549123.9683105x2.9977x21.9991x3, 平方误差为：2.176191667187105e-05

拟合的函数图像如下：

TT3次多项式拟合,平方误差=2.1762e-05543210-1-2-3-4-5 -1-0.500.5x11.52y 离散值拟合曲线第 6 页 共 13 页

实验5.1 常微分方程性态和R-K法稳定性试验

[试验目的]：

考察下面的微分方程右端项中函数y前面的参数对方程性态的影响〔它可使方程为好条件的或坏条件的〕和研究计算步长对R-K法计算稳定性的影响。 [实验题目]：

常微分方程初值问题

yyx1,0x1 y(0)1,其中，5050。其精确解为y(x)exx

[实验要求]：

〔1〕对于参数，分别去四个不同的数值：一个大的正值，一个小的正值，一个绝对值小的负值和一个绝对值大的负值。取步长h0.01，分别用经典R-K法计算，将四组计算结果画在同一张图上，进行比较并说明相应初值问题的性态。

〔2〕对于参数为一个绝对值不大的负值和两个计算步，一个计算步使参数h在经典R-K法的稳定域内，另一个步长在经典的R-K法的稳定域外。分别用经典R-K法计算并比较计算结果。取全域等距的10个点上的计算值，列表说明。

解:对于4阶R-K法 绝对稳定区为：2.785h0 这里，所以绝对稳定区为：2.785h0 〔1〕

对于h0.01，绝对稳定区：278.50 a h 98765y2 0.01 1 0.01 -1 0.01 -2 0.01 精确解数值解微分方程数值解43210 00.10.20.30.40.5x0.60.70.80.91

第 7 页 共 13 页

(2)

对于20，稳定区0h0.1391 a h

微分方程数值解,a=-20,h=0.011.110.960.80.70.60.50.420.30.20.1 00.10.20.30.40.5x0.60.70.80.9110 0yy-20 0.01 -20 0.15 微分方程数值解,a=-20,h=0.15 精确解数值解87精确解数值解 5430.10.20.30.4x0.50.60.70.80.9

x 数值解y1 y〔精确解〕 〔a=-20,h=0.01〕 y1-y 数值解y2 〔a=-20,h=0.15〕 y1-y 0.15 0.30 0.45 0.60 0.75 0.199787 0.302479 0.450123 0.600006 0.750000 0.1997 2.35E-06 0.302479 2.34E-07 0.450123 1.75E-08 0.600006 1.16E-09 0.750000 7.23E-11 0.900000 4.32E-12 h=0.15时，数值解不稳定。

1.525000 2.190625 3.049609 4.174463 5.6886 7.657969 1.325213 1.888146 2.599486 3.574457 4.914886 6.757969 0.90 0.900000 可见h=0.01时，数值解稳定

第 8 页 共 13 页

程序源代码

function testCharpt2_1

%对数值分析实验题第2章第1题进行分析

promps={'输入f为选择f(x)；输入h为选择h(x)；输入g为选择g(x)'}; result=inputdlg(promps,'请选择实验函数'); chooseFunction=char(result); switch chooseFunction case 'f'

f=inline('1./(1+25*x.^2)'); a=-1; b=1;

nameFuc='f(x)'; case 'h'

f=inline('x./(1+x.^4)'); a=-5; b=5

nameFuc='h(x)' case 'g'

f=inline('atan(x)'); a=-5; b=5

nameFuc='g(x)' end

% promps2={'n='};

% nNumble=inputdlg(promps2,'请输入分点数n'); nNumble=[2:11]

for i=1:length(nNumble)

x=linspace(a,b,nNumble(i)+1); y=feval(f,x); xx=a:0.1:b;

yy=lagrange(x,y,xx) figure

fplot(f,[a,b],'*') hold on

plot(xx,yy,'LineWidth',2) xlabel('x') ylabel('y')

legend(nameFuc,'lagrange(x)')

nameTitle=['多项式求值的振荡现象',' n=',num2str(nNumble(i))] title(nameTitle,'FontSize',14); grid on end

第 9 页 共 13 页

function yy=lagrange(x,y,xx) %s实现拉格朗日插值

%输入参数x，y分别为已知插值点的自变量和因变量 %输入参数xx为拟合点的自变量值

%输出参数yy为对应自变量xx的拟合值 xLength=length(x); xxLength=length(xx); for i1=1:xxLength yy(i1)=0;

for i2=1:xLength p=1;

for i3=1:xLength if(i2~=i3)

p=p*(xx(i1)-x(i3))/(x(i2)-x(i3)); end end

yy(i1)=yy(i1)+p*y(i2); end end

第 10 页 共 13 页

function testCharpt3_1()

%对数值分析实验题第3章第1题进行分析 %输入参数：自变量x，因变量y %输入参数：多项式拟合次数n clc clear

format long

x=[-1.0,-0.5,0.0,0.5,1.0,1.5,2.0]

y=[-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552] n=3 A=[];

for i=1:length(x)

A=[A;1 x(i) x(i)^2 x(i)^3] end

A2=A'*A;

a=inv(A2)*A'*y'%多项式的系数 % a=roundn(a,-6)

yy=a(1)+a(2)*x+a(3)*x.^2+a(4)*x.^3; r=(y-yy)*(y-yy)' %平方误差 clf

hold on plot(x,y,'or'); x2=-1:0.01:2;

y2=a(1)+a(2)*x2+a(3)*x2.^2+a(4)*x2.^3; plot(x2,y2,'LineWidth',2); legend('离散值','拟合曲线') xlabel('x'); ylabel('y');

title(['3次多项式拟合,平方误差=',num2str(r)],'FontSize',14); grid

on

第 11 页 共 13 页

function testCharpt5_1

%对数值分析实验题第3章第1题进行分析 %输入参数：参数a，步长h %精确解和数值解图形比照 %%第1问输入

a=[2 1 -1 -2]% 输入a的取值

h=[0.01 0.01 0.01 0.01]%输入h的取值 %%第2问输入

% a=[-20 -20]% 输入a的取值 % h=[0.01 0.15]%输入h的取值 %%

func=inline('1+(y-x).*a');%定义函数 for i=1:length(a)

x=0:h(i):1;%求解区间 y=x;

N=length(x); y(1)=1; for n=1:N-1

k1=func(a(i),x(n),y(n));

k2=func(a(i),x(n)+h(i)/2,y(n)+k1*h(i)/2); k3=func(a(i),x(n)+h(i)/2,y(n)+k2*h(i)/2); k4=func(a(i),x(n)+h(i),y(n)+k3*h(i)) ;

y(n+1)=y(n)+h(i)*(k1+2*k2+2*k3+k4)/6;%数值解 end

y0=exp(a(i)*x)+x;%精确解

% figure()%如果叠绘图去掉此句命令 plot(x,y0)

hold on plot(x,y,'*')

legend('精确解','数值解') xlabel('x'); ylabel('y');

title(['微分方程数值解,a=',num2str(a(i)),',h=',num2str(h(i))],'FontSize',14); grid on end

第 12 页 共 13 页