C#,数值计算——插值和外推,三次样条插值(Spline_interp)的计算方法与源程序

1 文本格式

using System;

namespace Legalsoft.Truffer
{
    /// <summary>
    /// 三次样条插值
    /// Cubic Spline Interpolation
    /// Cubic spline interpolation object. Construct with x and y vectors, and
    /// (optionally) values of the first derivative at the endpoints, then call
    /// interp for interpolated values
    /// </summary>
    public class Spline_interp : Base_interp
    {
        private double[] y2 { get; set; }

        public Spline_interp(double[] xv, double[] yv, double yp1 = 1.0e99, double ypn = 1.0e99) : base(xv, yv[0], 2)
        {
            this.y2 = new double[xv.Length];
            sety2(xv, yv, yp1, ypn);
        }

        public Spline_interp(double[] xv, double yv, double yp1 = 1.0e99, double ypn = 1.0e99) : base(xv, yv, 2)
        {
            this.y2 = new double[xv.Length];
            sety2(xv, y2, yp1, ypn);
        }

        /// <summary>
        /// This routine stores an array y2[0..n - 1] with second derivatives of the
        /// interpolating function at the tabulated points pointed to by xv, using
        /// function values pointed to by yv.If yp1 and/or ypn are equal to 1.0E99 or
        /// larger, the routine is signaled to set the corresponding boundary condition
        /// for a natural spline, with zero second derivative on that boundary;
        /// otherwise, they are the values of the first derivatives at the endpoints.
        /// </summary>
        /// <param name="xv"></param>
        /// <param name="yv"></param>
        /// <param name="yp1"></param>
        /// <param name="ypn"></param>
        public void sety2(double[] xv, double[] yv, double yp1, double ypn)
        {
            double[] u = new double[n - 1];
            if (yp1 > 0.99e99)
            {
                y2[0] = u[0] = 0.0;
            }
            else
            {
                y2[0] = -0.5;
                u[0] = (3.0 / (xv[1] - xv[0])) * ((yv[1] - yv[0]) / (xv[1] - xv[0]) - yp1);
            }
            for (int i = 1; i < n - 1; i++)
            {
                double sig = (xv[i] - xv[i - 1]) / (xv[i + 1] - xv[i - 1]);
                double p = sig * y2[i - 1] + 2.0;
                y2[i] = (sig - 1.0) / p;
                u[i] = (yv[i + 1] - yv[i]) / (xv[i + 1] - xv[i]) - (yv[i] - yv[i - 1]) / (xv[i] - xv[i - 1]);
                u[i] = (6.0 * u[i] / (xv[i + 1] - xv[i - 1]) - sig * u[i - 1]) / p;
            }
            double qn;
            double un;
            if (ypn > 0.99e99)
            {
                qn = un = 0.0;
            }
            else
            {
                qn = 0.5;
                un = (3.0 / (xv[n - 1] - xv[n - 2])) * (ypn - (yv[n - 1] - yv[n - 2]) / (xv[n - 1] - xv[n - 2]));
            }
            y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
            for (int k = n - 2; k >= 0; k--)
            {
                y2[k] = y2[k] * y2[k + 1] + u[k];
            }
        }

        /// <summary>
        /// Given a value x, and using pointers to data xx and yy, this routine returns
        /// an interpolated value y, and stores an error estimate dy. The returned
        /// value is obtained by mm-point polynomial interpolation on the subrange
        /// xx[jl..jl + mm - 1].
        /// </summary>
        /// <param name="jl"></param>
        /// <param name="x"></param>
        /// <returns></returns>
        /// <exception cref="Exception"></exception>
        public override double rawinterp(int jl, double x)
        {
            int klo = jl;
            int khi = jl + 1;
            double h = xx[khi] - xx[klo];
            //if (h == 0.0)
            if (Math.Abs(h) <= float.Epsilon)
            {
                throw new Exception("Bad input to routine splint");
            }
            double a = (xx[khi] - x) / h;
            double b = (x - xx[klo]) / h;
            double y = a * yy[klo] + b * yy[khi] + ((a * a * a - a) * y2[klo] + (b * b * b - b) * y2[khi]) * (h * h) / 6.0;
            return y;
        }
    }
}
 

