数值计算讨论组

欢迎来到数值计算讨论组。 本论坛旨在增强师生交流的时效性,主要讨论数值计算相关问题,包括初等数值分析、矩阵计算、非线性方程(组)数值解法、微分方程数值解法等方面的内容。


    [数] 三次样条插值的参考程序

    分享

    Admin
    Admin

    帖子数 : 46
    注册日期 : 10-05-13

    [数] 三次样条插值的参考程序

    帖子 由 Admin 于 周三 十二月 08, 2010 9:19 pm

    本实验项目的要求是让大家通过求解以节点处的一阶导数值为未知量的线性方程组,实现对自然边值条件下的三次样条函数的确定。这个问题看上很麻烦,但是只要把步骤弄清楚,参考我在课堂上给大家的程序,小心仔细的设计那些看似复杂的参数,求解三对角方程即可得到结果。大部分同学完成得很好。我们建议大家在学习完消元法求解线性方程组的知识后,再回头来看一下本程序中线性方程组的求解过程,尝试自己写一个求解三对角方程的程序或函数。
    示范程序如下:
    代码:
    %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    % spline_student.m
    % 本程序用于生成自然边界条件(第2种边界条件)下的样条曲线,并估计函数在某些点处的近似函数值
    % 其中x,y为同样长度的向量,代表离散数据。我们要在给定x,y这些数据的前提下,对其进行三次样条插值
    % 程序大体分为以下几个步骤:
    % 1. 声明(定义)所需的系数矩阵A,右端向量g以及需要的常向量mu,lambda,h
    % 2. 为系数矩阵和右端向量赋值
    % 3. 求解线性方程组
    % 4. 利用线性方程组求解得到的结果,计算样条表达式中的系数C.
    % 5. 绘制样条函数图象
    % 6. 计算给定点x_nu处的样条函数值,即被逼近的函数在该点处的近似值
    %
    % note: 本程序是通过求解节点处一阶导数值相关的线性方程组(三对角)实现的
    %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    % 清空内存,命令行,关闭所有的figure窗口
    % clear; clc; close all;
    clear; close all;

    x = [1 2 4 5 6 7];
    y = [1 3 4 2 3 4];

    n = length(x);
    A = zeros(n,n);
    g = zeros(n,1);
    mu = zeros(n-2,1);
    lambda = zeros(n-2,1);
    h = zeros(n-1,1);
    h(1) = x(2) - x(1);
    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);
        g(j) = 3*(mu(j-1)*(y(j+1) - y(j))/h(j)+lambda(j-1)*(y(j) - y(j-1))/h(j-1));
    end

    for j = 2:n-1
        A(j,j-1) = lambda(j-1);
        A(j,j) = 2;
        A(j,j+1) = mu(j-1);
    end
    A(1,1) = 2;A(1,2)=1;A(n,n)=2;A(n,n-1)=1;
    g(1) = 3*(y(2)-y(1))/h(1);
    g(n) = 3*(y(n)-y(n-1))/h(n-1);

    m = A\g;

    C = zeros(n-1,3);
    for j = 1:n-1
        C(j,1) = m(j);
        C(j,2) = (3*(y(j+1)-y(j))/h(j)-2*m(j)-m(j+1))/h(j);
        C(j,3) = (m(j+1) + m(j)-2*(y(j+1)-y(j))/h(j))/h(j)^2;
    end

    h_test = 0.01;
    x_test = x(1):h_test:x(n);
    m = length(x_test);
    y_test = zeros(m,1);
    j = 1;
    for i = 1:m
        if x_test(i) <= 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;
    for j = 1:n-1
        if x_nu > x(j) && x_nu <= x(j+1)
            break;
        else
            continue;
        end
    end
    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



      目前的日期/时间是周二 七月 17, 2018 7:57 am