C#，数值计算——函数计算，Ratfn的计算方法与源程序

1 文本格式

using System;

namespace Legalsoft.Truffer
{
    public class Ratfn
    {
        private double[] cofs { get; set; }
        private int nn { get; set; }
        private int dd { get; set; }

        public Ratfn(double[] num, double[] den)
        {
            this.cofs = new double[num.Length + den.Length - 1];
            this.nn = num.Length;
            this.dd = den.Length;

            for (int j = 0; j < nn; j++)
            {
                cofs[j] = num[j] / den[0];
            }
            for (int j = 1; j < dd; j++)
            {
                cofs[j + nn - 1] = den[j] / den[0];
            }
        }

        public Ratfn(double[] coffs, int n, int d)
        {
            this.cofs = Globals.CopyFrom(coffs);
            this.nn = n;
            this.dd = d;
        }

        public double get(double x)
        {
            double sumn = 0.0;
            double sumd = 0.0;
            for (int j = nn - 1; j >= 0; j--)
            {
                sumn = sumn * x + cofs[j];
            }
            for (int j = nn + dd - 2; j >= nn; j--)
            {
                sumd = sumd * x + cofs[j];
            }
            return sumn / (1.0 + x * sumd);
        }

        public static Ratfn pade(double[] cof)
        {
            int n = (cof.Length - 1) / 2;
            double[,] q = new double[n, n];
            double[,] qlu = new double[n, n];
            double[] x = new double[n];
            double[] y = new double[n];
            for (int j = 0; j < n; j++)
            {
                y[j] = cof[n + j + 1];
                for (int k = 0; k < n; k++)
                {
                    q[j, k] = cof[j - k + n];
                }
            }

            LUdcmp lu = new LUdcmp(q);
            lu.solve( y,  x);

            for (int j = 0; j < 4; j++)
            {
                lu.mprove(y,  x);
            }
            for (int k = 0; k < n; k++)
            {
                double sum = cof[k + 1];
                for (int j = 0; j <= k; j++)
                {
                    sum -= x[j] * cof[k - j];
                }
                y[k] = sum;
            }

            double[] num = new double[n + 1];
            double[] denom = new double[n + 1];
            num[0] = cof[0];
            denom[0] = 1.0;
            for (int j = 0; j < n; j++)
            {
                num[j + 1] = y[j];
                denom[j + 1] = -x[j];
            }
            return new Ratfn(num, denom);
        }

        public static Ratfn ratlsq(UniVarRealValueFun fn, double a, double b, int mm, int kk, ref double dev)
        {
            const int NPFAC = 8;
            const int MAXIT = 5;
            const double BIG = 1.0e99;
            const double PIO2 = 1.570796326794896619;

            int ncof = mm + kk + 1;
            int npt = NPFAC * ncof;
            double[] bb = new double[npt];
            double[] coff = new double[ncof];
            double[] ee = new double[npt];
            double[] fs = new double[npt];
            double[] wt = new double[npt];
            double[] xs = new double[npt];
            double[,] u = new double[npt, ncof];
            Ratfn ratbest = new Ratfn(coff, mm + 1, kk + 1);

            dev = BIG;
            for (int i = 0; i < npt; i++)
            {
                if (i < (npt / 2) - 1)
                {
                    double hth = PIO2 * i / (npt - 1.0);
                    xs[i] = a + (b - a) * Globals.SQR(Math.Sin(hth));
                }
                else
                {
                    double hth = PIO2 * (npt - i) / (npt - 1.0);
                    xs[i] = b - (b - a) * Globals.SQR(Math.Sin(hth));
                }
                fs[i] = fn.funk(xs[i]);
                wt[i] = 1.0;
                ee[i] = 1.0;
            }
            double e = 0.0;
            for (int it = 0; it < MAXIT; it++)
            {
                for (int i = 0; i < npt; i++)
                {
                    double power = wt[i];
                    bb[i] = power * (fs[i] + Globals.SIGN(e, ee[i]));
                    for (int j = 0; j < mm + 1; j++)
                    {
                        u[i, j] = power;
                        power *= xs[i];
                    }
                    power = -bb[i];
                    for (int j = mm + 1; j < ncof; j++)
                    {
                        power *= xs[i];
                        u[i, j] = power;
                    }
                }
                SVD svd = new SVD(u);
                svd.solve(bb,  coff);
                double devmax = 0.0;
                double sum = 0.0;
                Ratfn rat = new Ratfn(coff, mm + 1, kk + 1);
                for (int j = 0; j < npt; j++)
                {
                    ee[j] = rat.get(xs[j]) - fs[j];
                    wt[j] = Math.Abs(ee[j]);
                    sum += wt[j];
                    if (wt[j] > devmax)
                    {
                        devmax = wt[j];
                    }
                }
                e = sum / npt;
                if (devmax <= dev)
                {
                    ratbest = rat;
                    //ratbest.CopyFrom(rat);
                    dev = devmax;
                }
                //Console.Write(" ratlsq iteration= ");
                //Console.Write(it);
                //Console.Write("  max error= ");
                //Console.Write("{0,10}", devmax);
                //Console.Write("{0}", "\n");
            }
            return ratbest;
        }

    }
}
 