2 代码格式

using System;namespace Legalsoft.Truffer
{/// <summary>/// 三次样条插值/// Cubic Spline Interpolation/// Cubic spline interpolation object. Construct with x and y vectors, and/// (optionally) values of the first derivative at the endpoints, then call/// interp for interpolated values/// </summary>public class Spline_interp : Base_interp{private double[] y2 { get; set; }public Spline_interp(double[] xv, double[] yv, double yp1 = 1.0e99, double ypn = 1.0e99) : base(xv, yv[0], 2){this.y2 = new double[xv.Length];sety2(xv, yv, yp1, ypn);}public Spline_interp(double[] xv, double yv, double yp1 = 1.0e99, double ypn = 1.0e99) : base(xv, yv, 2){this.y2 = new double[xv.Length];sety2(xv, y2, yp1, ypn);}/// <summary>/// This routine stores an array y2[0..n - 1] with second derivatives of the/// interpolating function at the tabulated points pointed to by xv, using/// function values pointed to by yv.If yp1 and/or ypn are equal to 1.0E99 or/// larger, the routine is signaled to set the corresponding boundary condition/// for a natural spline, with zero second derivative on that boundary;/// otherwise, they are the values of the first derivatives at the endpoints./// </summary>/// <param name="xv"></param>/// <param name="yv"></param>/// <param name="yp1"></param>/// <param name="ypn"></param>public void sety2(double[] xv, double[] yv, double yp1, double ypn){double[] u = new double[n - 1];if (yp1 > 0.99e99){y2[0] = u[0] = 0.0;}else{y2[0] = -0.5;u[0] = (3.0 / (xv[1] - xv[0])) * ((yv[1] - yv[0]) / (xv[1] - xv[0]) - yp1);}for (int i = 1; i < n - 1; i++){double sig = (xv[i] - xv[i - 1]) / (xv[i + 1] - xv[i - 1]);double p = sig * y2[i - 1] + 2.0;y2[i] = (sig - 1.0) / p;u[i] = (yv[i + 1] - yv[i]) / (xv[i + 1] - xv[i]) - (yv[i] - yv[i - 1]) / (xv[i] - xv[i - 1]);u[i] = (6.0 * u[i] / (xv[i + 1] - xv[i - 1]) - sig * u[i - 1]) / p;}double qn;double un;if (ypn > 0.99e99){qn = un = 0.0;}else{qn = 0.5;un = (3.0 / (xv[n - 1] - xv[n - 2])) * (ypn - (yv[n - 1] - yv[n - 2]) / (xv[n - 1] - xv[n - 2]));}y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);for (int k = n - 2; k >= 0; k--){y2[k] = y2[k] * y2[k + 1] + u[k];}}/// <summary>/// Given a value x, and using pointers to data xx and yy, this routine returns/// an interpolated value y, and stores an error estimate dy. The returned/// value is obtained by mm-point polynomial interpolation on the subrange/// xx[jl..jl + mm - 1]./// </summary>/// <param name="jl"></param>/// <param name="x"></param>/// <returns></returns>/// <exception cref="Exception"></exception>public override double rawinterp(int jl, double x){int klo = jl;int khi = jl + 1;double h = xx[khi] - xx[klo];//if (h == 0.0)if (Math.Abs(h) <= float.Epsilon){throw new Exception("Bad input to routine splint");}double a = (xx[khi] - x) / h;double b = (x - xx[klo]) / h;double y = a * yy[klo] + b * yy[khi] + ((a * a * a - a) * y2[klo] + (b * b * b - b) * y2[khi]) * (h * h) / 6.0;return y;}}
}

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/192378.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

C++EasyX之井字棋

视频链接 井字棋 用EasyX和C实现井字棋小游戏 源码及注释 #include<graphics.h>char board_data[3][3] {{-,-,-},{-,-,-},{-,-,-}, };char current_piece O;//检测指定棋子的玩家是否获胜 bool CheckWin(char c) {// 检查每一行for (int i 0; i < 3; i){if (bo…

