XuQiang
上一次上机课我们做了P71例6这个问题,但是做出来的同学很少很少。好多同学都在机房忙着看基础知识,部分同学按照课本P70的程序提示来做,这里我 需要谈一下我的看法。 首先,我认为大家应该在课下提前做一些准备工作。我们授课课堂上布置了问题之后,同学们如果能在课下提前想一下程序的大体步骤或者简单写一些伪代码,那 么到了上机的时候效率必然比现在高。 其次,关于课本上的程序提示,我认为应该在自己编写出程序之后再去读。我认为写程序处理一个问题,应先用自己的想法把它编写调试实现出来(第一次实现时 程序一般不是最优的),然后再不断的改进优化(节省计算量+节省存储空间),直到满意为止。而课本上P70的程序提示中,很多地方用了我们一眼看不懂的 变量名字,它可能为了节省存储量,一会让该变量代替这个东西,一会让该变量代替另一个东西。所以说,只有我们把自己的程序实现出来,再去读课本提示时才 能看到它处理问题时的一些想法和妙用。读别人优秀的代码是一件好事情,甚至读伪代码都是好事情,可是如果不对比着自己的不足盲目的去读,只会造成比着葫 芦画葫芦,并不起很大的作用。 以上个人意见,仅供参考。 最后,附件里面贴出我给的两个Matlab程序,供大家参考。
spline_test.m内有较为详细的注释。程序没有优化到极致,尽量保持了一种阅读性较强的风格。
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% spline_test.m
% 本程序用于生成自然边界条件(第2种边界条件)下的样条曲线,并估计函数在某些点处的近似函数值
% 其中x,y为同样长度的向量,代表离散数据。我们要在给定x,y这些数据的前提下,对其进行三次样条插值
% 程序大体分为以下几个步骤:
% 1. 声明(定义)所需的系数矩阵A,右端向量d以及需要的常向量mu,lambda,h
% 2. 为系数矩阵和右端向量赋值
% 3. 求解线性方程组
% 4. 利用线性方程组求解得到的结果,计算样条表达式中的系数C,详见P65 (8.16)
% 5. 绘制样条函数图象
% 6. 计算给定点x_nu处的样条函数值,即被逼近的函数在该点处的近似值
%
% 2008.11.23 已调试,未优化
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% 清空内存,命令行,关闭所有的figure窗口
clear; clc; close all;
% 给定数据x,y,该处x,y可以在满足条件1时任意给定,只需它们是同样长度的向量即可,但是注意最后计算样条函数在
% x_nu处的函数值的时候,x_nu必需位于我们考虑的区间[x(1),x(n)]上。
% 条件1:x的每个分量按下标的顺序依次增大
x = [1 2 4 5 6 7];
y = [1 3 4 2 3 4];
% 取得x或者y的长度
n = length(x);
% 声明A,d,mu,lambda,h,并对它们所有的元素赋予初值0
% 关于它们的大小请参考课本P66 最上面的式子
A = zeros(n-2,n-2);
d = zeros(n-2,1);
mu = zeros(n-2,1);
lambda = zeros(n-2,1);
h = zeros(n-1,1);
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% 先给出第一个步长h(1),因为在下面的循环中会被用到
h(1) = x(2) - x(1);
% for循环给h,lambda,mu和右端向量d赋值
for j = 2:n-1
h(j) = x(j+1) - x(j);
lambda(j-1) = h(j)/(h(j) + h(j-1));
mu(j-1) = 1 - lambda(j-1);
d(j-1) = 6 * ((y(j+1) - y(j))/h(j) - (y(j) - y(j-1))/h(j-1))/(h(j) + h(j-1));
end
% 另一个for循环给A赋值,注意:A为三对角阵,所以赋值很容易实现,当然也可以用另外的方法。
for j = 1:n-3
A(j,j) = 2;
A(j,j+1) = lambda(j);
A(j+1,j) = mu(j+1);
end
A(n-2,n-2) = 2;
% ++++++++++++++++++++关于前面两个for循环的一点注释+++++++++++++++++++++++++
% 1.两个for循环没有合并在一起的原因在于,分开循环从语句来看比较明朗,如果合并在一起来赋值的话
% 可能需要if判断。两种赋值方式孰优孰劣,未比较过。我想这不是三次样条插值最关键的地方。
% 2.三对角矩阵A是稀疏矩阵,是否可以采用存储量更小的方式存储它?
% 当然可以,本程序中利用的存储方式,存储量大,一般三对角矩阵利用带状存储比较合适
% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% 求解以A为系数矩阵,d为右端向量的线性方程组AM=d
% 这里直接用A的逆乘d来得到结果,三对角方程的追赶法实现,还是留到后面让同学们自己来做吧
% 这里的M中每个分量仅代表第2个节点至倒数第2个节点处的二阶导数
M = A\d;
% 把第1个和最后1个节点处的二阶导数‘0’,即已知的自然边界条件放入M的首尾
% 生成我们下面要用到的M,现在的M的每个分量代表的就是每个节点处的二阶导数了
M = [0;M;0];
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% 生成样条表达式中的系数C,详见P65 (8.16)
C = zeros(n-1,3);
for j = 1:n-1
C(j,1) = (y(j+1) - y(j))/h(j) - (2*M(j) + M(j+1))*h(j)/6;
C(j,2) = M(j)/2;
C(j,3) = (M(j+1) - M(j))/(6*h(j));
end
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% 绘制样条函数图象部分
% 给定一个小步长h_test,一般0.01就可以
h_test = 0.01;
% x_test是是一个向量
% 它以x(1)与x(n)为第1个分量和最后一个分量,h_test为中间每个相邻分量之间的距离
x_test = x(1):h_test:x(n);
% 取得x_test的长度
m = length(x_test);
% 下面计算与x_test的每个分量处相应的样条函数值,存入y_test
y_test = zeros(m,1);
j = 1;
for i = 1:m
if x_test(i) <= x(j+1) % 判断x_test(i)位于哪个单元[x(j),x(j+1)]中
y_test(i) = y(j) + C(j,1)*(x_test(i) - x(j)) + C(j,2)*(x_test(i) - x(j))^2 + C(j,3)*(x_test(i) - x(j))^3;
else
j = j + 1;
y_test(i) = y(j) + C(j,1)*(x_test(i) - x(j)) + C(j,2)*(x_test(i) - x(j))^2 + C(j,3)*(x_test(i) - x(j))^3;
end
end
% 绘制插值节点和样条函数图象
plot(x,y,'*',x_test,y_test,'r')
legend('插值节点','三次样条',2)
axis auto
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% 计算x_nu = 3处的样条函数值部分
x_nu = 3;
% 判断x_nu所在的单元
for j = 1:n-1
if x_nu > x(j) && x_nu <= x(j+1)
break;
else
continue;
end
end
% x_nu位于第j个单元[x(j),x(j+1)]
% 计算x_nu处的样条函数值:y_nu
y_nu = y(j) + C(j,1)*(x_nu - x(j)) + C(j,2)*(x_nu - x(j))^2 + C(j,3)*(x_nu - x(j))^3;
% 打印y_nu
y_nu
上一次上机课我们做了P71例6这个问题,但是做出来的同学很少很少。好多同学都在机房忙着看基础知识,部分同学按照课本P70的程序提示来做,这里我 需要谈一下我的看法。 首先,我认为大家应该在课下提前做一些准备工作。我们授课课堂上布置了问题之后,同学们如果能在课下提前想一下程序的大体步骤或者简单写一些伪代码,那 么到了上机的时候效率必然比现在高。 其次,关于课本上的程序提示,我认为应该在自己编写出程序之后再去读。我认为写程序处理一个问题,应先用自己的想法把它编写调试实现出来(第一次实现时 程序一般不是最优的),然后再不断的改进优化(节省计算量+节省存储空间),直到满意为止。而课本上P70的程序提示中,很多地方用了我们一眼看不懂的 变量名字,它可能为了节省存储量,一会让该变量代替这个东西,一会让该变量代替另一个东西。所以说,只有我们把自己的程序实现出来,再去读课本提示时才 能看到它处理问题时的一些想法和妙用。读别人优秀的代码是一件好事情,甚至读伪代码都是好事情,可是如果不对比着自己的不足盲目的去读,只会造成比着葫 芦画葫芦,并不起很大的作用。 以上个人意见,仅供参考。 最后,附件里面贴出我给的两个Matlab程序,供大家参考。
spline_test.m内有较为详细的注释。程序没有优化到极致,尽量保持了一种阅读性较强的风格。
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% spline_test.m
% 本程序用于生成自然边界条件(第2种边界条件)下的样条曲线,并估计函数在某些点处的近似函数值
% 其中x,y为同样长度的向量,代表离散数据。我们要在给定x,y这些数据的前提下,对其进行三次样条插值
% 程序大体分为以下几个步骤:
% 1. 声明(定义)所需的系数矩阵A,右端向量d以及需要的常向量mu,lambda,h
% 2. 为系数矩阵和右端向量赋值
% 3. 求解线性方程组
% 4. 利用线性方程组求解得到的结果,计算样条表达式中的系数C,详见P65 (8.16)
% 5. 绘制样条函数图象
% 6. 计算给定点x_nu处的样条函数值,即被逼近的函数在该点处的近似值
%
% 2008.11.23 已调试,未优化
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% 清空内存,命令行,关闭所有的figure窗口
clear; clc; close all;
% 给定数据x,y,该处x,y可以在满足条件1时任意给定,只需它们是同样长度的向量即可,但是注意最后计算样条函数在
% x_nu处的函数值的时候,x_nu必需位于我们考虑的区间[x(1),x(n)]上。
% 条件1:x的每个分量按下标的顺序依次增大
x = [1 2 4 5 6 7];
y = [1 3 4 2 3 4];
% 取得x或者y的长度
n = length(x);
% 声明A,d,mu,lambda,h,并对它们所有的元素赋予初值0
% 关于它们的大小请参考课本P66 最上面的式子
A = zeros(n-2,n-2);
d = zeros(n-2,1);
mu = zeros(n-2,1);
lambda = zeros(n-2,1);
h = zeros(n-1,1);
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% 先给出第一个步长h(1),因为在下面的循环中会被用到
h(1) = x(2) - x(1);
% for循环给h,lambda,mu和右端向量d赋值
for j = 2:n-1
h(j) = x(j+1) - x(j);
lambda(j-1) = h(j)/(h(j) + h(j-1));
mu(j-1) = 1 - lambda(j-1);
d(j-1) = 6 * ((y(j+1) - y(j))/h(j) - (y(j) - y(j-1))/h(j-1))/(h(j) + h(j-1));
end
% 另一个for循环给A赋值,注意:A为三对角阵,所以赋值很容易实现,当然也可以用另外的方法。
for j = 1:n-3
A(j,j) = 2;
A(j,j+1) = lambda(j);
A(j+1,j) = mu(j+1);
end
A(n-2,n-2) = 2;
% ++++++++++++++++++++关于前面两个for循环的一点注释+++++++++++++++++++++++++
% 1.两个for循环没有合并在一起的原因在于,分开循环从语句来看比较明朗,如果合并在一起来赋值的话
% 可能需要if判断。两种赋值方式孰优孰劣,未比较过。我想这不是三次样条插值最关键的地方。
% 2.三对角矩阵A是稀疏矩阵,是否可以采用存储量更小的方式存储它?
% 当然可以,本程序中利用的存储方式,存储量大,一般三对角矩阵利用带状存储比较合适
% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% 求解以A为系数矩阵,d为右端向量的线性方程组AM=d
% 这里直接用A的逆乘d来得到结果,三对角方程的追赶法实现,还是留到后面让同学们自己来做吧
% 这里的M中每个分量仅代表第2个节点至倒数第2个节点处的二阶导数
M = A\d;
% 把第1个和最后1个节点处的二阶导数‘0’,即已知的自然边界条件放入M的首尾
% 生成我们下面要用到的M,现在的M的每个分量代表的就是每个节点处的二阶导数了
M = [0;M;0];
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% 生成样条表达式中的系数C,详见P65 (8.16)
C = zeros(n-1,3);
for j = 1:n-1
C(j,1) = (y(j+1) - y(j))/h(j) - (2*M(j) + M(j+1))*h(j)/6;
C(j,2) = M(j)/2;
C(j,3) = (M(j+1) - M(j))/(6*h(j));
end
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% 绘制样条函数图象部分
% 给定一个小步长h_test,一般0.01就可以
h_test = 0.01;
% x_test是是一个向量
% 它以x(1)与x(n)为第1个分量和最后一个分量,h_test为中间每个相邻分量之间的距离
x_test = x(1):h_test:x(n);
% 取得x_test的长度
m = length(x_test);
% 下面计算与x_test的每个分量处相应的样条函数值,存入y_test
y_test = zeros(m,1);
j = 1;
for i = 1:m
if x_test(i) <= x(j+1) % 判断x_test(i)位于哪个单元[x(j),x(j+1)]中
y_test(i) = y(j) + C(j,1)*(x_test(i) - x(j)) + C(j,2)*(x_test(i) - x(j))^2 + C(j,3)*(x_test(i) - x(j))^3;
else
j = j + 1;
y_test(i) = y(j) + C(j,1)*(x_test(i) - x(j)) + C(j,2)*(x_test(i) - x(j))^2 + C(j,3)*(x_test(i) - x(j))^3;
end
end
% 绘制插值节点和样条函数图象
plot(x,y,'*',x_test,y_test,'r')
legend('插值节点','三次样条',2)
axis auto
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% 计算x_nu = 3处的样条函数值部分
x_nu = 3;
% 判断x_nu所在的单元
for j = 1:n-1
if x_nu > x(j) && x_nu <= x(j+1)
break;
else
continue;
end
end
% x_nu位于第j个单元[x(j),x(j+1)]
% 计算x_nu处的样条函数值:y_nu
y_nu = y(j) + C(j,1)*(x_nu - x(j)) + C(j,2)*(x_nu - x(j))^2 + C(j,3)*(x_nu - x(j))^3;
% 打印y_nu
y_nu