三次B样条反求控制点

三次B样条反求控制点
三次B样条反求控制点

三次样条插值代码

2 三次样条插值程序 三次样条插值利用方案二(求解固支样条或压紧样条) 按照要求要起点和终点的一阶导数值已知, 可得关于01,,.....,n M M M 的严格对角占优势的三对角方程组 然后利用三对角法(追赶法)解此线性方程组。 (1)编写M 文件,并保存文件名scfit.m % x,y 分别为n 个节点的横坐标和纵坐标值组成的向量 % dx0和dxn 分别为S 的导数在x0和xn 处的值,即m 0和m n n=length(x)-1; h=diff(x); d=diff(y)./h; a=h(2:n-1); b=2*(h(1:n-1)+h(2:n)); c=h(2:n); u=6*diff(d); b(1)=b(1)-h(1)/2; u(1)=u(1)-3*(d(1)-dx0); b(n-1)=b(n-1)-h(n)/2; u(n-1)=u(n-1)-3*(dxn-d(n)); %追赶法部分 for k=2:n-1 temp=a(k-1)/b(k-1); b(k)=b(k)-temp*c(k-1); u(k)=u(k)-temp*u(k-1); end m(n)=u(n-1)/b(n-1); for k=n-2:-1:1 m(k+1)=(u(k)-c(k)*m(k+2))/b(k); end %求S K1,S K2,S K3,S K4 m(1)=3*(d(1)-dx0)/h(1)-m(2)/2; m(n+1)=3*(dxn-d(n))/h(n)-m(n)/2; for k=0:n-1 00 ()S x m '=()n n S x m '=0011111111212212n n n n n n M d M d M d M d μλμλ----??????????????????????=??????????????????????????

三次样条插值函数

沈阳航空航天大学 数学软件课程设计 (设计程序) 题目三次样条插值函数 班级 / 学号 学生姓名 指导教师

沈阳航空航天大学 课程设计任务书 课程名称数学软件课程设计 院(系)理学院专业信息与计算科学 班级学号姓名 课程设计题目三次样条插值函数 课程设计时间: 2010 年12月20日至2010 年12月31日 课程设计的内容及要求: 1.三次样条插值函数 给出函数在互异点处的值分别为。 (1)掌握求三次样条插值函数的基本原理; (2)编写程序求在第一边界条件下函数的三次样条插值函数; (3)在区间上取n=10,20,分别用等距节点对函数 作三次样条插值函数,利用(1)的结果画出插值函数的图形,并在该图形界面中同时画出的图形。 [要求] 1.学习态度要认真,要积极参与课程设计,锻炼独立思考能力; 2.严格遵守上机时间安排; 3.按照MATLAB编程训练的任务要求来编写程序; 4.根据任务书来完成课程设计论文; 5.报告书写格式要求按照沈阳航空航天大学“课程设计报告撰写规范”; 6.报告上交时间:课程设计结束时上交报告;

7.严谨抄袭行为。 指导教师年月日负责教师年月日学生签字年月日

沈阳航空航天大学 课程设计成绩评定单 课程名称数学软件课程设计 院(系)理学院专业信息与计算科学课程设计题目三次样条插值函数 学号姓名 指导教师评语: 课程设计成绩 指导教师签字 年月日

目录 一正文 (1) 1问题分析 (1) 1.1 题目 (1) 1.2 分析 (1) 2 研究方法原理 (1) 2.1 求三次样条插值多项式,算法组织 (1) 3 算例结果 (3) 二总结 (7) 参考文献 (8) 附录 (9) 源程序: (9) 程序1 (9) 程序2 (10) 程序3 (12) 程序 4 (12)

数值分析作业-三次样条插值