测试:性能测试

一、性能测试 性能测试是一种评估软件、系统或服务在特定条件下性能的过程。性能测试有助于确定系统的响应时间、吞吐量、可扩展性、稳定性和资源消耗等关键指标。 一、响应时间 响应时间&#xff08;Response Time&#xff09;是性能测试中的一个重要指标&#xff0c;用于衡…

算法学习系列(三):汉诺塔

目录&#xff1a; 引言一、问题描述二、问题求解三、测试四、附录&#xff08;所有代码&#xff09; 引言 这个汉诺塔问题就是一个典型的递归问题&#xff0c;这篇博客也算是上一篇的一个扩展吧&#xff0c;都是递归问题&#xff0c;这个问题太大&#xff0c;而且牵扯到的问题…

深度学习——第03章 Python程序设计语言(3.1 Python语言基础)

无论是在机器学习还是深度学习中&#xff0c;Python已经成为主导性的编程语言。而且&#xff0c;现在许多主流的深度学习框架&#xff0c;例如PyTorch、TensorFlow也都是基于Python。本课程主要是围绕“理论实战”同时进行&#xff0c;所以本章将重点介绍深度学习中Python的必备…

webGIS使用JS,高德API完成智慧校园项目打卡功能

代码实现&#xff1a; var map new AMap.Map(container,{center:[],//目标点的中心位置zoom:16,viewMode:3D,pitch:45,})//使用控件AMap.plugin([AMap.ToolBar,AMap.Scale,AMap.ControlBar,AMap.HawkEye,AMap.MoveAnimation],function(){map.addControl(new AMap.ToolBar({ pos…

【Java面试——JUC全局观、原子类、锁、集合类、线程池、工具类】

目录 3.3 JUC全局观 JUC框架包含几个部分? Lock框架和Tools哪些核心的类? JUC并发集合哪些核心的类? JUC原子类哪些核心的类? JUC线程池哪些核心的类? 3.4 JUC原子类 线程安全的实现方法有哪些? 什么是CAS? CAS使用示例&#xff0c;结合AtomicInteger给出示例? CAS会有…

Python遥感开发之快速判断TIF数据为空

Python遥感开发之快速判断TIF数据为空 前言&#xff1a;介绍一下如何使用python下的gdal读取tif数据的时候&#xff0c;快速判断该tif数据是否为空&#xff0c;如果为空的话就把当前的tif删掉。 如图所示&#xff0c;通过arcgis查看箭头指向的为空值。 仅通过文件的大小无法判…

NVMe Over Fabrics with iRDMA总结 - 1

1.0 Introduction简介 NVM Express* (NVMe*) drives are high-speed, low-latency, solid-state drives (SSDs), that connect over the server Peripheral Component Interconnect Express* (PCIe*) bus. NVM Express* (NVMe*) 硬盘是高速、低延迟的固态硬盘 (SSD),通过服务…

人工智能中的模型评估

1 概述 1.1 定义 人工智能&#xff08;AI&#xff09;模型评估是一个关键的过程&#xff0c;用于确定模型在特定任务上的性能和有效性。这个过程涉及使用各种技术和指标来衡量模型的准确度、可靠性、泛化能力以及其他重要特性。在不同的应用场景中&#xff0c;模型评估的具体…

Qt Creator 11.0.3同时使用Qt6.5和Qt5.14.2

Qt Creator 11.0.3同时使用Qt6.5和Qt5.14.2 概要方法1.打开Qt Creator中的Kit&#xff0c;这里我直接附上几张截图&#xff0c;不同的版本打开位置可能有所不同&#xff0c;总之最终目的是要打开构建套件&#xff08;Kit&#xff09;2.可以看到构建套件里面有包含了“构建套件K…

深度学习记录--计算图(前向后向传播)

什么是计算图&#xff1f; 从一个例子入手&#xff1a; 将函数J的计算用流程图表示出来&#xff0c;这样的流程图被称为计算图 简单来说&#xff0c;计算图是用来显示每个变量间的关系的一种图 两种传播方式 计算图有两种传播方式&#xff1a;前向传播 和 后向传播 什么是前…

