数值分析实验报告三次样条插值
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
y 1 (B(4）- B(2) * C(N 1, 2)) /(2 B（2) * C(N 1, 1))
8. 计算{M j} 对于 j 1, 2, ,N 1 （1）计算 Mj： y2 C(j , 2) C(j , 1) * y 1 （2） h x(j 1) x(j ) （3） C(j , 3) (y 2 y 1 ) / 6 / h
Sx(i)=vpa(Sx(i));%三样条插值函数表达式 end for i=1:n disp('S(x)='); fprintf('%s (%d,%d)\n',char(Sx(i)),X(i),X(i+1)); end for i=1:n if xi>=X(i)&&xi<=X(i+1) S=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(xi-X(i))^3; end end disp('xi S'); fprintf('%d,%d\n',xi,S); return
d=U(1,k)-x(1,i); V(1,k)=((C(i,3)*d+C(i,2))*d+C(i,1))*d+y(1,i); a=2; break end if a==1 V(1,k)=10^7; end if a==2 break end end endຫໍສະໝຸດ Baidu
代码 2 对于三次样条插值第一类边界条件
function S=Threch1(X,Y,dy0,dyn,xi) %请输入 1*n 矩阵 X、Y，X 中每个元素都是插值节点，Y 中每个元素都是插值节点对应的函数值，dy0 左端点处的一阶导数， dyn 右端点处的一阶导数，S 求得的三次样条插值函数的值 %第二章 P71 例一三次样条插值 n=length(X)-1; d=zeros(n+1,1); h=zeros(1,n-1); f1=zeros(1,n-1); f2=zeros(1,n-2); for i=1:n%求函数的一阶差商 h(i)=X(i+1)-X(i); f1(i)=(Y(i+1)-Y(i))/h(i); end for i=2:n%求函数的二阶差商 f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1)); d(i)=6*f2(i); end d(1)=6*(f1(1)-dy0)/h(1); d(n+1)=6*(dyn-f1(n-1))/h(n-1);%?赋初值 A=zeros(n+1,n+1); B=zeros(1,n-1); C=zeros(1,n-1); for i=1:n-1 B(i)=h(i)/(h(i)+h(i+1)); C(i)=1-B(i); end A(1,2)=1; A(n+1,n)=1; for i=1:n+1 A(i,i)=2; end for i=2:n A(i,i-1)=B(i-1); A(i,i+1)=C(i-1); end M=A\d; syms x; for i=1:n Sx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))... +M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3); digits(4);
n 计算{Mj}j 2, ,n 1)，且存放在数组 C(N-1,3)内。 1 及样条系数 cj .1 ,cj , 2 ,cj , 3(j 1,
XX
完成日期：
2017/11/11
1. 输入： N ,x ,x(i ),y(i ),y (i )(i 0, ,N ),B(i ),(i 1, , 4)(或计算） 2. h1 x(2) x(1), p y(2) y(1) 3. 对于 j 1, 2, ,N 1 （1） h2 x(j 1) x(j ) （2） q y(j 1) y(j ) （3） h0 h1 h2 （4）计算 λj： C(j , 1) h2 / h0
2) 1 C(j , 1) （5）计算 μj： C(j ,
（6）计算 βj： C(j , 3) 6 *（q/h2 p / h1 ) / h0
（7） h1 h2 （8） p q 4. 计算 β1 及 z1： C(1, 1) B(1) / 2,C(1, 2) B(3) / 2 5. 如果 N=2 则转 7 6. 计算 βj 及 zj 对于 j 1, 2, ,n - 1 对于 j 1, 2, ,N 1 （1） p 2 C(j , 2) * C(j 1, 1) （2） C(j , 1) i C(j 1, 1) （3） C(j , 2) zj (C(j , 3) C(j , 2) * C(j 1, 2)) / p 7. 计算 z n M n
C(j , 2) y 2 / 2
C(j , 1) (y(j 1) y(j ))/h (2 * y 2 y 1 ) * h / 6
（4） y 1 y 2 7.输出：数组 C(N-1,3) （二） 三次样条插值计算 已知 y f(x )函数(xi ,f(xi ))， （i 0, 1, ,n )(a x1 x 2 x n b )，及 边 界 条 件 参 数 B(1:4) ， 给 定 x U (i ) [a,b ],i 1, ,m ) 内 ， 要 求 计 算
对于三次样条插值第二类边界条件
function [S,Sx]=Spline32(X,Y,d2y0,d2yn,xi) %请输入 1*n 矩阵 X、Y，X 中每个元素都是插值节点，Y 中每个元素都是插值节点对应的函数值，dy0 左端点处的一阶导数， dyn 右端点处的一阶导数，S 求得的三次样条插值函数的值，Sx 求得的函数 %第二章 P71 例一三次样条插值 n=length(X)-1; d=zeros(n+1,1); h=zeros(1,n-1); f1=zeros(1,n-1); f2=zeros(1,n-2); for i=1:n%求一阶差商 h(i)=X(i+1)-X(i); f1(i)=(Y(i+1)-Y(i))/h(i); end for i=2:n%求二阶差商 f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1)); d(i)=6*f2(i); end d(1)=2*d2y0; d(n+1)=2*d2yn;%赋初值 A=zeros(n+1,n+1); B=zeros(1,n-1); C=zeros(1,n-1); for i=1:n-1 B(i)=h(i)/(h(i)+h(i+1)); C(i)=1-B(i); end A(1,2)=0; A(n+1,n)=0; for i=1:n+1 A(i,i)=2; end for i=2:n A(i,i-1)=B(i-1); A(i,i+1)=C(i-1); end M=A\d; syms x; for i=1:n Sx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))... +M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3); digits(4);
Sx(i)=vpa(Sx(i)); end for i=1:n if xi>=X(i)&&xi<=X(i+1) S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(xi-X(i))^3; end end S=[xi S(2)]; Return
%计算λj C(j,1)=h2/h0; %计算μj C(j,2)=1-C(j,1); %计算 dj C(j,3)=6*(q/h2-p/h1)/h0; h1=h2; p=q; end %计算β1 C(1,1)=B(1,1)/2; %计算 z1 C(1,2)=B(1,3)/2; if N==2 %计算 zn=Mn y1=(B(1,4)-B(1,2)*C(N-1,2))/(2-B(1,2)*C(N-1,1)); else for j=2:1:N-1 %计算βj p=2-C(j,2)*C(j-1,1); C(j,1)=C(j,1)/p; %计算 zj C(j,2)=(C(j,3)-C(j,2)*C(j-1,2))/p; end %计算 zn=Mn y1=(B(1,4)-B(1,2)*C(N-1,2))/(2-B(1,2)*C(N-1,1)); end %计算{Mj}及样条系数 for j=N-1:-1:1 %计算 Mj y2=C(j,2)-C(j,1)*y1; h=x(1,j+1)-x(1,j); C(j,3)=(y1-y2)/6/h; C(j,2)=y2/2; C(j,1)=(y(1,j+1)-y(1,j))/h-(2*y2+y1)*h/6; y1=y2; end n=size(x,2);m=size(U,2); for k=1:1:m if U(1,k)<x(1,1) a=1; end if U(1,k)>x(1,n) a=1; end if U(1,k)==x(1,n) V(1,k)=y(1,n); a=2; end end i=zeros(1,n-1); for j=1:1:n-1 i(1,j)=n-j; end for i=i if U(1,k)<x(1,i) a=1; else
S(U(i))(i 1, ,m )，存放在 V(i) (i 1, ,m )内，且由上算法计算出样条系数， 存放在数组 C(N-1,3)内。 1.输入： x ,x(i ),y(i ),(i 0, ,n ),n,及U (i ),(i 1， ，m ),m 2.对于 k 1, 2, ,m （1）如果 U(k)<x(1)则转（5） （2）如果 U(k)>x(1)则转（5） （3）如果 U(k)=x(1)则 V(k)←y(1),转（6） （4）对于 i n 1,n 2, , 1 1）如果 U(k)<x(i)则转 5） 2） d ← U(k) - x(i) 3）V (k ) ((C(i, 3) * d C(i , 2)) * d C(i , 1) * d y(i ) 4）转（6） 5）Continue （5）V (k ) 107 （6）Continue 7.输出 U (i ),V (i ),(i 1, 2, ,m ) 主要程序代码 根据书中算法得出如下代码 1，然而结果不对，我自己又做出了代码 2 代码 1：
function [U,V,C] =Spline(x,y,B,U) %请输入 1*n 矩阵 x、y，x 中每个元素都是插值节点，y 中每个元素都是插值节点对应的函数值，B 为 边界函数参数，U 为给定 x； %第二章 P71 例一三次样条插值 N=size(x,2);C=zeros(N-1,3); h1=x(1,2)-x(1,1);p=y(1,2)-y(1,1); for j=2:1:N-1 h2=x(1,j+1)-x(1,j); q=y(1,j+1)-y(1,1); h0=h1+h2;
山东师范大学数学科学学院实验报告
实验课程： 数值分析引论 实验项目： 三次样条插值
姓名： XXX 学号： 2015080401XX 班级： XXX 班 专业： 数学与应用数学 指导教师： 实验目的 1、三次样条插值法的 matlab 实现； 2、用三次样条插值法程序进行插值计算。 实验内容： 一、三次样条插值法 问题分析和算法设计 三次样条插值 （一） 计算三次样条插值函数系数 已知 y f(x )函数表(xi ,f(xi ))， （i 0, 1, ,n )，及边界条件参数 B(1:4), 本算法
8. 计算{M j} 对于 j 1, 2, ,N 1 （1）计算 Mj： y2 C(j , 2) C(j , 1) * y 1 （2） h x(j 1) x(j ) （3） C(j , 3) (y 2 y 1 ) / 6 / h
Sx(i)=vpa(Sx(i));%三样条插值函数表达式 end for i=1:n disp('S(x)='); fprintf('%s (%d,%d)\n',char(Sx(i)),X(i),X(i+1)); end for i=1:n if xi>=X(i)&&xi<=X(i+1) S=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(xi-X(i))^3; end end disp('xi S'); fprintf('%d,%d\n',xi,S); return
d=U(1,k)-x(1,i); V(1,k)=((C(i,3)*d+C(i,2))*d+C(i,1))*d+y(1,i); a=2; break end if a==1 V(1,k)=10^7; end if a==2 break end end endຫໍສະໝຸດ Baidu
代码 2 对于三次样条插值第一类边界条件
function S=Threch1(X,Y,dy0,dyn,xi) %请输入 1*n 矩阵 X、Y，X 中每个元素都是插值节点，Y 中每个元素都是插值节点对应的函数值，dy0 左端点处的一阶导数， dyn 右端点处的一阶导数，S 求得的三次样条插值函数的值 %第二章 P71 例一三次样条插值 n=length(X)-1; d=zeros(n+1,1); h=zeros(1,n-1); f1=zeros(1,n-1); f2=zeros(1,n-2); for i=1:n%求函数的一阶差商 h(i)=X(i+1)-X(i); f1(i)=(Y(i+1)-Y(i))/h(i); end for i=2:n%求函数的二阶差商 f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1)); d(i)=6*f2(i); end d(1)=6*(f1(1)-dy0)/h(1); d(n+1)=6*(dyn-f1(n-1))/h(n-1);%?赋初值 A=zeros(n+1,n+1); B=zeros(1,n-1); C=zeros(1,n-1); for i=1:n-1 B(i)=h(i)/(h(i)+h(i+1)); C(i)=1-B(i); end A(1,2)=1; A(n+1,n)=1; for i=1:n+1 A(i,i)=2; end for i=2:n A(i,i-1)=B(i-1); A(i,i+1)=C(i-1); end M=A\d; syms x; for i=1:n Sx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))... +M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3); digits(4);
n 计算{Mj}j 2, ,n 1)，且存放在数组 C(N-1,3)内。 1 及样条系数 cj .1 ,cj , 2 ,cj , 3(j 1,
XX
完成日期：
2017/11/11
1. 输入： N ,x ,x(i ),y(i ),y (i )(i 0, ,N ),B(i ),(i 1, , 4)(或计算） 2. h1 x(2) x(1), p y(2) y(1) 3. 对于 j 1, 2, ,N 1 （1） h2 x(j 1) x(j ) （2） q y(j 1) y(j ) （3） h0 h1 h2 （4）计算 λj： C(j , 1) h2 / h0
2) 1 C(j , 1) （5）计算 μj： C(j ,
（6）计算 βj： C(j , 3) 6 *（q/h2 p / h1 ) / h0
（7） h1 h2 （8） p q 4. 计算 β1 及 z1： C(1, 1) B(1) / 2,C(1, 2) B(3) / 2 5. 如果 N=2 则转 7 6. 计算 βj 及 zj 对于 j 1, 2, ,n - 1 对于 j 1, 2, ,N 1 （1） p 2 C(j , 2) * C(j 1, 1) （2） C(j , 1) i C(j 1, 1) （3） C(j , 2) zj (C(j , 3) C(j , 2) * C(j 1, 2)) / p 7. 计算 z n M n
C(j , 2) y 2 / 2
C(j , 1) (y(j 1) y(j ))/h (2 * y 2 y 1 ) * h / 6
（4） y 1 y 2 7.输出：数组 C(N-1,3) （二） 三次样条插值计算 已知 y f(x )函数(xi ,f(xi ))， （i 0, 1, ,n )(a x1 x 2 x n b )，及 边 界 条 件 参 数 B(1:4) ， 给 定 x U (i ) [a,b ],i 1, ,m ) 内 ， 要 求 计 算
对于三次样条插值第二类边界条件
function [S,Sx]=Spline32(X,Y,d2y0,d2yn,xi) %请输入 1*n 矩阵 X、Y，X 中每个元素都是插值节点，Y 中每个元素都是插值节点对应的函数值，dy0 左端点处的一阶导数， dyn 右端点处的一阶导数，S 求得的三次样条插值函数的值，Sx 求得的函数 %第二章 P71 例一三次样条插值 n=length(X)-1; d=zeros(n+1,1); h=zeros(1,n-1); f1=zeros(1,n-1); f2=zeros(1,n-2); for i=1:n%求一阶差商 h(i)=X(i+1)-X(i); f1(i)=(Y(i+1)-Y(i))/h(i); end for i=2:n%求二阶差商 f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1)); d(i)=6*f2(i); end d(1)=2*d2y0; d(n+1)=2*d2yn;%赋初值 A=zeros(n+1,n+1); B=zeros(1,n-1); C=zeros(1,n-1); for i=1:n-1 B(i)=h(i)/(h(i)+h(i+1)); C(i)=1-B(i); end A(1,2)=0; A(n+1,n)=0; for i=1:n+1 A(i,i)=2; end for i=2:n A(i,i-1)=B(i-1); A(i,i+1)=C(i-1); end M=A\d; syms x; for i=1:n Sx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))... +M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3); digits(4);
Sx(i)=vpa(Sx(i)); end for i=1:n if xi>=X(i)&&xi<=X(i+1) S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(xi-X(i))^3; end end S=[xi S(2)]; Return
%计算λj C(j,1)=h2/h0; %计算μj C(j,2)=1-C(j,1); %计算 dj C(j,3)=6*(q/h2-p/h1)/h0; h1=h2; p=q; end %计算β1 C(1,1)=B(1,1)/2; %计算 z1 C(1,2)=B(1,3)/2; if N==2 %计算 zn=Mn y1=(B(1,4)-B(1,2)*C(N-1,2))/(2-B(1,2)*C(N-1,1)); else for j=2:1:N-1 %计算βj p=2-C(j,2)*C(j-1,1); C(j,1)=C(j,1)/p; %计算 zj C(j,2)=(C(j,3)-C(j,2)*C(j-1,2))/p; end %计算 zn=Mn y1=(B(1,4)-B(1,2)*C(N-1,2))/(2-B(1,2)*C(N-1,1)); end %计算{Mj}及样条系数 for j=N-1:-1:1 %计算 Mj y2=C(j,2)-C(j,1)*y1; h=x(1,j+1)-x(1,j); C(j,3)=(y1-y2)/6/h; C(j,2)=y2/2; C(j,1)=(y(1,j+1)-y(1,j))/h-(2*y2+y1)*h/6; y1=y2; end n=size(x,2);m=size(U,2); for k=1:1:m if U(1,k)<x(1,1) a=1; end if U(1,k)>x(1,n) a=1; end if U(1,k)==x(1,n) V(1,k)=y(1,n); a=2; end end i=zeros(1,n-1); for j=1:1:n-1 i(1,j)=n-j; end for i=i if U(1,k)<x(1,i) a=1; else
S(U(i))(i 1, ,m )，存放在 V(i) (i 1, ,m )内，且由上算法计算出样条系数， 存放在数组 C(N-1,3)内。 1.输入： x ,x(i ),y(i ),(i 0, ,n ),n,及U (i ),(i 1， ，m ),m 2.对于 k 1, 2, ,m （1）如果 U(k)<x(1)则转（5） （2）如果 U(k)>x(1)则转（5） （3）如果 U(k)=x(1)则 V(k)←y(1),转（6） （4）对于 i n 1,n 2, , 1 1）如果 U(k)<x(i)则转 5） 2） d ← U(k) - x(i) 3）V (k ) ((C(i, 3) * d C(i , 2)) * d C(i , 1) * d y(i ) 4）转（6） 5）Continue （5）V (k ) 107 （6）Continue 7.输出 U (i ),V (i ),(i 1, 2, ,m ) 主要程序代码 根据书中算法得出如下代码 1，然而结果不对，我自己又做出了代码 2 代码 1：
function [U,V,C] =Spline(x,y,B,U) %请输入 1*n 矩阵 x、y，x 中每个元素都是插值节点，y 中每个元素都是插值节点对应的函数值，B 为 边界函数参数，U 为给定 x； %第二章 P71 例一三次样条插值 N=size(x,2);C=zeros(N-1,3); h1=x(1,2)-x(1,1);p=y(1,2)-y(1,1); for j=2:1:N-1 h2=x(1,j+1)-x(1,j); q=y(1,j+1)-y(1,1); h0=h1+h2;
山东师范大学数学科学学院实验报告
实验课程： 数值分析引论 实验项目： 三次样条插值
姓名： XXX 学号： 2015080401XX 班级： XXX 班 专业： 数学与应用数学 指导教师： 实验目的 1、三次样条插值法的 matlab 实现； 2、用三次样条插值法程序进行插值计算。 实验内容： 一、三次样条插值法 问题分析和算法设计 三次样条插值 （一） 计算三次样条插值函数系数 已知 y f(x )函数表(xi ,f(xi ))， （i 0, 1, ,n )，及边界条件参数 B(1:4), 本算法