数值计算方法作业 实验4.3 三次样条差值函数 实验目的: 掌握三次样条插值函数的三弯矩方法。 实验函数: dt e x f x t ? ∞ -- = 2 221)(π 实验内容: (1) 编程实现求三次样条插值函数的算法,分别考虑不同的边界条件; (2) 计算各插值节点的弯矩值; (3) 在同一坐标系中绘制函数f(x),插值多项式,三次样条插值多项式的曲线 比较插值结果。 实验4.5 三次样条差值函数的收敛性 实验目的: 多项式插值不一定是收敛的,即插值的节点多,效果不一定好。对三次样条插值函数如何呢?理论上证明三次样条插值函数的收敛性是比较困难的,通过本实验可以证明这一理论结果。 实验内容: 按照一定的规则分别选择等距或非等距的插值节点,并不断增加插值节点的个数。 实验要求: (1) 随着节点个数的增加,比较被逼近函数和三样条插值函数的误差变化情 况,分析所得结果并与拉格朗日插值多项式比较; (2) 三次样条插值函数的思想最早产生于工业部门。作为工业应用的例子,考

虑如下例子:某汽车制造商根据三次样条插值函数设计车门曲线,其中一 算法描述: 拉格朗日插值: 错误!未找到引用源。 其中错误!未找到引用源。是拉格朗日基函数,其表达式为:() ∏ ≠=--=n i j j j i j i x x x x x l 0) ()( 牛顿插值: ) )...()(](,...,,[.... ))(0](,,[)0](,[)()(1102101210100----++--+-+=n n n x x x x x x x x x x f x x x x x x x f x x x x f x f x N 其中????? ?? ?? ?????? --=--= --= -)/(]),...,[],...,[(]...,[..],[],[],,[)()(],[01102110x x x x x f x x x f x x x f x x x x f x x f x x x f x x x f x f x x f n n n n i k j i k j k j i j i j i j i 三样条插值: 所谓三次样条插值多项式Sn(x)是一种分段函数,它在节点Xi(a

(精选)三次样条插值的MATLAB实现

MATLAB 程序设计期中考查 在许多问题中,通常根据实验、观测或经验得到的函数表或离散点上的信息,去研究分析函数的有关特性。其中插值法是一种最基本的方法,以下给出最基本的插值问题——三次样条插值的基本提法: 对插值区间[]b a ,进行划分:b x x x a n ≤

三次样条插值方法的应用

CENTRAL SOUTH UNIVERSITY 数值分析实验报告

三次样条插值方法的应用 一、问题背景 分段低次插值函数往往具有很好的收敛性,计算过程简单,稳定性好,并且易于在在电子计算机上实现,但其光滑性较差,对于像高速飞机的机翼形线船体放样等型值线往往要求具有二阶光滑度,即有二阶连续导数,早期工程师制图时,把富有弹性的细长木条(即所谓的样条)用压铁固定在样点上,在其他地方让他自由弯曲,然后沿木条画下曲线,称为样条曲线。样条曲线实际上是由分段三次曲线并接而成,在连接点即样点上要求二阶导数连续,从数学上加以概括就得到数学样条这一概念。下面我们讨论最常用的三次样条函数及其应用。 二、数学模型 样条函数可以给出光滑的插值曲线(面),因此在数值逼近、常微分方程和偏微分方程的数值解及科学和工程的计算中起着重要的作用。 设区间[]b ,a 上给定有关划分b x x n =<<<=Λ10x a ,S 为[]b ,a 上满足下面条件的函数。 ● )(b a C S ,2∈; ● S 在每个子区间[]1,+i i x x 上是三次多项式。 则称S 为关于划分的三次样条函数。常用的三次样条函数的边界条件有三种类型: ● Ⅰ型 ()()n n n f x S f x S ''0'',==。 ● Ⅱ型 ()()n n n f x S f x S ''''0'''',==,其特殊情况为()()0''''==n n x S x S 。 ● Ⅲ型 ()()Λ3,2,1,0,0==j x S x S n j j ,此条件称为周期样条函数。

鉴于Ⅱ型三次样条插值函数在实际应用中的重要地位,在此主要对它进行详细介绍。 三、算法及流程 按照传统的编程方法,可将公式直接转换为MATLAB可是别的语言即可;另一种是运用矩阵运算,发挥MATLAB在矩阵运算上的优势。两种方法都可以方便地得到结果。方法二更直观,但计算系数时要特别注意。这里计算的是方法一的程序,采用的是Ⅱ型边界条件,取名为spline2.m。 Matlab代码如下: function s=spline2(x0,y0,y21,y2n,x) %s=spline2(x0,y0,y21,y2n,x) %x0,y0 are existed points,x are insert points,y21,y2n are the second %dirivitive numbers given. n=length(x0); km=length(x); a(1)=-0.5; b(1)=3*(y0(2)-y0(1))/(2*(x0(2)-x0(1))); for j=1:(n-1) h(j)=x0(j+1)-x0(j); end for j=2:(n-1) alpha(j)=h(j-1)/(h(j-1)+h(j)); beta(j)=3*((1-alpha(j))*(y0(j)-y0(j-1))/h(j-1)+alpha(j)*(y0(j+1)-y0(j))/h(j));

三次样条函数

计算方法实验报告 1、实验题目 三次样条插函数。 2、实验内容 三次样条插值是建立在Hermite 插值的基础上的。Hermite 插值是在一个区间上的插值,而三次样条插则是建立多个区间上插值,构造一个具有二阶光滑度的曲线,在求出给定点上对应的函数。本实验就是建立一个能根据三次样条插值函数求根的程序。 3、算法思想 给定一个区间,并把它分成n 等份,并且给出了每个结点对就的横坐标和纵坐标。利用程序输出给定插值点对应的值。横坐标设为:X 0, X 1, X 2, X 3, …X n 纵坐标为Y 0, Y 1, Y 2, …Y n ,设插点为u 。则令h k =X k+1-X k ,λk =1-+k k k h h h , μk =11--+k k k h h h , g k =3(1 11--+-+-k k k k k k k k h y y h y y λμ), 其中k=1,2,…,n-1 再根据第一类边界条件则可以确定公式6.16,再根据6.17解出方程中的m 向量,最后代入公式6.8求解。 4、源程序清单 #include #define N 21/*最大结点个数减一*/ void sanCi() { /*定义过程数据变量*/ float x[N],y[N],h[N]; /*横纵坐标及区间长度*/ float rr[N],uu[N],gg[N]; /*计算m 用的中间数组rr 、uu 、gg 分别对应:λ、μ、g 数组*/

float aa[N],bb[N],tt[N]; /*矩阵分解时用到的中间变量aa、bb、tt分别对应:α、β数组以及A=LU时中间矩阵*/ float mm[N]; /*最后要用到的系数m*/ int n,k,kv,chose; /* n为实际结点个数,k为下标,kv为最后确定k的值*/ float s,u; /*最后计算u对应的值*/ printf("请输入区间段数:"); scanf("%d",&n); /*输入结点个数*/ /*输入所有横坐标:*/ printf("输入所有横坐标:"); for(k=0; k<=n; k++) scanf("%f",&x[k]); /*输入对应纵坐标:*/ printf("输入对应纵坐标:"); for(k=0; k<=n; k++) scanf("%f",&y[k]); for(k=0; k

MATLAB三次样条插值之三弯矩法

MATLAB三次样条插值之三弯矩法 首先说这个程序并不完善,为了实现通用(1,2,…,n)格式解题,以及为调用追赶法程序,没有针对节点数在三个以下的情况进行分类讨论。希望能有朋友给出更好的方法。 首先,通过函数sanwanj得到方程的系数矩阵,即追赶法方程的四个向量参数,接下来调 用追赶法(在intersanwj函数中),得到三次样条分段函数系数因子,然后进行多项式合并 得到分段函数的解析式,程序最后部分通过判断输入值的区间自动选择对应的分段函数并计算 改点的值。附:追赶法程序chase %%%%%%%%%%%%%% function [newv,w,newu,newd]=sanwj(x,y,x0,y0,y1a,y1b)?%三弯矩样 条插值?%将插值点分两次输入,x0y0单独输入?% 边值条件a的二阶导数 y1a 和b 的二阶导数y1b,这里建议将y1a和y1b换成y2a和y2b,以便于和三转角代码相区别 ?n=length(x);m=length(y); if m~=n?error('x or y 输入有误,再来'); end?v=ones(n-1,1);u=ones(n-1,1);d=zeros(n-1,1);?w=2*o nes(n+1);?h0=x(1)-x0;?h=zeros(n-1,1); for k=1:n-1?h(k)=x(k+1)-x(k);?end v(1)=h0/(h0+h(1)); u(1)=1-v(1); d(1)=6*((y(2)-y(1))/h(1)-(y(1)-y0)/h0)/(h0+h(1));?% for k=2:n-1?v(k)=h(k-1)/(h(k-1)+h(k));?u(k)=1-v(k);?d(k)= 6*((y(k+1)-y(k))/h(k)-(y(k)-y(k-1))/h(k-1))/(h(k-1)+h(k)); end newv=[v;1];?newu=[1;u]; d0=6*((y(1)-y0)/h0-y1a)/h0; d(n)=6*(y1b-(y(n)-y(n-1))/h(n-1))/h(n-1); newd=[d0;d]; %%%%%%%%%%%% function intersanwj(x,y,x0,y0,y1a,y1b) %三弯矩样条插值?%第一部分?n=length(x);m=length(y); if m~=n?error('xory 输入有误,再来'); end?%重新定义h?h=zeros(n,1); h(1)=x(1)-x0; for k=2:n h(k)=x(k)-x(k-1);?end %sptep1调用三弯矩函数?[a,b,c,d]=sanwj(x,y,x0,y0,y1a,y1b);

实验四三次样条插值Word版

实验四三次样条插值的应用 一、问题描述 The upper portion of this noble beast is to be approximated using clamped cubic spline interpolants. The curve is drawn on a grid from which the table is constructed. Use Algorithm 3.5 to construct the three clamped cubic splines. 二、模型建立 三次样条插值 给定一个列表显示的函数 yi=y(xi),i=0,1,2,...,N-1。特别注意在xj和xj+1之间的一个特殊的区间。该区间的线性插值公式为:

(3.3.1)式和(3.3.2)式是拉格朗日插值公式(3.1.1)的特殊情况。 因为它是(分段)线性的,(3.3.1)式在每一区间内的二阶导数为零,在横坐标为xj处的二阶导数不定义或无限。三次样条插值的目的就是要得到一个内插公式,不论在区间内亦或其边界上,其一阶导数平滑,二阶导数连续。 做一个与事实相反的个假设,除yi的列表值之外,我们还有函数二阶导数y"的列表值,即一系列的yi"值,则在每个区间内,可以在(3.3.1)式的右边加上一个三次多项式,其二阶导数从左边的yj"值线性变化到右边的yj+1"值,这么做便得到了所需的连续二阶导数。如果还将三次多项式构造在xj和xj+1处为零,则不会破坏在终点xj和xj+1处与列表函数值yj和yj+1的一致性。 进行一些辅助计算便可知,仅有一种办法才能进行这种构造,即用 注意,(3.3.3)式和(3.3.4)式对自变量x的依赖,是完全通过A和B对x的线性依赖,以及C和D(通过A和B)对x的三次依赖而实现。可以很容易地验证,y"事实上是该插值多项式的二阶导数。使用ABCD的定义对x求(3.3.3)式的导数,计算dA/dx dB/dx dC/dx dD/dx,结果为一阶导数

三次样条插值的C程序(很全啊)

三次样条插值C/C++程序(自己整理的) 具体推导看书<<数值分析>> code: #include using namespace std; const int MAXN = 100; int n; double x[MAXN], y[MAXN]; //下标从0..n double alph[MAXN], beta[MAXN], a[MAXN], b[MAXN]; double h[MAXN]; double m[MAXN]; //各点的一阶导数; inline double sqr(double pa) { return pa * pa; } double sunc(double p, int i) { return (1 + 2 * (p - x[i]) / (x[i + 1] - x[i])) * sqr((p - x[i + 1]) / (x[i + 1] - x[i])) * y[i] + (1 + 2 * (p - x[i + 1]) / (x[i] - x[i + 1])) * sqr((p - x[i]) / (x[i + 1] - x[i])) * y[i + 1] + (p - x[i]) * sqr((p - x[i + 1]) / (x[i] - x[i + 1])) * m[i] + (p - x[i + 1]) * sqr((p - x[i]) / (x[i + 1] - x[i])) * m[i + 1]; } int main() { int i, j;

freopen("threeInsert.in", "r", stdin); scanf("%d", &n); for (i = 0; i <= n; i++) scanf("%lf%lf", &x[i], &y[i]); // scanf("%lf%lf", &m[0], &m[n]); for (i = 0; i <= n - 1; i++) h[i] = x[i + 1] - x[i]; //第一种边界条件 //alph[0] = 0; alph[n] = 1; beta[0] = 2 * m[0]; beta[n] = 2 * m[n]; //第二种边界条件 alph[0] = 1; alph[n] = 0; beta[0] = 3 * (y[1] - y[0]) / h[0]; beta[n] = 3 * (y[n] - y[n - 1] / h[n - 1]); for (i = 1; i <= n - 1; i++) { alph[i] = h[i - 1] / (h[i - 1] + h[i]); beta[i] = 3 * ((1 - alph[i]) * (y[i] - y[i - 1]) / h[i - 1] + alph[i] * (y[i + 1] - y[i]) / h[i]); } a[0] = - alph[0] / 2; b[0] = beta[0] / 2; for (i = 1; i <= n; i++) { a[i] = - alph[i] / (2 + (1 - alph[i]) * a[i - 1]); b[i] = (beta[i] - (1 - alph[i]) * b[i - 1]) / (2 + (1 - alph[i]) * a[i - 1]); } m[n + 1] = 0; for (i = n; i >= 0; i--) { m[i] = a[i] * m[i + 1] + b[i]; } scanf("%lf", &xx); for (i = 0; i < n; i++) { if (xx >= x[i] && xx <= x[i + 1]) break ; } printf("%lf\n", sunc(xx, i));

三次样条程序

function [a b f]=spline3(x,y,flag,vl,vr) %三次样条插值函数 %(x,y)为插值节点,xx为插值点; %flag表端点边界条件类型: %flag=1:第一类边界条件(端点一阶导数给定); %flag=2:第二类边界条件(端点二阶导数给定); %flag=3:第三类边界条件; %vl,vr表左右端点处的在边界条件值。 %样条函数为:Si(x)=yi+bi*(x-xi)+ci*(x-xi)^2+di*(x-xi)^3 %b,c,d分别为各子区间上的系数值 %yy表插值点处的函数值. if length(x)~=length(y) error('输入数据应成对!'); end n=length(x); a=zeros(n-1,1); b=zeros(n-1,1); dx=a;dy=a; A=zeros(n);B=zeros(n,1); for i=1:n-1 dx(i)=x(i+1)-x(i); dy(i)=y(i+1)-y(i); end for i=2:n-1 A(i,i-1)=dx(i-1)/(dx(i-1)+dx(i)); A(i,i)=2; A(i,i+1)=dx(i)/(dx(i-1)+dx(i)); B(i,1)=6*(dy(i)/dx(i)-dy(i-1)/dx(i-1))/(dx(i)+dx(i-1)); end %---------------------------------% %端点一阶导数条件% if flag==1 A(1,1)=2; A(1,2)=1; A(n,n-1)=1; A(n,n)=2; B(1,1)=6*(dy(1)/dx(1)-vl)/dx(1); B(n,1)=6*(vr-dy(n-1)/dx(n-1))/dx(n-1); c=A\B; end %---------------% %端点二阶导数条件% if flag==2 A(1,1)=2; A(n,n)=2; B(1,1)=2*vl; B(n,1)=2*vr; c=A\B; end

三次样条插值多项式matlab

三次样条插值多项式 ——计算物理实验作业四 陈万物理学2013级 主程序: clear,clc; format rat x = [1,4,9,16,25,36,49,64]; y = [1,2,3,4,5,6,7,8]; f1 = ; fn = 1/16; [a,b,c,d,M,S] = spline(x,y,f1,fn); 子程序1: function [a,b,c,d,M,S]=spline(x,y,f1,fn) % 三次样条插值函数 % x是插值节点的横坐标 % y是插值节点的纵坐标 % u是插值点的横坐标 % f1是左端点的一阶导数 % fn是右端点的一阶导数 % a是三对角矩阵对角线下边一行 % b是三对角矩阵对角线 % c是三对角矩阵对角线上边一行 % S是插值点的纵坐标

n = length(x); h = zeros(1,n-1); deltay = zeros(1,n); miu = zeros(1,n-1); lamda = zeros(1,n-1); d = zeros(1,n-1); for j = 1:n-1 h(j) = x(j+1)-x(j); deltay(j) = y(j+1)-y(j); end % 得到h矩阵 for j = 2:n-1 sumh = h(j-1) + h(j); miu(j) = h(j-1) / sumh; lamda(j) = h(j) / sumh; d(j) = 6*( deltay(j)/h(j)-(deltay(j-1)/h(j-1)))/sumh; end % 根据第一类边界条件,作如下规定 lamda(1) = 1; d(1) = 6*(deltay(1)/h(1)-f1)/h(1); miu(1) = 1; d(n) = 6*(fn-deltay(n-1)/h(n-1))/h(n-1);

三次样条插值MATLAB程序及结果展示

23、汽车门曲线三次样条插值曲线相关程序以及结果 原始数据点: x = 0:10; %取自变量为1,2,3, (10) y = [2.51 3.30 4.04 4.70 5.22 5.54 5.78 5.40 5.57 5.70 5.80]; %输入因变量y的值 xx = linspace(min(x),max(x),200); %在x的上下界之间取200个插值节点 pp = csape(x,y,'comlete',[0.8,0.2]); %分段三次样条插值,边界条件为左右端点的一阶导数为0.8和0.2 yy = ppval(pp,xx);%计算200个插值节点对应的y值 plot(x,y,'ko',xx,yy,'k') %画出给定的11个点以及插值函数的图像

24、飞鸟外形上部自然边界条件的三次样条插值曲线相关程序以及结果 原始数据如下: x =[0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0 7.0 8.0 9.2 10.5 11.3 11.6 12.0 12.6 13.0 13.3]; y = [1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25];

xx = linspace(min(x),max(x),200); pp = csape(x,y,'second'); %分段三次样条插值,边界条件为左右端点的二阶导数为0,也称为自然边界条件 yy = ppval(pp,xx); plot(x,y,'ko',xx,yy,'k')

MATLAB三次样条插值之三转角法

非常类似前面的三弯矩法,这里的sanzhj函数和intersanzhj作用相当于前面的sanwanj和intersanwj,追赶法程序通用,代码如下。 %%%%%%%%%%%%%%%%%%% function [newu,w,newv,d]=sanzhj(x,y,x0,y0,y1a,y1b) %三转角样条插值 %将插值点分两次输入,x0 y0单独输入 %边值条件a的一阶导数y1a 和b的一阶导数y1b n=length(x);m=length(y); if m~=n error('x or y 输入有误,再来'); end v=ones(n-1,1); u=ones(n-1,1); d=zeros(n-1,1); w=2*ones(n-1,1); h0=x (1)-x0; h=zeros(n-1,1); for k=1:n-1 h(k)=x(k+1)-x(k); end v(1)=h0/(h0+h

(1));u(1)=1-v (1);d(1)=3*(v (1)*(y (2)-y (1))/h (1)+u (1)*((y (1)-y0))/h0);%for k=2:n-1 v(k)=h(k-1)/(h(k-1)+h(k)); u(k)=1-v(k); d(k)=3*(v(k)*(y(k+1)-y(k))/h(k)+u(k)*(y(k)-y(k-1))/h(k-1)); end d(1)=d (1)-u (1)*y1a; d(n-1)=d(n-1)-v(n-1)*y1b; newv=v(1:n-2,: ); newu=u(2:n-1,: ); %%%%%%%%%%%% function intersanzhj(x,y,x0,y0,y1a,y1b)

紧压三次样条插值程序讲解

求压紧三次样条函数的函数程序 function S=liti05_7(X,Y,dx0,dxn) %Input -X is the 1xn abscissa vector % -Y is the 1xn ordinate vector % -dx0=S'(x0) first derivative boundary condition % -dxn=S'(xn) first derivative boundary condition %Output -S:rows of S are the confficients,indescending order,for the cubic interplants N=length(X)-1; %求当前问题的规模数--x的最大下标 H=diff(X); %求X的差分 h0 h1…… h n-1 D=diff(Y)./H; %求 y 对 x 的一阶差商 d0 d1…… d n-1 A=H(2:N-1); %为求线性方程组系数矩阵上次对角元做准备 B=2*(H(1:N-1)+H(2:N)); %求线性方程组系数矩阵主对角元做准备 C=H(2:N-1); %求线性方程组系数矩阵上次对角元 U=6*diff(D); %为求线性方程组常数列向量做准备 %压紧样条端点约束 B(1)=B(1)-H(1)/2; %3h0/2+2h1 U(1)=U(1)-3*(D(1)-dx0); %修改 u(1) B(N-1)=B(N-1)-H(N)/2; % 2h(N-2)+3*h(N-1)/2 U(N-1)=U(N-1)-3*(dxn-D(N)); %修改 u(N-1) A %输出线性方程组系数矩阵上次对角元向量 B %输出线性方程组系数矩阵主对角元向量 C %输出线性方程组系数矩阵下次对角元向量 U %输出线性方程组常数列向量 % 下面开始解三对角线行方程组 % 首先转化为上三角线性方程组 for k=2:N-1 temp=A(k-1)/B(k-1); B(k)=B(k)-temp*C(k-1); U(k)=U(k)-temp*U(k-1); end %回代求解 M(N)=U(N-1)/B(N-1); for k=N-2:-1:1 M(k+1)=(U(k)-C(k)*M(k+2))/B(k); end %计算端点x0 x n上的二阶导数 M(1)=3*(D(1)-dx0)/H(1)-M(2)/2; M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2; M % 输出各节点上的二阶导数 for k=0:N-1 %计算第k个多项式的系数

Matlab程序三次样条插值函数

已知一组数据点,编写一程序求解三次样条插值函数满足 并针对下面一组具体实验数据 求解,其中边界条件为. 解:Matlab计算程序为: clear clc x=[0.25 0.3 0.39 0.45 0.53] y=[0.5000 0.5477 0.6245 0.6708 0.7280] n=length(x); for i=1:n-1 h(i)=x(i+1)-x(i); end for i=1:n-2 k(i)=h(i+1)/(h(i)+h(i+1)); u(i)=h(i)/(h(i)+h(i+1)); end for i=1:n-2 gl(i)=3*(u(i)*(y(i+2)-y(i+1))/h(i+1)+k(i)*(y(i+1)-y(i))/h(i)); end g0=3*(y(2)-y(1))/h(1); g00=3*(y(n)-y(n-1))/h(n-1); g=[g0 gl g00]; g=transpose(g) k1=[k 1]; u1=[1 u]; Q=2*eye(5)+diag(u1,1)+diag(k1,-1) m=transpose(Q\g) syms X; for i=1:n-1 p1(i)=(1+2*(X-x(i))/h(i))*((X-x(i+1))/h(i))^2*y(i); p2(i)=(1-2*(X-x(i+1))/h(i))*((X-x(i))/h(i))^2*y(i+1); p3(i)=(X-x(i))*((X-x(i+1))/h(i))^2*m(i); p4(i)=(X-x(i+1))*((X-x(i))/h(i))^2*m(i+1);

相关文档
最新文档