使用dirmap命令行时报错,提示缺少gevent模块

记得以前是可以的&#xff0c;可能是时间长了重装了系统&#xff0c;引起的。 修复方法。升级pip&#xff0c;然后重新下载安装gevent模块。 具体&#xff1a; python -m pip install --upgrade pip 使用下面命令解决下载慢的问题。 pip config set global.index-url http…

Drools 7 ObjectDataCompiler

ObjectDataCompiler 是用于将一个java 对象填充到template 内生成drl。 实际是转为String对象&#xff0c;父类为DataProviderCompiler要注意对象值为空的处理&#xff0c;如果填充遇到null&#xff0c;会导致rule很奇怪。 测试java 代码如下 public class ObjectDataCompil…

【WPF.NET开发】WPF.NET桌面应用开发概述

本文内容 为何从 .NET Framework 升级使用 WPF 进行编程标记和代码隐藏输入和命令控件布局数据绑定图形和动画文本和版式自定义 WPF 应用 Windows Presentation Foundation (WPF) 是一个与分辨率无关的 UI 框架&#xff0c;使用基于矢量的呈现引擎&#xff0c;构建用于利用现…

LeetCode [中等] 二叉树中序—二叉搜索树中第K小的元素

二叉搜索树中第K小的元素 二叉搜索树具有如下性质&#xff1a; 结点的左子树只包含小于当前结点的数。 结点的右子树只包含大于当前结点的数。 所有左子树和右子树自身必须也是二叉搜索树。 二叉树的中序遍历即按照访问左子树——根结点——右子树的方式遍历二叉树&#x…

图片处理OpenCV IMDecode模式说明【生产问题处理】

OpenCV IMDecode模式说明【生产问题处理】 1 前言 今天售后同事反馈说客户使用我们的图片处理&#xff0c;将PNG图片处理为JPG图片之后&#xff0c;变为了白板。 我们图片处理使用的是openCV来进行处理 2 分析 2.1 图片是否损坏&#xff1a;非标准PNG头部 于是&#xff0c;马…

SHAP(六):使用 XGBoost 和 HyperOpt 进行信用卡欺诈检测

SHAP&#xff08;六&#xff09;&#xff1a;使用 XGBoost 和 HyperOpt 进行信用卡欺诈检测 本笔记本介绍了 XGBoost Classifier 在金融行业中的实现&#xff0c;特别是在信用卡欺诈检测方面。 构建 XGBoost 分类器后&#xff0c;它将使用 HyperOpt 库&#xff08;sklearn 的 …

【U8+】用友U8删除固定资产卡片,提示:当前卡片不是本月录入的卡片,不能删除。

【问题描述】 用友U8软件&#xff0c;参照已有账套新建账套的时候&#xff0c;选择结转期初余额。 例如&#xff1a;参照已有账套的2022年新建2023年的账套。 结转期初的时候勾选了固定资产模块&#xff0c; 建立成功后登录23年新的账套后&#xff0c;删除固定资产卡片&#xf…

基于Eclipse+SDK+ADT+DDMS的安卓开发环境完整搭建过程

基于EclipseSDKADTDDMS的安卓开发环境完整搭建过程 1 基本概念2 SDK安装3 Eclipse安装4 ADT插件安装4.1 在线安装&#xff08;太慢不建议选择&#xff09;4.2 离线安装&#xff08;建议选择&#xff09; 5 配置SDK6 集成安装7 创建安卓虚拟设备8 创建并启动安卓虚拟机8 关于DDM…

node.js出现version `GLIBC_2.27‘ not found的解决方案

大家好,我是爱编程的喵喵。双985硕士毕业,现担任全栈工程师一职,热衷于将数据思维应用到工作与生活中。从事机器学习以及相关的前后端开发工作。曾在阿里云、科大讯飞、CCF等比赛获得多次Top名次。现为CSDN博客专家、人工智能领域优质创作者。喜欢通过博客创作的方式对所学的…