2 代码格式

using System;namespace Legalsoft.Truffer
{public class Ratfn{private double[] cofs { get; set; }private int nn { get; set; }private int dd { get; set; }public Ratfn(double[] num, double[] den){this.cofs = new double[num.Length + den.Length - 1];this.nn = num.Length;this.dd = den.Length;for (int j = 0; j < nn; j++){cofs[j] = num[j] / den[0];}for (int j = 1; j < dd; j++){cofs[j + nn - 1] = den[j] / den[0];}}public Ratfn(double[] coffs, int n, int d){this.cofs = Globals.CopyFrom(coffs);this.nn = n;this.dd = d;}public double get(double x){double sumn = 0.0;double sumd = 0.0;for (int j = nn - 1; j >= 0; j--){sumn = sumn * x + cofs[j];}for (int j = nn + dd - 2; j >= nn; j--){sumd = sumd * x + cofs[j];}return sumn / (1.0 + x * sumd);}public static Ratfn pade(double[] cof){int n = (cof.Length - 1) / 2;double[,] q = new double[n, n];double[,] qlu = new double[n, n];double[] x = new double[n];double[] y = new double[n];for (int j = 0; j < n; j++){y[j] = cof[n + j + 1];for (int k = 0; k < n; k++){q[j, k] = cof[j - k + n];}}LUdcmp lu = new LUdcmp(q);lu.solve( y,  x);for (int j = 0; j < 4; j++){lu.mprove(y,  x);}for (int k = 0; k < n; k++){double sum = cof[k + 1];for (int j = 0; j <= k; j++){sum -= x[j] * cof[k - j];}y[k] = sum;}double[] num = new double[n + 1];double[] denom = new double[n + 1];num[0] = cof[0];denom[0] = 1.0;for (int j = 0; j < n; j++){num[j + 1] = y[j];denom[j + 1] = -x[j];}return new Ratfn(num, denom);}public static Ratfn ratlsq(UniVarRealValueFun fn, double a, double b, int mm, int kk, ref double dev){const int NPFAC = 8;const int MAXIT = 5;const double BIG = 1.0e99;const double PIO2 = 1.570796326794896619;int ncof = mm + kk + 1;int npt = NPFAC * ncof;double[] bb = new double[npt];double[] coff = new double[ncof];double[] ee = new double[npt];double[] fs = new double[npt];double[] wt = new double[npt];double[] xs = new double[npt];double[,] u = new double[npt, ncof];Ratfn ratbest = new Ratfn(coff, mm + 1, kk + 1);dev = BIG;for (int i = 0; i < npt; i++){if (i < (npt / 2) - 1){double hth = PIO2 * i / (npt - 1.0);xs[i] = a + (b - a) * Globals.SQR(Math.Sin(hth));}else{double hth = PIO2 * (npt - i) / (npt - 1.0);xs[i] = b - (b - a) * Globals.SQR(Math.Sin(hth));}fs[i] = fn.funk(xs[i]);wt[i] = 1.0;ee[i] = 1.0;}double e = 0.0;for (int it = 0; it < MAXIT; it++){for (int i = 0; i < npt; i++){double power = wt[i];bb[i] = power * (fs[i] + Globals.SIGN(e, ee[i]));for (int j = 0; j < mm + 1; j++){u[i, j] = power;power *= xs[i];}power = -bb[i];for (int j = mm + 1; j < ncof; j++){power *= xs[i];u[i, j] = power;}}SVD svd = new SVD(u);svd.solve(bb,  coff);double devmax = 0.0;double sum = 0.0;Ratfn rat = new Ratfn(coff, mm + 1, kk + 1);for (int j = 0; j < npt; j++){ee[j] = rat.get(xs[j]) - fs[j];wt[j] = Math.Abs(ee[j]);sum += wt[j];if (wt[j] > devmax){devmax = wt[j];}}e = sum / npt;if (devmax <= dev){ratbest = rat;//ratbest.CopyFrom(rat);dev = devmax;}//Console.Write(" ratlsq iteration= ");//Console.Write(it);//Console.Write("  max error= ");//Console.Write("{0,10}", devmax);//Console.Write("{0}", "\n");}return ratbest;}}
}

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

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

相关文章

c++之xml的创建，增删改查

c++之xml的创建，增删改查

c之xml的创建&#xff0c;增删改查 1.创建写入2.添加3.删除4.修改&#xff1a; 1.创建写入 #include <stdio.h> #include <typeinfo> #include "F:/EDGE/tinyxml/tinyxml.h" #include <iostream> #include <string> #include <Winsock2.…
阅读更多...
HarmonyOS分布式文件系统开发指导

HarmonyOS分布式文件系统开发指导

分布式文件系统概述 分布式文件系统&#xff08;hmdfs&#xff0c;HarmonyOS Distributed File System&#xff09;提供跨设备的文件访问能力&#xff0c;适用于如下场景&#xff1a; 两台设备组网&#xff0c;用户可以利用一台设备上的编辑软件编辑另外一台设备上的文档。平板…
阅读更多...
实现Vue3 readonly，教你如何一步步重构

实现Vue3 readonly，教你如何一步步重构

本文通过实现readonly方法&#xff0c;一步步展示重构的流程。 前言 readonly接受一个对象&#xff0c;返回一个原值的只读代理。 实现 Vue3 中readonly方法&#xff0c;先来看一下它的使用。 <script setup> import { readonly } from "vue";let user {n…
阅读更多...
无人机航迹规划MATLAB：七种优化算法（DBO、LO、SWO、COA、LSO、KOA、GRO）求解无人机路径规划

无人机航迹规划MATLAB：七种优化算法（DBO、LO、SWO、COA、LSO、KOA、GRO）求解无人机路径规划

一、七种算法&#xff08;DBO、LO、SWO、COA、LSO、KOA、GRO&#xff09;简介 1、蜣螂优化算法DBO 蜣螂优化算法&#xff08;Dung beetle optimizer&#xff0c;DBO&#xff09;由Jiankai Xue和Bo Shen于2022年提出&#xff0c;该算法主要受蜣螂的滚球、跳舞、觅食、偷窃和繁…
阅读更多...
php 插入排序算法实现

php 插入排序算法实现

插入排序是一种简单直观的排序算法&#xff0c;它的基本思想是将一个数据序列分为有序区和无序区&#xff0c;每次从无序区选择一个元素插入到有序区的合适位置&#xff0c;直到整个序列有序为止 5, 3, 8, 2, 0, 1 HP中可以使用以下代码实现插入排序算法&#xff1a; functi…
阅读更多...
Scala---方法与函数

Scala---方法与函数

一、Scala方法的定义 有参方法&无参方法 def fun (a: Int , b: Int) : Unit {println(ab) } fun(1,1)def fun1 (a: Int , b: Int) ab println(fun1(1,2)) 注意点&#xff1a; 方法定义语法 用def来定义可以定义传入的参数&#xff0c;要指定传入参数的类型方法可以写返…
阅读更多...
西门子精彩触摸屏SMART LINE V4 面板使用U盘下载项目程序的具体方法示例

西门子精彩触摸屏SMART LINE V4 面板使用U盘下载项目程序的具体方法示例

西门子精彩触摸屏SMART LINE V4 面板使用U盘下载项目程序的具体方法示例 WinCC flexible SMART V4 SP1 软件针对SMART LINE V4 面板新增了使用U盘下载项目功能。 注意：“使用U盘下载项目”功能仅支持触摸屏OS版本为V4.0.1.0 及以上的设备。 使用U盘下载项目的步骤可参考以下内…
阅读更多...
Skybox天空盒子的更换教程_unity基础开发教程

Skybox天空盒子的更换教程_unity基础开发教程

Skybox天空盒子的更换 Skybox的下载与导入更换SkyboxSkybox属性自定义 Skybox的下载与导入 打开资源商店 搜索FREE Skybox 这里是我使用的是这一款资源&#xff0c;点击添加至我的资源 打开包管理器Package Manager Packages选择My Assets 搜索Sky 选择刚刚添加的天空盒子 点…
阅读更多...
高质量实时渲染笔记

高质量实时渲染笔记

文章目录 Real-time shadows1 自遮挡问题2 解决阴影detach问题&#xff1f;3 Aliasing4 近似积分5 percentage closer soft shadows(PCSS)percenta closer filtering(PCF)PCSS的思想 6 Variance Soft Shadow Mapping (VSSM)步骤Moment Shadow Mapping 7 Distance field shadow …
阅读更多...
Leetcode刷题详解—— 图像渲染

Leetcode刷题详解—— 图像渲染

1. 题目链接&#xff1a;733. 图像渲染 2. 题目描述&#xff1a; 有一幅以 m x n 的二维整数数组表示的图画 image &#xff0c;其中 image[i][j] 表示该图画的像素值大小。 你也被给予三个整数 sr , sc 和 newColor 。你应该从像素 image[sr][sc] 开始对图像进行 上色填充 。…
阅读更多...
PHP使用文件缓存实现html静态化

PHP使用文件缓存实现html静态化

<?php // 动态生成的内容 $content "<html><body><h1>time:".date("Y-m-d H:i:s")."</h1></body></html>"; // 静态文件保存路径和文件名 $staticFilePath "file.html"; if(file_exists($s…
阅读更多...
3D Gaussian Splatting文件的压缩【3D高斯泼溅】

3D Gaussian Splatting文件的压缩【3D高斯泼溅】

在上一篇文章中&#xff0c;我开始研究高斯泼溅&#xff08;3DGS&#xff1a;3D Gaussian Splatting&#xff09;。 它的问题之一是数据集并不小。 渲染图看起来不错。 但“自行车”、“卡车”、“花园”数据集分别是一个 1.42GB、0.59GB、1.35GB 的 PLY 文件。 它们几乎按原样…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023