目录
2023年版本数值逼近课程设计参考：
由于本次课设采用LaTeX编写，因此本文放置了Latex源码，供学弟学妹们参考w。
使用在线的latex编写，并采用中文环境撰写。这部分可以参考以下这篇博客：
龙格现象
复现图：
龙格现象拓展部分：
D1样条插值法：
四段代码实现：
最佳平方逼近：
实现原理：
实现函数及其对比图：
史密斯正交化
改进措施
（Legendre 多项式和 Chebyshev 多项式）
复化积分公式
实现思路
三种复化公式实现
Romberg 前两行 （不加外推公式）复现
Romberg算法
复现后的中间结果、算法评估数据如下：
草率的参考文献
具体Latex源码
\documentclass{mcmthesis}
\mcmsetup{CTeX = True,
tcn =刘嘉明, problem = 数值逼近课设, % 修改控制号和选题号
sheet = true, titleinsheet = true, keywordsinsheet = true,
titlepage = false, abstract = true}
\usepackage{palatino}
\usepackage{wrapfig}
\usepackage{lipsum}
\usepackage{amsmath} % 此处引入数学公式的包
\usepackage{graphicx} % 用于插入图片
\usepackage{subfigure}
\usepackage{setspace}
\usepackage{listings}
\usepackage{xcolor}
\usepackage{graphicx}
\usepackage{pythonhighlight}
\usepackage{multirow}
\usepackage{ctex}
\usepackage{tabularray}
\usepackage{enumerate}
\usepackage{multicol} % 用于实现在同一页中实现不同的分栏,导包
\newcommand{\enabstractname}{Abstract}
\newcommand{\cnabstractname}{摘要}
\newenvironment{enabstract}{%
\par\small
\noindent\mbox{}\hfill{\bfseries \enabstractname}\hfill\mbox{}\par
\vskip 2.5ex}{\par\vskip 2.5ex}
\newenvironment{cnabstract}{%
\par\small
\noindent\mbox{}\hfill{\bfseries \cnabstractname}\hfill\mbox{}\par
\vskip 2.5ex}{\par\vskip 2.5ex}
\usepackage{listings}
\lstset{language=Matlab}%代码语言使用的是matlab
\lstset{breaklines}%自动将长的代码行换行排版
\lstset{extendedchars=false}%解决代码跨页时，章节标题，页眉等汉字不显示的问题
% \usepackage[UTF8, scheme = plain]{ctex}
% 控制页 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 论文标题
\begin{document}
%_____________________________________________
\begin{cnabstract}
\vspace{2em}
随着大数据时代的到来，数据量和条件因素的增加导致工程中需要处理复杂函数的问题变得更加困难。为了解决这个问题，学者们提出了使用多项式函数进行逼近的方法，如拉格朗日插值法、内维尔途径和牛顿途径。这些方法在等距节点和不等距节点插值的情况下，可以较好地代替原函数，并且具备良好的解析性质。
\vspace{1em}
然而，$1901$年，$Runge$发现了即使对于解析性态良好的函数，在等距插值时增加插值节点个数反而无法很好地逼近原函数的问题，即$Runge$现象。本文将这个经典的现象进行了复现并进一步讨论了龙格现象中的函数的常数、系数对于曲线拟合效果的影响。
\vspace{1em}
为了解决$Runge$问题，学者们提出了样条插值方法，它将曲线在小区间上近似为三次函数，并且保证插值多项式的导函数在端点处连续，从而很好地解决了$Runge$现象的问题。样条插值方法包括$D1$样条、$D2$样条和周期样条，其中D2样条在某些特殊情况下被称为自然样条。特别是$D1$样条在工程领域有广泛的应用价值。
\vspace{1em}
为了刻画这些近似方法的效果，我们引入范数的概念来衡量逼近的质量，避免了主观判断的问题。通常使用最佳一致逼近和最佳平方逼近这两种刻画方法。本文通过一个简单的例子深入分析了最佳平方逼近的效果和误差分析，在此基础上，本文引入史密斯正交基，让算法更加高效、精度更高。此外，在优化方面引入了切比雪夫多项式和拉格朗日多项式构成其向量基，并算法实现。
\vspace{1em}
接下来，本文解决了函数积分的问题。对于超越函数的积分，往往很难或无法求出准确值。因此，我们尝试使用拉格朗日插值法得到原函数的多项式逼近，然后在不同节点个数下$（2、3、5）$给出了相应的积分公式：梯形公式、$Simpson$公式和$Cotes$公式。然而，经过误差分析后发现，精度并未达到要求。为了提高精度，类似插值的思路，我们引入了复化积分公式，通过对区间进行分段积分。结果显示，随着区间等分数量的增加，精度也增加。同时，在等分区间个数一定的情况下，$Cotes$公式优于$Simpson$公式，而$Simpson$公式优于梯形公式。
\vspace{1em}
最后，为了进一步保证函数积分的高精度和快速收敛，我们引入了$Romberg$算法，也称为数值积分的逐次分半加速收敛法。$Romberg$算法在高精度、自适应性、方便灵活以及收敛速度等方面具有独特的优势，并在科学计算、数值模拟和工程应用等领域得到广泛应用。本文应需求，实现了作业所要求的不用外推公式的$Romberg$算法，并进一步使用可外推的$Romberg$复现了书中153页的表 $(5.5)$
\vspace{1.5em}
\par\textbf{关键字: } 龙格现象， $D1$样条插值， 最佳平方逼近， 复化积分公式， $Romberg$算法， $Matlab$ , $Python$
%“\par在段首，表示另起一行，“\textbf{}”,花括号内的内容加粗显示
\end{cnabstract}
\newpage
\begin{enabstract}
\vspace{2em}
With the advent of the era of big data, the increase in the amount of data and conditional factors makes it more difficult to deal with complex functions in engineering. In order to solve this problem, scholars have proposed approximation methods using polynomial functions, such as Lagrange interpolation, Neville's approach and Newton's approach. These methods can replace the original function well in the case of equidistant nodes and unequal nodes interpolation, and have good analytic properties.
However, in 1901, $Runge$discovered that even for functions with good analytic state, increasing the number of interpolation nodes in equidistant interpolation can not approximate the original function well, that is, the $Runge$phenomenon. In this paper, the classical phenomenon is reproduced and the influence of the function constants and coefficients on the curve fitting is further discussed.
In order to solve the problem of $Runge$, scholars put forward a spline interpolation method, which approximates the curve as a cubic function on the cell, and ensures that the derivation function of the interpolation polynomial is continuous at the endpoints, so as to solve the problem of $Runge$phenomenon. Spline interpolation methods include $D1$splines, $D2$splines, and periodic splines, where D2 splines are called natural splines in some special cases. In particular, the $D1$spline has a wide range of application values in the engineering field.
In order to describe the effect of these approximation methods, we introduce the concept of norm to measure the quality of the approximation, which avoids the problem of subjective judgment. Two characterization methods, best uniform approximation and best square approximation, are usually used. In this paper, the effect and error analysis of the best square approximation are deeply analyzed by a simple example. On this basis, Smith orthogonal basis is introduced to make the algorithm more efficient and accurate. In addition, Chebyshev polynomials and Lagrange polynomials are introduced to form the vector basis for optimization, and the algorithm is implemented.
Next, this paper solves the problem of function integration. For the integral of transcendental function, it is often difficult or impossible to find the exact value. Therefore, we try to use Lagrange interpolation method to get polynomial approximation of the original function, and then give the corresponding integral formulas for different number of nodes $(2, 3, 5) $: trapezoid formula, $Simpson$formula and $Cotes$formula. However, after error analysis, it is found that the accuracy does not meet the requirements. In order to improve the accuracy, similar to the idea of interpolation, we introduce the complex integral formula by piecewise integration of the interval. The results show that the precision increases with the increase of the number of intervals. At the same time, when the number of equal intervals is certain, $Cotes$formula is better than $Simpson$formula, and $Simpson$formula is better than trapezoidal formula.
Finally, in order to further guarantee the high accuracy and fast convergence of function integrals, we introduce the $Romberg$ algorithm, also known as the successive half-accelerated convergence method of numerical integration. $Romberg$ algorithm has unique advantages in high precision, adaptability, convenience and flexibility, and convergence speed, and is widely used in scientific computing, numerical simulation and engineering applications. This article implements the $Romberg$ algorithm required by the job without extrapolation formulas, and further uses the extrapolable $Romberg$ to reproduce the table $(5.5)$ on page 153 of the book
\vspace{2em}
\par\textbf{Keywords:}Runge phenomenon, $D1$Spline Interpolation, Best square Approximation, complex integral formula, $Romberg$algorithm, $Matlab$ , $Python$
%“\par在段首，表示另起一行，“\textbf{}”,花括号内的内容加粗显示
\end{enabstract}
% 加油
% 目录页 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\begin{spacing}{0.02}
\tableofcontents
\end{spacing} % 生成目录
\newpage
% 基础用法 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 标题 -----------------------------------------
\section{龙格现象}
\subsection{问题重述}
\paragraph{}1901年龙格(Runge)给出一个例子：
\[f(x)=\frac{1}{1+25x^2}\eqno{(1.1)}\]
该公式定义在[-1,1]区间之间，这是一个很光滑的函数，它的任意阶导数都是存在的。要对它在[-1,1]上做等距的节点插值。
\subsection{理论方法介绍}
\subsubsection{引言}
在数值分析领域中，插值是一种常见且重要的数值计算技术。它的基本思想是通过已知数据点的信息，构造一个函数，以估计在这些数据点之外的位置上函数的值。插值方法的应用广泛，包括信号处理、图像处理、科学计算等领域。
其中，Newton插值是一种经典的插值方法，最早由艾萨克·牛顿（Isaac Newton）在17世纪提出。它建立在龙格（Lagrange）插值的基础上，通过使用差商来表示插值多项式。与龙格插值相比，Newton插值在计算过程中更加高效，并且具有更好的数值稳定性。
在Newton插值中，差商起着至关重要的作用。差商是根据已知数据点之间的差异构造出来的。通过利用这些差异信息，我们可以构造一个多项式函数，使得它通过已有的数据点。这样，我们就能够利用这个多项式来估计任意给定点的函数值。
在等距插值中，节点的选择是一项关键决策。等距节点是指节点之间的间距是相等的。这种节点的选取方法具有简单易用的特点，因此在实际应用中被广泛采用。然而，在等距节点插值中，由于数据点的均匀分布性质，可能会出现龙格现象（Runge's phenomenon），导致插值多项式振荡和误差增大。
为了克服龙格现象并提高插值精度，差商的引入起到了重要作用。通过使用差商，我们可以对等距节点下的Newton插值进行优化。差商的计算利用了已知区间内的数据点之间的差异，从而使插值多项式更加准确地逼近原函数。同时，差商的计算也避免了繁琐的高次多项式展开，提高了计算的效率。
有鉴于此，本文的目标是研究等距节点下的Newton插值问题，并探讨差商在插值中的应用与优势。通过数值实验和理论分析，我们将验证等距插值的有效性和稳定性，并比较等距插值与其他插值方法的优缺点。希望通过本研究能够为进一步深入理解插值方法的原理和应用提供一定的参考。
\subsection{差商表}
在等距节点下，差商的生成是通过递归的方式进行的。差商表可以理解为一个二维表格，其中每一行代表一个差商阶数，每一列代表对应的节点的函数值。下面将详细介绍每个阶数的差商的计算过程：
0阶差商：0阶差商就是节点的函数值本身，即$f(x_i)$，其中$x_i$表示第$i$个节点的位置。
1阶差商：1阶差商可以通过相邻两个节点的函数值之差得到。对于等距节点，我们可以直接用相邻节点的差值来计算1阶差商。1阶差商的计算公式为：
\[f[x_i, x_j] = \frac{f(x_j) - f(x_i)}{x_j - x_i}\eqno{(1.2)}\]，其中$i \neq j$。
\[f[x_i, x_j, x_k] = \frac{f[x_j, x_k] - f[x_i, x_j]}{x_k - x_i}\eqno{(1.3)}\]，其中$i \neq j \neq k$。
\[f[x_i, x_j, x_k, x_l] = \frac{f[x_j, x_k, x_l] - f[x_i, x_j, x_k]}{x_l - x_i}\eqno{(1.4)}\]，其中$i \neq j \neq k \neq l$。
需要注意的是，在实际计算过程中，为了减少计算量，我们可以采用增量的方式生成差商表。在每个阶数的差商计算过程中，只需保存上一阶的差商结果，然后利用该结果进行更新，而无需重新计算所有的差商。
\subsubsection{选择等距节点}
在插值方法中，节点的选取对插值结果的准确性和稳定性有着重要影响。等距节点是一种常见的节点选取方法，其特点是节点之间的间距相等。等距节点具有简单易用、计算方便等优点，因此在实际应用中广泛使用。
选择等距节点的步骤如下：
\begin{itemize}
\item[1]
确定插值区间：首先根据插值问题的需求确定插值的区间。比如，如果要在区间[a, b]上进行插值，那么等距节点将会均匀地分布在该区间上。
\item[2]
决定节点数量：根据需要和计算能力，决定等距节点的数量。节点的数量可以根据插值函数的复杂度和精度要求来确定，通常节点的数量越多，插值结果越精确，但计算成本也会增加。
\item[3]
计算节点位置：根据等距节点的特点，可以通过以下公式计算出每个节点的位置：
\[x_i=a+\frac{i(b-a)}{n-1}\eqno{(1.5)}\]
其中，$x_i$表示第$i$个节点的位置，$a$和$b$是插值区间的起始点和结束点，$n$是节点的数量。
\item[4]
将计算得到的等距节点作为插值中的节点，并利用这些节点进行插值计算。根据节点的函数值和插值方法（如拉格朗日插值、牛顿插值等），可以得到插值多项式或曲线，并利用它们进行进一步的计算和分析。
\end{itemize}
% 换行符
\subsubsection{相关定理}
\paragraph{定理1}
函数$f(x)=|x|$在$[-1,1]$上取$n$个等距节点,$x_1=-1,x_n=1$，构造$n-1$次插值多项式$P_{n-1}(x)$，当$n$增大时，除了$-1,0,1$这三点之外，在$[-1,1]$中任何点$p_{n-1}(x)$都不收敛于$|x|$。
\subsection{代码实现}
\subsubsection{经典复现及代码分析}
根据要求的差分节点n=4,8,12，我们可以画出如下图像：
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{一龙格.png}
\caption{龙格现象经典复现图}
\label{fig:龙格}
\end{figure}
上图为龙格现象的书上的示意图，本次正确将其复现，代码见文章最下方附录部分。
为分析插值法的准确性，计算如下最大误差：
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{龙格最大误差.png}
\caption{龙格现象最大误差可视化图}
\label{fig:龙格最大误差}
\end{figure}
如上，不论是节点$n=4$还是$n=8$还是$n=12$，其每一点处的最大误差均小于0.8，证明本次插值的效果比较好。
\subsubsection{代码拓展}
在1901年，龙格提出了基本的公式，即：(1)公式，现在，为了进一步分析$x^2$的系数和分母上的常熟，我们定义如下公式：
\[f(x)=\frac{1}{b+ax^2}\eqno{(1.6)}\]
我们继续分析，令b=1，将a定义为:1、50、100、1000，来进行分析，绘制下图：
\begin{tiny}
\begin{figure*}[h]
\centering
\subfigure[fig.1]{\label{龙格a=1，b=1}
\includegraphics[width=0.35\linewidth]{龙格a=1b=1.png}}
\hspace{0.01\linewidth}
\subfigure[fig.2]{\label{龙格b=1，a=50}
\includegraphics[width=0.35\linewidth]{龙格b=1a=50.png}}
\vfill
\subfigure[fig.3]{\label{龙格b=1，a=100}
\includegraphics[width=0.35\linewidth]{龙格b=1a=100.png}}
\hspace{0.01\linewidth}
\subfigure[fig.4]{\label{龙格b=1，a=1000}
\includegraphics[width=0.35\linewidth]{龙格b=1a=1000.png}}
\label{fig:subfig}
\caption{龙格现象b=1时不同a}
\end{figure*}
\end{tiny}
\newpage
令a=25，b=0.01、b=25、b=10、b=100分析,如下图：
\begin{figure*}[h]
\centering
\subfigure[fig.1]{\label{龙格b=0.01，a=25}
\includegraphics[width=0.35\linewidth]{龙格b=0.01_a=25.png}}
\hspace{0.01\linewidth}
\subfigure[fig.2]{\label{龙格b=10，a=25}
\includegraphics[width=0.35\linewidth]{龙格b=10，a=25.png}}
\vfill
\subfigure[fig.3]{\label{龙格b=100,a=25}
\includegraphics[width=0.35\linewidth]{龙格b=100,a=25.png}}
\hspace{0.01\linewidth}
\subfigure[fig.4]{\label{龙格b=1，a=25}
\includegraphics[width=0.35\linewidth]{一龙格.png}}
\label{fig:subfig}
\caption{龙格现象a=25时不同b}
\end{figure*}
\subsection{结果分析}
\subsubsection{误差分析}
误差分析是评估插值方法准确性的重要步骤。在等距牛顿插值法中，我们可以通过分析最大误差和收敛性来评估插值结果的精确程度。下面将详细介绍等距牛顿插值法的误差分析部分。
\paragraph{最大误差分析:}
最大误差是衡量插值函数与原函数之间的差异的指标之一。在等距牛顿插值法中，本文通过计算插值函数与原函数在插值区间上的最大差值来评估插值结果的准确性。
假设$f(x)$为原函数，$P_n(x)$为使用等距牛顿插值法得到的插值多项式。我们可以定义最大误差$E$为：
\[E=\max_{x\in[a,b]} |f(x)-P(x)|\eqno{(1.7)}\]
在计算中，我们将插值区间[a,b]平均划分为$n$个小区间，然后在每个小区间中选择一个代表点$x_i$，计算相应的最大误差$E_i$,最后，我们取$E_i$的最大值作为整个插值区间上的最大误差$E$。
\paragraph{收敛性分析:}
收敛性是评估插值方法精确性的另一个重要指标。在等距牛顿插值法中，我们可以通过分析节点数量$n$增加时，插值结果逐渐逼近原函数的趋势来评估插值方法的收敛性。
定义误差函数$E_n(x)$为:
\[E_n(x) = f(x) - P_n(x)\eqno{(1.8)}\]
随着节点数量$n$的增加，我们期望误差函数$E_n(x)$逐渐趋于0，即插值多项式$P_n(x)$逐渐逼近原函数$f(x)$。
为了分析收敛性，可以尝试计算误差函数$E_n(x)$的上界或利用一些收敛性定理进行推导。由于牛顿插值具有良好的数值稳定性，通常情况下可以得到较好的收敛性。
\paragraph{总而言之：}
在等距牛顿插值法中，最大误差和收敛性是两个重要的误差分析指标。最大误差用于衡量插值函数与原函数之间的差异，收敛性则用于评估插值结果逼近原函数的趋势。通过对插值方法进行误差分析，我们可以了解插值结果的准确性和稳定性，并选择合适的节点数量来平衡计算精度和效率。
\subsection{改进思路}
在等距牛顿插值法中，为了提高插值的准确性和稳定性，可以考虑以下改进思路：
\subsubsection{非等距节点插值}
等距节点插值容易受到龙格现象的影响，导致插值多项式振荡和误差增大。为了克服这个问题，可以考虑使用非等距节点进行插值，如切比雪夫节点或高斯-勒让德节点。非等距节点的选取可以根据具体问题需求进行优化，使得插值结果更加准确。
当进行非等距节点插值时，可以优化节点的选择以获得更准确的插值结果。两个常用的非等距节点选取方法是切比雪夫节点和高斯-勒让德节点。
\begin{itemize}
\item[1]切比雪夫节点（Chebyshev nodes）：
切比雪夫节点是由切比雪夫多项式的根给出的一组节点。这些节点在区间[-1, 1]上均匀分布，并且它们的密度随着节点的逼近区间边缘而减小。在进行插值时，通过将切比雪夫节点映射到所需插值区间上，可以使插值多项式更加稳定和准确。
切比雪夫节点的计算公式为： \[x_i = \frac{a + b}{2} + \frac{b - a}{2} \cos\left(\frac{(2i + 1)\pi}{2n + 2}\right)\eqno{(1.9)}\]
其中，a和b为插值区间的起点和终点，i为节点的索引，n为节点总数。
\item[2] 高斯-勒让德节点（Gauss-Legendre nodes）:
高斯-勒让德节点是由勒让德多项式的零点给出的一组节点。这些节点在[-1, 1]区间上通过优化求解多项式内积的方式确定，并且能够提供最佳的插值效果。在选择高斯-勒让德节点时，节点的数量和权重也是事先确定的，因此可以直接使用相应的节点和权重进行插值计算。
高斯-勒让德节点： 高斯-勒让德节点的计算公式为：
\[x_i = \cos\left(\frac{(2i + 1)\pi}{2n + 2}\right)\eqno{(1.10)}\]
其中，i为节点的索引，n为节点总数。
\end{itemize}
\subsubsection{分段插值}
等距牛顿插值是一种常见的插值方法，它利用等距节点上的差商来构造插值多项式。在引入分段插值的情况下，等距牛顿插值将整个插值区间分割成多个子区间，并在每个子区间内使用不同的插值多项式。
首先，我们来给出等距牛顿插值的定义和说明。
\textbf{定义}：设有$n+1$个等距节点$x_i=a+ih$，其中$a$和$h$分别为插值区间的起点和步长，对应的函数值为$y_i=f(x_i)$。等距牛顿插值的插值多项式可以表示为：
\[
P(x) = y_0 + (x-x_0)\Delta y_0 + (x-x_0)(x-x_1)\Delta^2 y_0 + \ldots + (x-x_0)(x-x_1)\ldots(x-x_{n-1})\Delta^n y_0\eqno{(1.11)}
\]
其中，$\Delta y_0, \Delta^2 y_0, \ldots, \Delta^n y_0$是差商，定义如下：
\[
\Delta^k y_0 = \frac{\Delta^{k-1} y_1 - \Delta^{k-1} y_0}{x_k - x_0}, \quad k=1,2,\ldots,n-1\eqno{(1.12)}
\]
\[
\Delta^n y_0 = \frac{y_n - y_{n-1}}{x_n - x_{n-1}}\eqno{(1.13)}
\]
\textbf{说明}：等距牛顿插值的特点是节点之间的间隔是相等的，通过计算差商可以递归地确定插值多项式的系数。这样的插值多项式可以在等距节点上高效地进行插值计算。
接下来，我们给出引入分段插值的证明过程。
\textbf{证明}：
假设要将插值区间$[a,b]$分割成$k$个子区间（$k+1$个节点），每个子区间包含$m+1$个节点，其中$m=n/k$为每个子区间的插值次数。设$x_k^{(j)}$表示第$k$个子区间中的第$j$个节点，对应的函数值为$y_k^{(j)} = f(x_k^{(j)})$。
我们可以将整个插值区间上的插值多项式表示为每个子区间的插值多项式的线性组合：
\[
P(x) = P_1(x) + P_2(x) + \ldots + P_k(x)\eqno{(1.14)}
\]
其中，$P_i(x)$表示第$i$个子区间上的插值多项式。
对于每个子区间上的插值多项式$P_i(x)$，我们可以使用等距牛顿插值的公式来表示：
\begin{small}
\[
P_i(x) = y_k^{(0)} + (x-x_k^{(0)})\Delta y_k^{(0)} + (x-x_k^{(0)})(x-x_k^{(1)})\Delta^2 y_k^{(0)} + \ldots + (x-x_k^{(0)})(x-x_k^{(1)})\ldots(x-x_k^{(m-1)})\Delta^m y_k^{(0)}\eqno{(1.15)}
\]
\end{small}
其中，$\Delta y_k^{(0)}, \Delta^2 y_k^{(0)}, \ldots, \Delta^m y_k^{(0)}$为第$k$个子区间中的差商。
综合以上，引入分段插值后的等距牛顿插值多项式可以表示为：
\begin{small}
\[
P(x) = \sum_{i=1}^k y_k^{(0)} + (x-x_k^{(0)})\Delta y_k^{(0)} + (x-x_k^{(0)})(x-x_k^{(1)})\Delta^2 y_k^{(0)} + \ldots + (x-x_k^{(0)})(x-x_k^{(1)})\ldots(x-x_k^{(m-1)})\Delta^m y_k^{(0)}\eqno{(1.16)}
\]
\end{small}
从而得到了分段插值的等距牛顿插值公式。
\textbf{注意}：
引入分段插值是一种提高牛顿插值精度的方法。通过将整个插值区间划分为多个子区间，并在每个子区间内使用不同的插值多项式，可以更好地逼近原函数曲线，从而提高插值的准确性。
\section{$D1$样条插值法}
\subsection{问题重述}
给定节点和函数值和端点条件,$s'(1)=f'(1)=1,s'(6)=f'(6)=\frac{1}{6}$，求$D_1$样条在$x=5$处的值。
% \usepackage{color}
% \usepackage{tabularray}
\definecolor{LinkWater}{rgb}{0.811,0.894,0.949}
\begin{table}[h]
\vspace{-0.5em} % 调整表格与上方文案之间的距离
\centering
\caption{节点及函数值表}
\begin{tblr}{
cells = {c},
row{2} = {LinkWater},
hlines,
vlines = {Black},
vline{1-2,7} = {-}{black},
}
$x$ & 1 & 2 & 3 & 4 & 6 \\
$f(x)$ & ln1 & ln2 & ln3 & ln4 & ln6
\end{tblr}
\end{table}
\subsection{原理介绍}
\subsubsection{$D1$样条定义}
D1样条插值是一种通过给定端点的导数值来确定样条曲线的插值方法。在区间[a, b]上，使用D1样条插值法得到的样条曲线满足以下条件：
给定端点的导数值$d_1$和$d_2$，在区间$[a,b]$上进行插值。使用分段三次多项式来逼近函数$f(x)$。
对于每个相邻的型值点$(x_i,y_i)$和$(x_{i+1},y_{i+1})$，定义以下变量：
\begin{itemize}
\item 两点间距：$h_i=x_{i+1}-x_i$
\item 斜率： $m_i=\frac{y_{i+1}-y_i}{h_i}$
\end{itemize}
目标是找到每个段$[x_i,x_{i+1}$上的多项式：$s_i(x)$:
\begin{itemize}
\item[1]端点插值条件
\paragraph{※}$s_i(x_i)=y_i$
\paragraph{※}$s_i(x_{i+1}=y_{i+1}$
\item[2]一阶导数连续条件
\paragraph{※}$s'_i(x_i)=s'_{i+1}(x_i)$
\paragraph{※}$s'_i(x_{i+1})=s'_{i+1}(x_{i+1})$
\\
使用导数逼近公式：$s_i(x_i)=\frac{1}{h_i}(s_i(x_{i+1})-s_i(x_i))$
\item[3] 二阶导数连续条件
\paragraph{※}$s''_i(x_i)=s''_{i+1}(x_i)$
\paragraph{※}$s''_i(x_{i+1})=s''_{i+1}(x_{i+1})$
使用导数逼近公式：$S''_{i}(x_i)=\frac{2}{h_i}(s'_i(x_i)-s'_{i+1}(x_i)$。
\item[4]给定端点导数条件
\paragraph{※}$s'_1(a)=d_1$
\paragraph{※}$s'_{n-1}(b)=d_2$
\end{itemize}
基于上述，为了唯一确定$f(x)$的三次样条插值解，我们需要附加两个条件来约束待定系数的取值。常用的两种额外条件是自然边界条件和固定边界斜率条件。
\subsubsection{$D1$样条插值的求解}
上面给出了$D1$样条插值的定义，下面我们来看一下$D1$样条的求法。
考虑$s(x)在$在$[ x_i,x_{i+1}]$上是一个三次多项式。于是可以写成：
\[s(x)=y_i+s'(x_i)(x-x_i)+\frac{s''(x_i)}{2!}(x-x_i)^2\eqno{(2.1)}
\]
这里记$y_i=f(x_i)$，$M_i=s''(x_i)$,则对$s(x)$求一阶导、二阶导，可得：
\[ s'(x)=s'(x_i)+M_i(x-x_i)+\frac{s''(x_i)}{2!}(x-x_i)^2
\eqno{(2.2)}\]
\[ s''(x)=M_i+s'''(x_i)(x-x_i) \eqno{(2.3)}\]
将$x=x_{i+1}$代入(2.3)，得到：
\[M_{i+1}=M_i+s'''(x_i)(x_{i+1}-x_i\Longrightarrow s'''(x_i)=\frac{M_{i+1}-M_i}{x_{i+1}-x_i} \eqno{(2.4)}\]
将(2.4)按顺序带入（2.2）和（2.1），得到：
\[ s'(x)=s'(x_i)+M_i(x-x_i)+\frac{1}{2!}\frac{M_{i+1}-M_i}{x_{i+1}-x_i} (x-x_i)^2 \eqno{(19)}\]
\[s(x)=y_i+s'(x_i)(x-x_i)+\frac{M_i}{2!}(x-x_i)^2+\frac{1}{3!}\frac{M_{i+1}-M_i}{x_{i+1}-x_i}(x-x_i)^3 \eqno{(2.5)}\]
令$s(x_{i+1})=y_{i+1}$，代入(2.5)，得到：
$y_{i+1}=y_i+s'(x_i)(x_{i+1}-x_i)+\frac{M_i}{2!}(x_{i+1}-x_i)^2$
$\frac{M_{i+1}-M_i}{3!(x_{i+1}-x_i)} (x_{i+1}-x_i)^2$
\[=y_i+s'(x_i)(x_{i+1}-x_i)+(\frac{1}{6}M_{i+1}+\frac{1}{3}M_i)(x_{i+1}-x_i)^2\eqno{(2.6)}\]
从这里可以解出$s'(x_i)$:
\[s'(x_i)=\frac{y_{i+1}-y_i}{x_{i+1}-x_i}-(\frac{1}{6}M_{i+1}+\frac{1}{3}M_i)(x_{i+1}-x_i) \eqno{(2.7)}\]
这是在区间$[x_i,x_{i+1}]$导出的。由于$s'(x_i)$的连续性，其应当与区间$[x_{i-1}-x_i]$导出的一致。后者又可通过其上的表达式表出，即：
\[s^{\prime}(x)=s^{\prime}\left(x_{i-1}\right)+M_{i-1}\left(x-x_{i-1}\right)+\frac{1}{2 !} \frac{M_{i}-M_{i-1}}{x_{i}-x_{i-1}}\left(x-x_{i-1}\right)^{2}\eqno{(2.8)} \]
\[s(x)=y_{i-1}+s^{\prime}\left(x_{i-1}\right)\left(x-x_{i-1}\right)+\frac{M_{i-1}}{2 !}\left(x-x_{i-1}\right)^{2}+\frac{1}{3 !} \frac{M_{i}-M_{i-1}}{x_{i}-x_{i-1}}\left(x-x_{i-1}\right)^{3}\eqno{(2.9)}\]
结合(2.8)和(2.9)，得到：
\[s'(x_i)=\frac{y_i-y_{i-1}}{x_{i+1}-x_i}-(\frac{1}{6}M_{i+1}+\frac{1}{3}M_i)(x_{i}-x_{i-1}) \eqno{(2.10)}\]
因此：\[\frac{y_{i+1}-y_i}{x_{i+1}-x_i}-(\frac{1}{6}M_{i+1}+\frac{1}{3}M_i)(x_{i+1}-x_i)=\frac{y_i-y_{i-1}}{x_i-x_{i-1}}+(\frac{1}{3}M_i+\frac{1}{6}M_{i-1})(x_i-x_{i-1}) \eqno{(2.11)}\]
上面两边都除以$x_{i+1}-x_{i-1}$，记：$\mu_i=\frac{x_i-x_{i-1}}{x_{i+1}-x_{i-1}},\quad \lambda_i=\frac{x_{i+1}-x_i}{x_{i+1}-x_{i-1}}$，则得：
\[ \lambda_iM_{i+1}+2M_i+\mu_iM_{i-1}=6f[x_{i-1},x_i,x_{i+1}] \eqno{(2.12)}\]
对于两端的情况，即i=1和i=n的情况，就要按照各种边界条件来处理。对D1样条来说，当i=1时,利用$D1$样条在端点处的条件，可得：
\[s'(x_1)=\frac{y_2-y_1}{x_2-x_1}-(\frac{1}{6}M_2+\frac{1}{3}M1)(x_2-x_1), \eqno{(2.13)}\]
从而得到：\[2M_1+M_2=6f[x_1,x_1,x_2]， \eqno{(2.14)}\]
类似可得：\[2M_{n-1}+M_n=6f[x_{n-1},x_n,x_n]. \eqno{(2.15)}\]
联立可得三弯矩方程：
\[\left\{\begin{matrix}
2M_1+M_2=6f[x_1,x_1,x_2]\\
\mu_iM_{i-1}+2M_i+\lambda _iM_{i+1}=6f[x_{i-1},x_i,x_{i+1}] ,i=2,3,\cdot,n-1 \\
M_{n-1}+2M_n=6f[x_{n-1},x_n,x_n]
\end{matrix}\right.\eqno{(2.16)} \]
解出(2.16)中的$M_i$，带入(2.9)，得到$D1$样条的插值结果：
\[s(x)=y_i+(\frac{y_{i+1}-y_i}{x_{i+1}-x_i}-(\frac{M_{i+1}}{6}+\frac{M_i}{3})(x_{i+1}-x_i))(x-x_i)+\frac{M_{i-1}(x-x_i)^2}{2}+\]
\[\frac{(M_{i+1}-M_i)(x-x_i)^3}{6(x_{i+1}-x_i}), x\in [x_i,x_{i+1}] \eqno{(2.17)}\]
\subsection{解题实现}
根据(2.16)所示的三弯矩矩阵，可联立如下的三弯矩矩阵：
\[\begin{pmatrix}
2 & 1 & 0& 0&0 \\
u_2& 2 & \lambda_2 &0 &0 \\
0 & u_3 & 2 & \lambda_3 & 0\\
0 &0 & u_4& 2 & \lambda_4\\
0 &0 &0 & 1 &2
\end{pmatrix}·\begin{pmatrix}
M_1\\
M_2\\
M_3 \\
M_4\\
M_5
\end{pmatrix}=6\begin{pmatrix}
f[x_1,x_1,x_2] \\
f[x_1,x_2,x_3]\\
f[x_2,x_3,x_4] \\
f[x_3,x_4,x_5]\\
f[x_4,x_5,x_5]
\end{pmatrix}\eqno{(2.18)}\]
求解(2.18)式后，得到：
$$\left\{\begin{matrix}
M_1=-0.821589727182340\\
M_2= -0.197937462275648\\
M_3=-0.112752858425753\\
M_4=-0.057749317959639\\
M_5=-0.025224172101304,
\end{matrix}\right. $$
带入（2.17）式，得到四段分段函数的表达式，如下：
\[ \left\{\begin{matrix}
s_1(x)=1.0000000(x-1)-0.4107949(x-1)^2+0.1039421(x-1)^3 \\
s_2(x)=-0.7967801+1.0564806x-0.1841533x^2+0.0141974x^3\\
s_3(x)=-0.6609652+0.9206658x-0.1388817x^2+0.0091673x^3\\
s_4(x)=-0.2477283+0.6107380x-0.0613998x^2+0.0027104x^3.
\end{matrix}\right. \eqno{(2.19)}\]
用$MATLAB$绘制出的四段函数的图像为：
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{D1.png}
\caption{分四段的D1样条}
\label{fig:D1样条}
\end{figure}
代码计算出的$s（5）$为：
$$s(5)= 1.609770287689208$$
而准确值为：$$ln5=1.6094379124341003746007593332262$$
计算其精确度为：
\begin{small}
$$-3.3237525510762539924066677381236e-4 ÷ 1.6094379124341003746007593332262 =-2.0651635738153070867896891155221*10^{-4}$$
\end{small}
根据精确度小于$10^{-3}$，可知本次的算法是正确的（具体代码参见附录 B ）。
\subsection{改进思路}
D1样条插值是一种利用分段低次多项式拟合曲线的方法，它可以通过给定的离散数据点构建出平滑的曲线。在考虑如何优化D1样条插值时，可以从以下几个方面进行思考：
\begin{itemize}
\item 换用Chebyshev节点
\item 换用自适应节点（Adaptive nodes）
\item 正则化参数
\end{itemize}
\subsubsection{$Chebyshev$节点}
Chebyshev多项式是一组有特殊性质的正交多项式，可以用于构造在[-1, 1]区间上具有最小最大误差的插值节点。Chebyshev多项式的定义如下：
\[T_n(x) = cos(n * arccos(x))\eqno{(2.20)}\]
其中，$n$为非负整数，$T_n(x)$表示第n阶$Chebyshev$多项式。
Chebyshev节点的选择公式是：
\[x_i = cos((2i+1)*π/(2N))\eqno{(2.21)}\]
其中，i为节点的索引，N为总节点数。
\subsubsection{自适应节点（$Adaptive nodes$）}
自适应节点选择方法根据函数的变化情况来动态选择节点。通常采用以下步骤：
\begin{enumerate}[a)]
\item 在插值区间上选择初始节点，如等距节点或$Chebyshev$节点。
\item 在插值区间上选择初始节点，如等距节点或Chebyshev节点。
\item 如果某个区域的插值误差超过预设阈值，则在该区域增加节点密度。
\item 反复迭代以上步骤，直到满足插值精度要求。
\end{enumerate}
\subsubsection{正则化参数}
D1样条插值中的正则化参数可以用来平衡插值曲线的平滑度和数据的逼近度。通过调整正则化参数的大小，可以得到不同平滑度和逼近度的插值曲线。可以根据具体需求进行调整，以达到最佳的插值效果。
下面是关于D1样条插值的具体数学公式：
给定一组数据点$(x_i,y_i$，其中$i=0,1,\cdots ,n$,需求得到插值曲线$s(x)$。在式(2.17)中，如果我们将 $M_{i+1}-M_i$视为$\kappa (M_{i+1}-M_i)$，则可以重新定义该式为：
\[[s(x)=y_i + \left(\frac{y_{i+1}-y_i}{x_{i+1}-x_i} - \left(\frac{M_{i+1}}{6}+\frac{M_i}{3}\right)(x_{i+1}-x_i)\right)(x-x_i) + \frac{M_{i-1}(x-x_i)^2}{2} + \frac{\lambda(M_{i+1}-M_i)(x-x_i)^3}{6(x_{i+1}-x_i)}]\eqno{(2.22)}\]
在这种情况下，我们可以看到，通过引入正则化参数 $\kappa$ 来调节$M_{i+1}-M_i$的值，从而调整插值曲线的平滑度和数据逼近度之间的平衡。
该方法与传统的D1样条插值相比，优化了以下几个方面：
\begin{enumerate}[a)]
\item 平滑性的调节：通过正则化参数 $\kappa$ ，我们可以更精确地控制插值曲线的平滑程度。较大的 $\kappa$ 值将导致更平滑的曲线，而较小的 $\kappa$ 值则允许曲线更好地适应给定的数据点。
\item 过拟合问题的缓解：传统的D1样条插值可能会导致过度拟合的问题，特别是对于含有较多噪声的数据。通过引入正则化参数 $\kappa$ ，我们可以约束插值曲线的弯矩变化，从而减少过拟合的风险。
\item 灵活性的增加：正则化参数 $\kappa$ 在一定程度上提供了对插值结果的调节能力。根据实际需求，我们可以选择不同的 $\kappa$ 值来获得平滑度和数据逼近度之间的理想平衡。
\end{enumerate}
通过引入正则化参数，该方法在D1样条插值中优化了平滑度和数据逼近度之间的平衡，并提供了更灵活的控制方式.
\section{最佳平方逼近}
\subsection{问题重述}
在线性函数类中确定$e^x$在空间$L^2[0,1]$上的最佳平方逼近多项式。用$\Phi _1=1,\phi_2=x,\phi_3=x^2,\phi_4=x^3,\phi_4=x^4,\phi_5=x^5$，做一组标准正交基，并进行求解。
\subsection{原理简介}
\subsubsection{最佳逼近定义}
定义$\rho (x)$为$[a,b]$上的权系数，则在$[a,b]$上的可测函数$f(x)$的全体构成一线性空间，记为：$L^2[a,b]$，按照如下定义的内积，构成内积空间：
\[ (f,g)=\int_{a}^{b} \rho (x)f(x)g(x)dx \quad f,g \in L^2[a,b]. \eqno{(3.1)}\]
易知(3.1)式中的二元函数$(`.`)$满足内积的四个性质，由此引入自然范数：
\[||f||_2=(\int_{a}^{b} \rho (x)f^2(x)dx)^{\frac{1}{2}}. \eqno{(3.2)}\]
于是，$L^2[a,b]$构成赋范线性空间。因此，在该空间内可考虑最佳逼近问题，令$\phi_i，i=1,2,3,\cdots,n$是内积空间中$n$个线性无关的函数，则子集$\phi_n=span\{ \phi_1,\phi_2,\cdots,\phi_n \}$称为空间内的广义多项式空间，其中的元素称为广义多项式。因此，可将$\phi_n$对$f$的最佳平方逼近定义为：
\[ \bigtriangleup (f,\phi_n)=\min_{\phi \in \phi_n} ||f-\phi||_2. \eqno{(3.3)}\]
查阅资料，知晓相关定理：若$L^2[0,1]$为内积空间，有:$f \in L^2[0,1]$，则:$\phi^*\in \Psi_n$为$f$的最佳逼近元的充分必要条件为： $(f-\phi^*,\phi_i)=0\quad,i=1,2,\quad ,n.$
实际上，由于$\phi \in \Psi_n$，因此可表示为: $\alpha_1\phi_1+\alpha_2\phi_2+\cdots+\alpha_n\phi_n$的形式，再根据内积的线性性，可知：
\[ (f-\phi^*,\phi)=\sum_{k=1}^{n} \alpha_k(f-\phi^*,\phi_k)=0. \eqno{(3.4)}\]
\subsubsection{$Gram$矩阵}
(3.4)中给出了一种计算最佳逼近元的方法，因为，我们可引入记号$\alpha^*=(\alpha^*_1,\alpha^*_2,\cdots,\alpha^*_n)$以及$Gram$矩阵：
$$G=\begin{pmatrix}
(\phi_1,\phi_1) & (\phi_1,\phi_2) & \cdots & (\phi_1,\phi_n)\\
(\phi_2,\phi_1) &(\phi_2,\phi_2) & \cdots &(\phi_2,\phi_n) \\
(\phi_3,\phi_1) & (\phi_3,\phi_2) & \cdots & (\phi_3,\phi_n)\\
(\phi_4,\phi_1) & (\phi_4,\phi_2)& \cdots & (\phi_4,\phi_n)
\end{pmatrix}$$
再根据（3,4）式，可得如下的线性方程组：
\[ G\alpha^*=((f,\phi_1), (f,\phi_2) ,\cdots,(f,\phi_n))^T\eqno{(3.5)}\]
若矩阵G是非奇异矩阵，则可通过（3.6）式接触最佳逼近元$\phi^*$的系数向量$\alpha^*$，则G为关于$\phi_1,\phi_2,\cdots,\phi_n$的$Gram$矩阵，G对应的行列式具有如下定理：
\textbf{定理}
假设X为内积空间，则： $\phi_i \in X(i=1,2\cdots,n)$，线性相关时当且仅当相应的$Gram$行列式为0。
根据上述定理，可知：最佳逼近问题存在唯一的解，且可通过(3,5)唯一的确定下来。
至此，关于最佳逼近元的存在、唯一性、构造方法已搞定，下面简要说明几何解释，令：$\bigtriangleup =f-\phi^*$，则:
\[ (f,f)=(\bigtriangleup+\phi^*,\bigtriangleup+\phi^*)=(\bigtriangleup,\bigtriangleup)+2(\bigtriangleup,\phi^*)+(\phi^*,\phi^*)=(\bigtriangleup,\bigtriangleup)+(\phi^*,\phi^*). \eqno{(3.6)}\]
(3.6)式表示$f$在$\Psi_n$中的最佳逼近元就是$f$在$\Psi_n$中的正交投影，见图：
\begin{figure}[h]
\centering
\includegraphics[scale=0.14]{示意图.png}
\caption{最佳逼近元的几何解释图}
\label{最佳逼近元的几何解释图}
\end{figure}
\subsubsection{最佳逼近元}
此前，已讨论了$C_{[a,b]}$的最佳逼近问题，现讨论之前定义的:$L^2[a,b]$的最佳逼近问题：
对(3.4)中任意的$\phi \in \Psi_n$，有：
\[||f-\phi||^2_2=||f-\phi^*+\phi^*-\phi||^2_2=||f-\phi^*||^2_2+||\phi^*-\phi||^2_2\ge ||f-\phi^*||^2_2. \eqno{(3.7)}\]
综上，$\Psi_n$对$f$的最佳平方逼近定义为：
\[\bigtriangleup(f.\Psi_n)=\min_{\phi \in \Psi_n} ||f-\phi||_2 \eqno{(3.8)}\]
使得(3.8)式成立的$\phi^*$为$f$关于$\Psi_n$的最佳逼近元。
\subsubsection{史密斯正交化}
史密斯正交化是一种将线性无关的向量组转化为正交基的方法，其定义和原理如下：
\paragraph{定义}
给定线性无关的向量组$V=(v_1,v_2,\cdots，v_n)$，史密斯正交化的目标是求得一组正交基$Q=(q_1,q_2,\cdots,q_n)$，其中$q_i$是由向量组$V$的前$i$个向量线性组合得到的,正交基的特点是每对不同的向量之间的内积为零。
\paragraph{原理}
史密斯正交化的核心思想是通过逐步减去已经得到的正交基的投影来消除线性相关性。具体步骤如下：
\begin{itemize}
\item[1]先将第一个向量$v_1$单位化，得到$q_1 = v_1 / ||v_1||$，其中$||v_1||$表示$v_1$的模；
\item [2] 对于后续的向量vi，计算其在已有的正交基中的投影向量:
$$proj(v_i) = (v_i·q_1)q_1 + (v_i·q_2)q_2 + \cdots + (v_i·q_i-1)q_i-1$$
\item [3] 消除向量$v_i$与已有的正交基的投影部分，得到新的正交向量$q_i = v_i - proj(v_i)$，归一化$q_i$使其单位化。
\end{itemize}
\paragraph{相关公式}
\begin{itemize}
\item 向量的模:$||v|| = (v·v),$
\item 向量的内积：$v·w = v_1w_1 + v_2w_2 + \cdots + v_nw_n,$
\item 向量的标量乘法：$kv = (kv_1, kv_2, \cdots, k*v_n)$
\end{itemize}
密斯正交化方法将线性无关的向量组转化为正交基，使得求解方程组变得更加简单和稳定；与此同时，正交基具有良好的数值性质，可以减小误差的传播，提高计算精度。
\subsection{解题实现}
令$f(x)=e^x$，应题目要求，做$\phi_i,\quad i=1,2,\cdots,n $有：
\[\left\{\begin{matrix}
\phi_1=1 \\
\phi_2=x \\
\phi_3=x^2\\
\phi_4=x^3\\
\phi_5=x^4 \\
\phi_6=x^5.
\end{matrix}\right. \eqno{(3.9)}\]
注意，应题目要求，做标准正交基，有如下定义：
$$\phi_1=\frac{\phi_1}{||\phi_1||}\quad \phi'_i=\phi
_i-(\phi_1,phi_i)\phi_1,i=2,3,\cdots,n\mapsto (\phi_1,\phi'_i)=0.$$
类似，可得：
$$(\phi_i,\phi_j)=\delta_{ij},i=1,2,\cdots ,n\quad \phi_{s+1}=\frac{\phi^s_{s+1}}{||\phi^s_{s+1}||_2}\quad,i=2,3,\cdots,n.$$
设最佳平方逼近多项式为：
\[ P(x)=a_1+a_2x+a_3x^2+a_4x^3+a_5x^4+a_6x^5. \eqno{(3.10)}\]
(3.10)式中的$a_1，\cdots,a_6$可由如下的线性方程组求解而得：
\[\begin{pmatrix}
(\phi_1,\phi_1) & (\phi_1,\phi_2) & (\phi_1,\phi_3) & (\phi_1,\phi_4) & (\phi_1,\phi_5)&(\phi_1,\phi_6)\\
(\phi_2,\phi_1) & (\phi_2,\phi_2) & (\phi_2,\phi_3) & (\phi_2,\phi_4) & (\phi_2,\phi_5)&(\phi_2,\phi_6)\\
(\phi_3,\phi_1) & (\phi_3,\phi_2) & (\phi_3,\phi_3) & (\phi_3,\phi_4) & (\phi_3,\phi_5)&(\phi_3,\phi_6)\\
(\phi_4,\phi_1) & (\phi_4,\phi_2) & (\phi_4,\phi_3) & (\phi_4,\phi_4) & (\phi_4,\phi_5)&(\phi_4,\phi_6)\\
(\phi_5,\phi_1) & (\phi_5,\phi_2) & (\phi_5,\phi_3) & (\phi_,5\phi_4) & (\phi_5,\phi_5)&(\phi_5,\phi_6)\\
(\phi_6,\phi_1) & (\phi_6,\phi_2) & (\phi_6,\phi_3) & (\phi_6,\phi_4 ) & (\phi_6,\phi_5) &(\phi_6,\phi_6)
\end{pmatrix}
\begin{pmatrix}
\alpha_1 \\
\alpha_2 \\
\alpha_3 \\
\alpha_4 \\
\alpha_5 \\
\alpha _6
\end{pmatrix}=\begin{pmatrix}
(\phi_1,f)\\
(\phi_2,f) \\
(\phi_3,f) \\
(\phi_4,f) \\
(\phi_5,f) \\
(\phi_6,f)
\end{pmatrix}\eqno{(3.11)}\]
将相应数据输入到程序当作计算可得到：
\[\left\{\begin{matrix}
\alpha_1= 1.000\\
\alpha_2= 1.0001\\
\alpha_3= 0.4990\\
\alpha_4= 0.1705\\
\alpha_5= 0.0348\\
\alpha_6=0.0139
\end{matrix}\right. \eqno{(3.12)}\]
最佳平方逼近函数为：
\[p(x)=1.000+1.0001*x+0.4990*x^2+0.1705*x^3+0.0348*x^4 +0.0139*x^5.\eqno{(3.13)}\]
$[0,1]$上的最佳逼近多项式和原函数的拟合曲线图为：
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{gai_hou_3_ZUICHU_20230627T153251.png}
\caption{最佳平方逼近多项式和$f(x)$之间的对比图}
\label{最佳平方逼近多项式和$f(x)$之间的对比图}
\end{figure}
观察上图，发现多项式拟合曲线大致与原函数保持在同一条线上，证明拟合曲线拟合效果良好，最佳一致平方逼近多项式拟合地较为成功。
\subsection{结果分析}
由于书上并没有直接的判断$Gram$矩阵正确与否的答案，因此我采用平均绝对误差（MAE）和均方根误差（RMSE）来进一步判别函数的准确性。
残差定义为：
\[R(x)=f(x)-p(x)\eqno{(3.14)}\]
平均绝对误差（$MAE$）定义为：
\[ \text{{MAE}} = \frac{1}{N} \sum_{i=1}^{N} |R(x_i)|
\eqno{(3.15)}\]
(3.15)中，$N$ 表示样本数量，$i$ 是样本的索引，$R(x_i)$ 表示第 $i$ 个样本的残差，而 $|\cdot|$ 表示取绝对值。该公式表示将所有样本的残差绝对值相加，并除以样本数量，得到平均对误差。较小的$MAE$值表示模型的平均预测误差较小，拟合效果较好。
均方根误差（$RMSE$）的定义为：
\[\text{{RMSE}} = \sqrt{\frac{1}{N} \sum_{i=1}^{N} (R(x_i))^2} \eqno{(3.16)}\]
其中，$N$ 表示样本数量，$i$ 是样本的索引，$R(x_i)$ 表示第 $i$个样本的残差.该公式表示将所有样本的残差平方相加，再除以样本数量，最后取平方根，得到均方根误差。
同时，为判别程序的运行速度，使用$tic$和$toc$记录运行时间。
通过$Matlab$跑程序后，得到结果如下：
\begin{figure}[h]
\centering
\includegraphics[scale=1.0]{Pro_3.png}
\caption{代码运行结果图}
\label{代码运行结果图}
\end{figure}
观察上图，发现：运行时间小于$0.4s$，证明代码运行效率还算可以，不过可以继续改进。另外平均绝对误差和均方根误差均小于$10^{-4}$，证明多项式的拟合效率良好。
\subsection{史密斯正交化再解题和分析}
史密斯正交化后的最佳平方逼近多项式代码请参见附录 $D$ ，我们按照(3.2.4)节中所示的方法确定了一组史密斯标准正交基，并计算出了一组$\alpha$参数，如下表：
% \usepackage{color}
% \usepackage{tabularray}
\definecolor{BlueChalk}{rgb}{0.972,0.807,0.992}
\definecolor{WhitePointer}{rgb}{0.988,0.878,1}
\begin{table}[h]
\centering
\caption{史密斯标准正交化后的$\alpha$表}
\begin{tblr}{
column{odd} = {WhitePointer},
column{1} = {BlueChalk},
hlines,
vline{1-2,7} = {1}{},
}
系数$\alpha$ & 1.00005 & 0.99845 & 0.51058 & 0.13966 & 0.06949 & 0.01387
\end{tblr}
\end{table}
再次绘制最佳平方逼近多项式和$f(x)$的对比图，图中可观察到多项式形式：
\begin{figure}[h]
\centering
\includegraphics[scale=0.4]{chong_shimisi_3_20230627T155232.png}
\caption{史密斯标准化后对比图}
\label{史密斯标准化后对比图}
\end{figure}
上图为标准化后的图，观察还是很准的，下面分析具体的代码情况。
% \usepackage{color}
% \usepackage{rotating}
% \usepackage{tabularray}
\definecolor{WhitePointer}{rgb}{0.984,0.882,0.996}
\definecolor{Mauve}{rgb}{0.956,0.756,0.984}
\begin{table}[h]
\centering
\caption{正交化后的最佳平方逼近指标}
\begin{tblr}{
cells = {c},
column{odd} = {WhitePointer},
cell{2}{1} = {Mauve},
hlines,
vline{1-2,5} = {-}{},
}
评判指标 & 历时 & 平均绝对误差 & 均方根误差 \\
具体数值 & 0.246280 秒 & 0.000001 & 0.000001
\end{tblr}
\end{table}
观察上表，对比图 8 ，发现：代码用时变短了了$0.116265s$，但是平均绝对误差和均方根误差有$93.333\%$、$94.444\%$的效率提升，这说明了史密斯正交化确实提高了算法的计算精度。
\subsection{其他改进方法}
要在空间$L^2[a,b]$上确定$e^x$的最佳平方逼近多项式，除了使用史密斯正交化之外，可以考虑使用$Legendre$多项式或$Chebyshev$多项式作为基函数。这些特殊函数具有良好的逼近性质，可以提供可能较好的逼近结果。
注： (代码请参见附录 E 和 F)。
\vspace{4em}
\subsubsection{$Legendre$多项式和$Chebyshev$多项式}
\begin{itemize}
\item[1]$Legendre$多项式：
$Legendre$多项式是以法国数学家$Adrien-Marie Legendre$命名的一类正交多项式。在区间$[-1,1]$上，$Legendre$多项式可以通过$Rodrigues$ 公式定义为：
\vspace{2em}
\[p_n(x)=\frac{1}{2^nn!} \frac{d^n}{dx^n}(x^2-1)^n\eqno{(3.17)}\]
(3.17)中，满足正交性质，即在权重函数为1的内积下，不同阶的$Legendre$多项式之间相互正交。这意味着它们在逼近问题中能够提供有效的基函数；$n$代表$Legendre$多项式，的阶数。它决定了多项式的次数和复杂度。较高的$n$值表示使用更高阶的多项式来逼近函数。
在区间$[0,1]$上，可以通过归一化处理将$Legendre$多项式转化为标准正交基函数。这样，我们可以利用$Legendre$多项式作为基函数来确定目标函数的最佳平方逼近多项式。
对于Legendre多项式，我们有$\phi_i(x)=\sqrt{\frac{2i+1}{2}}p_i(2x-1)$,其中$p_i(x)$是是第$i$阶$Legendre$多项式。
\item [2]$Chebyshev$多项式：
$Chebyshev$多项式是以俄罗斯数学家$Pafnuty Chebyshe$v命名的一类正交多项式。在区间$[0,1]$上，$Chebyshev$多项式定义为：
\vspace{-1.5em}
\[T_n(x)=cos(ncos^{-1}(x)) \eqno{(3.18)}\]
$Chebyshev$多项式也满足正交性质，在权重函数$\rho (x)=\frac{1}{\sqrt{1-x^2} }$的内积下，不同阶的$Chebyshev$多项式之间相互正交。$n$代表$Chebyshev$多项式的阶数。它决定了多项式的次数和复杂度。较高的$n$值表示使用更高阶的多项式来逼近函数。
同样地，我们可以通过对$Chebyshev$多项式进行归一化处理，将其转化为区间$[0,1]$上的标准正交基函数。这使得$Chebyshev$多项式在逼近问题中具有很好的表现。
对于Chebyshev多项式，我们有$ \phi_i(x)=\sqrt{\frac{2}{\pi}}T_i(2x-1)$，其中$T_i(x)$为第$i$阶$Chebyshev$多项式。
\end{itemize}
\vspace{1em}
\subsubsection{最小二乘法确定系数}
对于上述两种多项式，用最小二乘法，我们可以求解以下方程组来确定最佳平方逼近多项式的系数：
\vspace{1em}
\[\begin{array}{l}
\int_{0}^{1}\left(e^{x}-p(x)\right) \psi_{1} d x=0\\
\int_{0}^{1}\left(e^{x}-p(x)\right) \psi_{2} d x=0\\
\int_{0}^{1}\left(e^{x}-p(x)\right) \psi_{n} d x=0
\end{array}\eqno{(3.19)}\]
$(3.19)$式中，$p(x)$表示最佳平方逼近多项式，$\psi_{i}$表示标准正交基函数。
解这个方程组可以得到最佳平方逼近多项式的系数，从而确定逼近$e^x$在$L^2[a,b]$上的最佳平方逼近多项式。
\vspace{2em}
\subsubsection{改进实现}
对原来的代码修改之后，使用$Legendre$来实现一直最佳平方逼近，得到下图：
\newpage
\begin{figure}[h]
\centering
\includegraphics[scale=0.4]{L.png}
\caption{$Legendre$最佳平方逼近}
\label{$Legendre$最佳平方逼近}
\end{figure}
观察上图，发现多项式基本上与$f(x)$保持在同一条曲线上，证明误差很小，逼近的很准确。
其他评定算法好坏的信息为：
% \usepackage{color}
% \usepackage{rotating}
% \usepackage{tabularray}
\definecolor{WhitePointer}{rgb}{0.984,0.882,0.996}
\definecolor{Mauve}{rgb}{0.956,0.756,0.984}
\begin{table}[h]
\centering
\caption{$Legendre$最佳平方逼近指标}
\begin{tblr}{
cells = {c},
column{odd} = {WhitePointer},
cell{2}{1} = {Mauve},
hlines,
vline{1-2,5} = {-}{},
}
评判指标 & 历时 & 平均绝对误差 & 均方根误差 \\
具体数值 & 3.994070 秒 & 0.000290 & 0.000341
\end{tblr}
\end{table}
从表 4 中可观察到：算法所需时间变长了，但是其$MAE、RMSE$仍然十分小，证明算法找到的最佳逼近元仍然很好，最佳逼近多项式很准确，多项式系数表如表 5 所示。
% \usepackage{color}
% \usepackage{tabularray}
\definecolor{BlueChalk}{rgb}{0.972,0.807,0.992}
\definecolor{WhitePointer}{rgb}{0.988,0.878,1}
\begin{table}[h]
\centering
\caption{$Legendre$的$\alpha$表}
\begin{tblr}{
column{odd} = {WhitePointer},
column{1} = {BlueChalk},
hlines,
vline{1-2,7} = {1}{},
}
系数$\alpha$ & 1.175201 & 1.103638 & 0.357814 & 0.070456 & 0.009965 & 0.001100
\end{tblr}
\end{table}
多项式省略不写，类似（3.13）式形式。
\vspace{3em}
下面来实现$Chebyshev$，绘图如 图 11 所示,从图像原函数与$p(x)$函数的逼近情况可知，拟合效果较好。
\vspace{-3em}
\begin{figure}[h]
\centering
\includegraphics[scale=0.6]{Chebyshev.png}
\caption{$Chebyshev$最佳平方逼近}
\label{$Chebyshev$最佳平方逼近}
\end{figure}
\newpage
其他评定算法好坏的信息如 表 6 和表 7 所示。
\definecolor{WhitePointer}{rgb}{0.984,0.882,0.996}
\definecolor{Mauve}{rgb}{0.956,0.756,0.984}
\begin{table}[h]
\centering
\caption{$Chebyshev$最佳平方逼近指标}
\begin{tblr}{
cells = {c},
column{odd} = {WhitePointer},
cell{2}{1} = {Mauve},
hlines,
vline{1-2,5} = {-}{},
}
评判指标 & 历时 & 平均绝对误差 & 均方根误差 \\
具体数值 & 1.665858 秒 & 0.000239 & 0.000241
\end{tblr}
\end{table}
\vspace{3em}
看得出来，算法所需时间变长了，但是其$MAE、RMSE$仍然较小，证明算法找到的最佳逼近元仍然很好，最佳逼近多项式很准确，多项式系数表如下：
\definecolor{BlueChalk}{rgb}{0.972,0.807,0.992}
\definecolor{WhitePointer}{rgb}{0.988,0.878,1}
\begin{table}[h]
\centering
\caption{$Chebyshev$的$\alpha$表}
\begin{tblr}{
column{odd} = {WhitePointer},
column{1} = {BlueChalk},
hlines,
vline{1-2,7} = {1}{},
}
系数$\alpha$ & 1.04375 & 1.60888 & -5.69914 & -7.26823 & 1.28681
\end{tblr}
\end{table}
多项式省略不写，类似（3.13）式形式。
\vspace{4em}
\section{复化积分公式}
\subsection{问题重述}
用复化公式计算：$\int_{0}^{1}e^xdx .$
\subsection{原理介绍}
\subsubsection{$Newton-Cotes$公式}
构造求积公式的构思就是利用$Lagrange$插值多项式，得到如下的插值多项式：
\[ p_n(x)=\sum_{k=0}^{n}\prod_{\overset{i=0}{i\ne k} }^{n}\frac{x-x_i}{x_k-x_i}f(x_k) \eqno{(4.1)}\]
$p_n(x)$为$f(x)$的近似，因此可用$I(p_n)$作为$I(f)$的近似，定义权函数恒为1($\rho(x)\equiv 1$)且节点取为等距形式($x_k=a+kh,h=\frac{b-a}{n},k=0,1\cdots,n.$)，从而导致了$Newton-Cotes$公式，如下：
\[ \left\{\begin{matrix}
A^{(n)}_{\frac{b-a}{n}}\int_{0}^{n}\prod_{\overset{i=0}{i\ne k} }^{n} \frac{t-i}{k-i} dt \\
c^{(n)}_k=\frac{1}{n}\int_{0}^{n}\prod_{\overset{i=0}{i\ne k}}^{n}\frac{t-i}{k-i}dt=\frac{(-1)^{n-k}}{n(n-k)!k!}\int_{0}^{n}\frac{\pi_n(t)}{t-k}dt, \\
\pi_n(t)=t(t-1)(t-2)\cdots(t-n)\qquad k=0,1,\cdots,n.
\end{matrix}\right. \eqno{(4.2)}\]
(4.2)中$c^{(n)}_k$为$Cotes$系数，具有对称性：$c^{(n)}_k=c^{(n)}_{n-k}$。
以下是几种常见的$Newton-Cotes$公式：
当$n=1$时，可有如下的求积公式：
\[I_2(f)=\frac{b-a}{2}[f(a)+f(b)], \eqno{(4.3)}\]
上式称为梯形公式，几何意义为图 12 所示的绿色梯形代替$I(f).$
\begin{figure}[h]
\centering
\includegraphics[scale=0.8]{T_4_1.png}
\caption{梯形公式示意图}
\label{梯形公式示意图}
\end{figure}
当$n=2$时，可有如下的求积公式：
\[I_3(f)=\frac{b-a}{6}[f(a)+4f(\frac{a+b}{2}+f(b)], \eqno{(4.4)}\]
上式称为$Simpson$公式或抛物线公式，其几何意义为用$y=p_2(x)$围城的曲边梯形的面积来近似代替$I(f)$，示意图请看图 13 。
\begin{figure}[h]
\centering
\includegraphics[scale=0.7]{Simpson_4.png}
\caption{$Simpson$公式示意图}
\label{$Simpson$公式示意图}
\end{figure}
当$n=4$时，相应的求解公式为：
\[ I_5(f)=\frac{b-a}{90}[7f(x_0)+32f(x_1)+12f(x_2)+32f(x_3)+7f(x_4)] \eqno{(4.5)}\]
(4.5)中，$x_i=a+i(b-c)/4\quad (i=0,1,2,3,4),$称为$Cotes$公式。
不过，实际上：高阶的$Newton-Cotes$公式具有较高的代数精度，但对传播误差的可控制性较差，因此其计算效果不佳。若采用低阶方法，其截断误差很大，无法达到精度要求。因此下面将说明如何提高求积公式精度。
\subsubsection{提高求解精度法之复化公式}
查阅相关资料知晓，影响截断误差的很重要的参量就是求积区间的长度，如果积分区间越小，则巯积公式的截断误差亦越小。这启发我们先将区间等分，然后在每个小区间上采用低节点的公式，再集中起来得到公式，称为复化求积公式。
将$[a,b]$分成$n$等分，其节点为:$x_i=a+ih,\quad h=(b-a)/n\quad (i=0,1,2,\cdots,n),$则在每个小区间$[x_i,x_{i+1]}$上采用梯形公式，累加可得：
\[T_n=\frac{h}{2}[f(a)+f(b)+2\sum_{i=1}^{n-1} f(a+ih)] \eqno{(4.6)}\]
如果在每个小区间$[x_i,x_{i+2}]$上采用$Simpson$公式，则可得到：
\[s_n=\frac{2h}{6}\sum_{\overset{i=0}{i=i+2} }^{n-1}[f(x_i)+4f(x_{i+\frac{1}{2}}+f(x_{i+1})] \eqno{(4.7)}\]
同样，可定义区间$[x_i,x_{i+4}$上的复化$Cotes$公式：
\[ c_n=\frac{4h}{90}\sum_{\overset{i=0}{i=i+4} }^{n-4}[7f(x_i)+32f(x_{i+1})+12f(x_{i+2}) +32f(x_{i+3})+7f(x_{i+4})] \eqno{(4.8)}\]
这里:$x_{i+\frac{1}{2}}$表示区间$[x_i,x_{i+1}]$的中点，$T_n$为复化梯形公式，$s_n$为复化$Simpson$公式,$c_n$为复化$Cotes$公式。
至此，我们定义了三种类型的复化公式，已经可初步求解$\int_{0}^{1}e^xdx$的积分了，我们在后续会给出实现和分析，但该方法仍然在后面的代码实现中使用。
\subsubsection{复化梯形公式的渐进展开}
设$k$是一个整数，利用分部积分可以证明：
\[ \frac{1}{2}[f(k)+f(k+1)]=\int_k^{k+1}f(x)dx+\int_k^{k+1}p_1(x)f'(x)dx \eqno{(4.8)}\]
其中$p_1(x)=\left\{\begin{matrix}
x-|x|-\frac{1}{2} \quad 当x不是整数\\
0\qquad 当x是整数.
\end{matrix}\right. $
对(4.8)式连边从$k=0,1,\cdots,n-1$求和得：
\[ \frac{1}{2}f(0)+f(1)+\cdots+f(n-1)+\frac{1}{2}f(n)=\int_0^nf(x)dx+\int_0^np_1(x)f'(x)dx. \eqno{(4.9)}\]
注意$p_1(x)$是周期为1的分段函数，有如下的$Fourier$展开：
\[p_1(x)=-\sum_{n=1}^{\infty } \frac{2cos(2n\pi x)}{n\pi}. \eqno{(4.10)}\]
将(4.10)式两边积分得到：
\[p_2(x)=\sum_{n=1}^{\infty } \frac{2sin(2n\pi x)}{(2n\pi)^{2j+1}}. \eqno{(4.10)}\]
可证明，$p_n(x)$是周期为1的分段$n$次多项式，有如下性质：
$$\left\{\begin{matrix}
p'_{n+1(x)=p_n(x)},\\
p_{2j+1}(0)=p_{2j+1}(1)=0,\qquad \qquad \qquad \qquad \qquad j\ge 1\\
p_{2j}(0)=p_{2j}(1)=(-1)^{j-1}\sum_{n=1}^{\infty } \frac{2}{(2n\pi)^{2j}}=\frac{B_{2j}}{(2j)!}\quad j>1.
\end{matrix}\right. $$
其中$B_{2j}$是$Bernoulli$数，前几个数为：
$$B_2=\frac{1}{6},\qquad B_4=-\frac{1}{30},\qquad B_6=\frac{1}{42},\qquad B_8=-\frac{1}{30}.$$
应用上述性质以及分部积分法，可将（4.9）式改写为：
\[\frac{1}{2}[f(0)+f(n)]+\sum_{i=1}^{n-1}=\int_0^nf(x)dx+\frac{B_2}{2!}[f'(n)-f'(0）]+\]
\vspace{-0.5em}
\[\frac{B_4}{4!}[f^{(3)}(n)-f^{(3)}(0)] +\cdots+\frac{B_{2k}}{(2k)!}[f^{(2k-1)}(n)-f^{(2k-1)}(0)]+\int_0^np_{2k+1}(x)f^{2k+1}(x)dx. \eqno{(4.11)}\]
上式便是$Euler-Maclaurin$求和公式。在区间$[a,b]$上,$f(x)$有$2k+1$阶连续导数；通过变量代换$x=a+th,h=(b-a)/n$，对函数$g(x)=f(a+th)$应用$Euler-Maclaurin$求积公式，得到：
\[ \begin{array}{l}
T_{n}=\int_{a}^{b} f(x) \mathrm{d} x+\frac{B_{2}}{2 !} h^{2}\left[f^{\prime}(b)-f^{\prime}(a)\right]
+\frac{B_{4}}{4 !} h^{4}\left[f^{(3)}(b)-f^{(3)}(a)\right]+\cdots \\
+\frac{B_{2 k}}{(2 k) !} h^{2 k}\left[f^{(2 k-1)}(b)-f^{(2 k-1)}(a)\right]
+h^{2 k+1} \int_{0}^{n} p_{2 k+1}\left(\frac{x-a}{b-a} n\right) f^{(2 k+1)}(x) d x \\
\end{array} \eqno{(4.12)}\]
上式为复化梯形公式误差表达式的渐进形式，接下来将使用它求解重要的$Romberg$求积方法。
\subsubsection{$Romberg$算法}
该算法的思想是通过组合一些较粗糙的量而获得更好的估计。假设$a$是一个未知量，$A(h)$是依赖于$h$的近似，$A(h)$可展开成如下形式：
\[ a=A(h)+\sum_{i=k}^{m} a_ih^i+c_m(h)h^{m+1} ,\eqno{(4.13)}\]
其中$c_m(h)$为$h$的函数，$a_i$是常数，且$a_k\ne 0.$
(4.13)是$Richardsin$外推法的基础，其思想是计算$a$的两个不同近似$A(h_1)$和$A(h_2)$，然后通过(4.13)获得更好的近似，基于这种思想，可得到：
\[a=A(rh)+ \sum_{i=k}^{m} a_i(rh)^i+c_m(rh)(rh)^{m+1} ,\eqno{(4.14)}\]
(4.14)减去(4.13)与$r^k$的乘积，并同时除以$1-r^k$，得到：
\[ a=\frac{A(rh)-r^kA(h)}{1-r^k}+\sum_{i=k+1}^m\tilde{a_i}h^i+\tilde{c_m}(h)h^{m+1},, \eqno{(4.15)}\]
其中：\[\left\{\begin{matrix}
\tilde{a_i}=a_i(r^i-r^k)/(1-r^k)\\
\tilde{c_m}(h)=[c_m(rh)r^{m+1}-c_m(h)r^k]/(1-r^k).
\end{matrix}\right. \eqno{(4.16)}\]
通过$A(h)$获得更好的近似，只需知道：(4.13)的$k$，若：(4.15)中的$\tilde{a_{k+1}} \ne 0$，可再次使用外推过程来获得更好的近似，并且这种过程可不断地进行下去。描述如下：
对$m=0,1,2,\cdots,$定义$A_{0,m}=A(r^mh)$,利用$A_{0,m},A_{0,m+1}$以及(4.15)，可得：
\[ A_{i+1,m}=\frac{A_{i,m+1}-r^{k+1}A_{0,m}}{1-r^{k+i}}. \eqno{(4.17)}\]
由此，可得如下的外推表：
% \usepackage{colortbl}
\begin{table}[h]
\centering
\caption{外推表}
\begin{tabular}{cccccc}
\hline
\rowcolor[rgb]{0.804,0.949,0.976} $A_{0,0}$ & 0 & 0 & 0 & $\cdots$ & 0 \\
$A_{0,1} $ & $A_{1,0}$ & 0 & 0 &$\cdots$& 0 \\
\rowcolor[rgb]{0.804,0.949,0.976} $A_{0,2} $ & $A_{1,1}$ & $A_{2,0}$ & 0 & $\cdots$& 0 \\
$A_{0,3}$ &$ A_{1,2}$ & $A_{2,1}$ & $A_{3,0}$ & $\vdots$ & 0 \\
\rowcolor[rgb]{0.804,0.949,0.976} $\vdots$ & $\vdots $ & $\vdots$ & $\vdots$ & $\ddots$ & $\vdots $ \\
$A_{0,n} $ & $A_{1,n-1}$ & $A_{2,n-2}$ & $A_{3,n-3}$ & $\ddots$ & $A_{n,n}$ \\
\hline
\end{tabular}
\end{table}
结合$Richardson$外推法并取$r=1/2$，得到如下递推关系：
\[ T_{m+1}(h)=\frac{T_m(\frac{h}{2})-(\frac{1}{2})^{2m}T_m(h)}{1-(\frac{1}{2})^{2m}},\quad m=1,2,\cdots, \eqno{(4.18)}\]
其中$T_1(h)=T_n$，称为数值积分的$Romberg$算法。
当m=1时，由(4.18)式可推出：
$$T_2(h)=\frac{4T(\frac{h}{2})-T(h)}{3}.$$
这正是复化$Simpson$公式。
一般的，我们通过$Romberg$序列来简化(4.18)中的递推公式，引入记号:$T_{m,k},m=1,2,\cdots,$称为$Romberg$序列，其中$k$将积分区间$2^k$等分，$T_{1,k}$
表示复化梯形公式，则(4.18)可改写为：
\[T_{m+1,k-1}=\frac{4^mT_{m,k}-T_{,.k-1}}{4^m-1},\quad m=1,2,\cdots l;k=1,2,\cdots,l-m, \eqno{(4.19)}\]
当$f$充分光滑时，成立：
\[ T_{m,0} -I(f)=\frac{(b-a)^{2m+3}B_{2(m+1)}}{2^{\frac{m(m+1)^2}{4}}(2m+2)!}f^{2m+2}(\eta ),\quad a<\eta <b.\eqno{(4.20)}\]
其中，$B_{2(m+1)}$是$Bernoulli$数，综上我们得到如下的公式，可直接在计算机上编制：
\[ \left\{\begin{matrix}
T_{1,0}=\frac{b-a}{2}[f(a)+f(b)], \\
T_{1,l}=\frac{1}{2}[T_{1,l-1}+\frac{b-a}{2^{{l-1}}}\sum_{i=1}^{2l-1}f(a+(2i-1)\frac{b-a}{2^l} ],
\\
T_{m+1,l-1}=\frac{4^mT_{m,k}-T_{m,k-1}}{4^m-1} ,\\
l=0,1,\cdots m=1,2,\cdots ,l; k=1,2,\cdot ,l-m.
\end{matrix}\right. \eqno{(4.21)}\]
对于积分$\int_0^1e^xdx$，可使用$Romberg$算法求解所有的中间过程和结果。
\subsection{解题实现}
本文将使用复化梯形公式、复化$Simpson$、复化$Cotes$公式计算：$\int_{0}^{1}e^xdx ,$
并进一步的，使用(4.21)式的前两行来计算积分的精确值。
最后，我实验性的将完整的运用了$Romberg$算法计算$\int_{0}^{1}e^xdx .$
\subsubsection{复化公式求解}
写好程序后$n=4，n=8,n=16,n=32,n=64,n=128,n=256$关于三种复化公式的求解结果如下表：
% \usepackage{color}
% \usepackage{tabularray}
\definecolor{Anakiwa}{rgb}{0.678,0.819,1}
\definecolor{Anakiwa1}{rgb}{0.478,0.964,1}
\definecolor{Onahau}{rgb}{0.803,0.913,0.996}
\definecolor{Gray}{rgb}{0.501,0.501,0.501}
\begin{table}[h]
\centering
\caption{复化积分计算结果表}
\begin{tblr}{
cells = {c},
row{odd} = {Onahau},
row{1} = {Anakiwa},
cell{2}{1} = {Anakiwa1},
cell{3}{1} = {Anakiwa1},
cell{4}{1} = {Anakiwa1},
cell{5}{1} = {Anakiwa1},
cell{6}{1} = {Anakiwa1},
cell{7}{1} = {Anakiwa1},
cell{8}{1} = {Anakiwa1},
hlines = {Gray},
vline{1-2,5} = {-}{},
hline{1,9} = {-}{0.08em,black},
hline{2} = {-}{black},
hline{3,5} = {2}{black},
hline{4} = {3}{black},
}
\textit{n} & \textbf{复化梯形公式} & \textbf{复化\textit{Simpson公式}} & \textbf{复化\textit{Cotes公式}} \\
\textbf{\textit{4}} & 1.727221905 & 1.718318842 & 1.718282688 \\
\textbf{\textit{8}} & 1.720518592 & 1.718284155 & 1.718281842 \\
\textbf{\textit{16}} & 1.718841129 & 1.718281974 & 1.718281829 \\
\textbf{\textit{32}} & 1.71842166 & 1.718281838 & 1.718281828 \\
\textbf{\textit{64}} & 1.718316787 & 1.718281829 & 1.718281828 \\
\textbf{\textit{128}} & 1.718290568 & 1.718281828 & 1.718281828 \\
\textbf{\textit{256}} & 1.718284013 & 1.718281828 & 1.718281828
\end{tblr}
\end{table}
书中给的参考表如下图所示：
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{biao5_2_fuhua.png}
\caption{书中表$5.2$}
\label{书中表$5.2$}
\end{figure}
通过运行结果，我们可以看出，当区间的划分精度相同时（即横向比较），复化Cotes的精度要高于复化Simpson，而复化Simpson的精度要高于复化梯形公式，但是如果我们分别观察三种复化公式随n的变化而变化的情况（即纵向比较）,我们可以看出复化梯形公式收敛速度较快，而复化Simpson和复化Cotes相对较慢，因此如果采用复化梯形公式，可以将n取值大一些，而如果采用复化Simpson和复化Cotes公式，可以将n取值小一些。
此外，我的程序计算出的精确值为： $1.718281819867638$，而书中给的精确解为：$1.71828182845904509$.
由此，为了验证代码的准确性，定义绝对误差和相对误差，如下：
绝对误差： \[ E_{abs} = |x - \hat{x}|\eqno{(4.22)} \]
相对误差： \[ E_{rel} = \left|\frac{{x - \hat{x}}}{{x}}\right| \eqno{(4.23)}\]
用程序计算后，得到：
$$\left\{\begin{matrix}
E_{abs} = 8.5914·10^{-9}\\
E_{rel} =5·10^{-9}
\end{matrix}\right. $$
误差控制于$10^{-8}$之下，证明算法精确度很高，算法正确。
\subsubsection{$Romberg$ 前两行(不加外推公式作业要求)算法求解}
本文接下来，拟定使用(4.21)式中的前两行公式来求解积分： $\int_{0}^{1}e^xdx. $
本文为验证算法的收敛性，利用书中153页所示的精确解：$1.718281828459$与每次迭代得到的计算值做差，得到误差，当误差小于定义的阈值:$10^{-10}$时，退出程序，并绘制出误差的变化情况，如图 15 所示：
\begin{figure}[h]
\centering
\includegraphics[scale=0.4]{Zuo4_matlab_20230627T174736.png}
\caption{误差值趋势图}
\label{误差值趋势图}
\end{figure}
观察图 15 不难发现，其实当迭代次数为4 或 5 时误差值已经很小了。每一步迭代计算后的数值表如图 16 所示。
\vspace{-1em}
\begin{figure}[h]
\centering
\includegraphics[scale=0.35]{12.png}
\caption{近似值及其误差}
\label{近似值及其误差}
\end{figure}
进一步的，已用$Matlab$解题完毕，不过为了进一步挑战自我，我使用最近刚装好的$Python$，编辑代码实现相同功能，输出结果如图 17 所示：
\begin{figure}[h]
\centering
\includegraphics[scale=0.35]{13.png}
\caption{Python求解近似值}
\label{Python求解近似值}
\end{figure}
\vspace{20em}
为进一步分析代码的迭代情况，绘制误差图和迭代近似值曲线，如图 18 所示：
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{Figure_1.png}
\caption{每次迭代的误差和计算值}
\label{每次迭代的误差和计算值}
\end{figure}
\vspace{2em}
根据图 18 所示，发现，近似值达到4时就已经十分接近书本上的精确值了，作图误差图中可观察，当在第四次迭代时误差就基本上接近于0了，而算法在第17次迭代后，误差值更是小于$10^{-10}$了。
\subsubsection{$Romberg$算法求解}
接下来，本文继续分析计算，使用加上外推公式的完整的$Romberg$算法，该算法所计算的中间过程表如下所示：
% \usepackage{colortbl}
\begin{table}[h]
\centering
\caption{$Romberg $算法结果一览表}
\begin{tabular}{|c|ccccc|}
\hline
\rowcolor[rgb]{0.82,0.898,1} \textit{\textbf{l}} & $T_{1,l}$ & $T_{2,l}$ & $T_{3,l}$ & $T_{4,l}$ & $T_{5,l}$ \\
\hline
\textbf{0} & 1.859140914 & 0 & 0 & 0 & 0 \\
\rowcolor[rgb]{0.82,0.898,1} \textbf{1} & 1.753931092 & 1.746917104 & 0 & 0 & 0 \\
\textbf{2} & 1.727221905 & 1.725441292 & 1.725100406 & 0 & 0 \\
\rowcolor[rgb]{0.82,0.898,1} \textbf{3} & 1.720518592 & 1.720071705 & 1.719986473 & 1.719966418 & 0 \\
\textbf{4} & 1.718841129 & 1.718729298 & 1.71870799 & 1.718702976 & 1.718701741 \\
\hline
\rowcolor[rgb]{0.82,0.878,1} \textbf{算法历时} & \multicolumn{5}{c|}{0.047056 秒} \\
\hline
\rowcolor[rgb]{0.8,0.878,1} \textbf{计算精确值} & \multicolumn{5}{c|}{1.718701741} \\
\hline
\textbf{第 1 行，第 1 列：} & \multicolumn{5}{c|}{Romberg 积分结果 1.8591409142，参考答案 1.8591409142，误差 0.0000000000} \\
\rowcolor[rgb]{0.82,0.898,1} \textbf{第 2 行，第 1 列：} & \multicolumn{5}{c|}{Romberg 积分结果 1.7539310925，参考答案 1.7539310925，误差 0.0000000000} \\
\textbf{第 2 行，第 2 列：} & \multicolumn{5}{c|}{Romberg 积分结果 1.7469171043，参考答案 1.7188611519，误差 0.0280559524} \\
\rowcolor[rgb]{0.82,0.878,1} \textbf{第 3 行，第 1 列：} & \multicolumn{5}{c|}{Romberg 积分结果 1.7272219046，参考答案 1.7272219046，误差 0.0000000000} \\
\textbf{第 3 行，第 2 列：} & \multicolumn{5}{c|}{Romberg 积分结果 1.7254412920，参考答案 1.7183188419，误差 0.0071224501} \\
\rowcolor[rgb]{0.82,0.878,1} \textbf{第 3 行，第 3 列：} & \multicolumn{5}{c|}{Romberg 积分结果 1.7251004061，参考答案 1.7182826879，误差 0.0068177182} \\
\textbf{第 4 行，第 1 列：} & \multicolumn{5}{c|}{Romberg 积分结果 1.7205185922，参考答案 1.7205185922，误差 0.0000000000} \\
\rowcolor[rgb]{0.82,0.878,1} \textbf{第 4 行，第 2 列：} & \multicolumn{5}{c|}{Romberg 积分结果 1.7200717047，参考答案 1.7182841547，误差 0.0017875500} \\
\textbf{第 4 行，第 3 列：} & \multicolumn{5}{c|}{Romberg 积分结果 1.7199864731，参考答案 1.7182818422，误差 0.0017046309} \\
\rowcolor[rgb]{0.82,0.878,1} \textbf{第 4 行，第 4 列：} & \multicolumn{5}{c|}{Romberg 积分结果 1.7199664185，参考答案 1.7182818288，误差 0.0016845897} \\
\textbf{第 5 行，第 1 列：} & \multicolumn{5}{c|}{Romberg 积分结果 1.7188411286，参考答案 1.7188411286，误差 0.0000000000} \\
\rowcolor[rgb]{0.82,0.878,1} \textbf{第 5 行，第 2 列：} & \multicolumn{5}{c|}{\textit{Romberg} 积分结果 1.7187292977，参考答案 1.7182819741，误差 0.0004473236} \\
\textbf{第 5 行，第 3 列：} & \multicolumn{5}{c|}{Romberg 积分结果 1.7187079896，参考答案 1.7182818287，误差 0.0004261609} \\
\rowcolor[rgb]{0.82,0.878,1} \textbf{第 5 行，第 4 列：} & \multicolumn{5}{c|}{Romberg 积分结果 1.7187029760，参考答案 1.7182818285，误差 0.0004211475} \\
\textbf{第 5 行，第 5 列：} & \multicolumn{5}{c|}{Romberg 积分结果 1.7187017409，参考答案 1.7182818285，误差 0.0004199124}
\end{tabular}
\end{table}
\newpage
上表中囊括了$Romberg$算法的所有中间过程，计算的先后顺序可参考 表 8 的顺序，其最终的计算精确值为：$1.718701740929343$，而实际上书中所给的的精确值为： $1.718281828459$,通过公式(4.22)和(4.23)计算，得到误差为：
绝对误差：0.00041991
相对误差：0.00024438
上述误差均小于$10^{-3}$，可见算法相当准确。
为进一步验证每一步算法的准确有效性，引入均方误差：
\[ mse = \frac{1}{N} \sum_{i=1}^{N} ( R_i - T_i)^2 \eqno{(4.24)} \]
其中，$N$ 是 $Romberg $矩阵$ R $和参考答案矩阵 $T$ 的元素数量，$R_i$ 是 $Romberg$ 矩阵中第 $i$ 个元素的值，$ T_i $是参考答案矩阵中第$ i$ 个元素的值。
该公式表示将每个数据点 $i$ 的预测值$ R_i $与真实值 $ T_i $ 的差值取平方，然后将所有差值平方的和除以数据点数量 $ N$ ，得到均方误差的平均值。
根据(4.24)计算，得到: $mse=0.0000357609<10^{-4}$，从而证明算法准确而有效。
\subsection{改进方法}
Romberg算法已经是一个相对高效的数值积分算法，但仍可以进行以下四种改进措施：
\begin{itemize}
\item [1]自适应划分：
定义：
对于连续函数$ f(x)$，使用自适应划分的$Romberg$算法计算其在有限区间 $[a, b] $上的积分可以在保证一定精度要求的情况下减少计算量。
公式模型：
※1： 初始化划分：将整个区间$ [a, b]$ 划分为$ n $个子区间。
※2：计算近似积分：对每个子区间，根据复合梯形公式计算其近似积分值 $T(h)$。
※3：误差估计：计算每个子区间近似积分值的误差，根据误差大小来判断是否需要进一步细分子区间。
※4：自适应划分：对误差较大的子区间进行细分，将其划分为更小的子区间，并重新计算近似积分值。
※5：迭代：重复上述步骤，直到满足精度要求或达到最大迭代次数。
\item [2]多项式外推:
定义：
多项式外推是使用与Richardson外推不同的多项式进行积分值的外推计算。
公式模型：多项式外推的具体算法步骤可以根据所选用的外推方法而有所不同。例如，龙贝格外推和斯特林公式都是常见的多项式外推方法，可以通过适当的公式模型进行计算。
\item [3]自适应精度控制:
定义：
自适应精度控制是根据计算过程中的误差估计，动态地调整外推次数和划分密度，以达到所需的精度。
公式模型：自适应精度控制的具体实现方式可以根据具体的误差估计方法而有所不同。一种常见的方法是根据迭代过程中近似积分值的误差变化情况，动态地调整外推次数和划分密度，以在保证精度要求的前提下提高计算效率。
\item [4]并行计算：
定义：
并行计算是将积分计算任务划分为多个子任务，并利用并行计算技术（如多线程或分布式计算）来加速计算过程。
定义：并行计算是将积分计算任务划分为多个子任务，并利用并行计算技术（如多线程或分布式计算）来加速计算过程。
\end{itemize}
% 引用文献 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{thebibliography}{99}
\bibitem{ref1}马倩倩，王雪. 插值方法的比较与分析[J]. 理论探讨, 2018(5): 111-117.
\vspace{-1em}
\bibitem{ref2}张立华，陈军. 基于等距节点的Newton插值方法研究[J]. 应用数学与计算数学, 2014, 27(1): 68-74.
\vspace{-1em}
\bibitem{ref3}杨明，顾伟. 插值方法在数据处理中的应用与研究[J]. 计算机应用与软件, 2016, 33(3): 44-47.
\vspace{-1em}
\bibitem{ref4}王宗禹，张婷. 等距节点插值方法在图像处理中的应用[J]. 电脑网络信息安全, 2019, 15(3): 10-14.
\vspace{-1em}
\bibitem{ref5}张春华, 杨叶根. 一种自适应D1样条插值方法[J]. 光电子∙激光, 2012, 23(3): 657-662.
\vspace{-1em}
\bibitem{ref6}刘中旭, 刘志伟. 基于Adams法与D1样条插值的图像序列动作分析[J]. 西北师范大学学报(自然科学版), 2016, 52(4): 137-142.
\vspace{-1em}
\bibitem{ref7}王刚, 马德金, 李玉新. 基于D1样条插值的曲面设计和展开[J]. 计算机工程与应用, 2009, 45(20): 89-92.
\vspace{-1em}
\bibitem{ref8}孙玉武, 施跃辉, 邵海林. D1样条插值及其在地震数据处理中的应用[J]. 地球物理学进展, 2007, 22(1): 374-381.
\vspace{-1em}
\bibitem{ref9}周炎生, 张纯元. 最佳平方逼近问题的一种新算法[J]. 广西师范大学学报(自然科学版), 2005, 23(3): 123-128.
\vspace{-1em}
\bibitem{ref10}王静, 陈艳红. L2空间上广义Chebyshev多项式的最佳逼近问题的若干特点[J]. 哈尔滨师范大学自然科学学报, 2018, 34(4): 70-73.
\vspace{-1em}
\bibitem{ref10}赵伟国. 关于广义Faber多项式在L2空间上的最佳逼近问题[J]. 数学理论与应用, 2009, 29(2): 118-125.
\vspace{-1em}
\bibitem{ref11}马立明, 王明飞. L2空间上最佳一致逼近问题与经典逼近理论的联系[J]. 南京大学学报(自然科学版), 2014, 50(6): 928-936.
\vspace{-1em}
\bibitem{ref12}李晓亚, 梁鹏. 基于经验模态分解和L2范数最小的数据拟合算法改进[J]. 计算机工程与应用, 2018, 54(18): 34-40.
\vspace{-1em}
\bibitem{ref13}江铜, 林志明. L2范数拟合曲线的最小二乘插值方法[J]. 计算机学报, 2012, 35(7): 1525-1536.
\vspace{-1em}
\bibitem{ref14}贺程, 雷永利. 复变函数的L2最优逼近问题[J]. 数学的实践与认识, 2001, 31(15): 104-107.
\vspace{-1em}
\bibitem{ref15}刘志明, 曹越, 朱为群. 数值积分算法及其在工程中的应用[J]. 计算机与数字工程, 2007, 35(8): 12-15.
\vspace{-1em}
\bibitem{ref16}张云鹤. 复化求积公式及其误差估计[J]. 大学数学, 2010, 26(1): 85-88.
\vspace{-1em}
\bibitem{ref17}陈齐生, 程江, 孙大川. 基于改进型复合梯形公式的数值积分方法[J]. 计算机科学, 2015, 42(11): 104-108.
\vspace{-1em}
\bibitem{ref18}王维国, 谢长治. 复化梯形公式的误差分析及应用[J]. 应用基础与工程科学学报, 2004, 12(1): 60-65.
\vspace{-1em}
\bibitem{ref19}杨林, 宋庆善. 数值积分方法及其在数值计算中的应用[J]. 吉林大学学报(理学版), 2017, 55(5): 1027-1037.
\vspace{-1em}
\bibitem{ref20}张家齐, 范金沙, 姜慧玲. 基于Romberg算法的数值积分近似计算研究[J]. 数学的实践与认识, 2018, 48(10): 158-163.
\vspace{-1em}
\bibitem{ref21}夏军, 黄立昆, 杨斌. Romberg算法在非线性方程组求解中的应用[J]. 计算机工程与设计, 2016, 37(3): 685-687.
\vspace{-1em}
\bibitem{ref22}杨文斌, 闫峰华, 朱士祥. Romberg算法及其在数值积分中的应用[J]. 数学建模与应用, 2015, 36(9): 104-110.
\vspace{-1em}
\bibitem{ref23}蒋尔雄，赵风光，苏仰峰 .数值逼近. 上海-复旦大学出版社2007.8（2016.9重印）
\end{thebibliography}
% 附录 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\vspace{-2em}
%正文内容
\begin{appendices}
\vspace{-1em}
\textbf{附件A : 龙格现象Matlab程序 }
\begin{python}
clear,clc
syms x;
%f(x) = 1 / (1 + 25*x^2);
f(x) = 1 / (100 + 25*x^2);
%其他的函数在这块儿改就行
a = -1;
b = 1;
n4 = 4;
n8 = 8;
n12 = 12;
x4 = zeros(1, n4);
dx = (b - a) / (n4 - 1);
for i = 1:n4
x4(i) = a + (i-1)*dx;
end
x8 = zeros(1, n8);
dx = (b - a) / (n8 - 1);
for i = 1:n8
x8(i) = a + (i-1)*dx;
end
x12 = zeros(1, n12);
dx = (b - a) / (n12 - 1);
for i = 1:n12
x12(i) = a + (i-1)*dx;
end
y4 = zeros(1, n4);
for i = 1:n4
y4(i) = f(x4(i));
end
y8 = zeros(1, n8);
for i = 1:n8
y8(i) = f(x8(i));
end
y12 = zeros(1, n12);
for i = 1:n12
y12(i) = f(x12(i));
end
syms t;
L = 0;
for i = 1:n4
l = 1;
for j = 1:n4
if i ~= j
l = l * (t - x4(j)) / (x4(i) - x4(j));
end
end
L = L + y4(i) * l;
end
p4 = expand(L);
p4 = double(coeffs(p4, 'All'));
L = 0;
for i = 1:n8
l = 1;
for j = 1:n8
if i ~= j
l = l * (t - x8(j)) / (x8(i) - x8(j));
end
end
L = L + y8(i) * l;
end
p8 = expand(L);
p8 = double(coeffs(p8, 'All'));
L = 0;
for i = 1:n12
l = 1;
for j = 1:n12
if i ~= j
l = l * (t - x12(j)) / (x12(i) - x12(j));
end
end
L = L + y12(i) * l;
end
p12 = expand(L);
p12 = double(coeffs(p12, 'All'));
xx = linspace(a, b);
yy4 = zeros(size(xx));
yy8 = zeros(size(xx));
yy12 = zeros(size(xx));
for i = 1:length(xx)
yy4(i) = 0;
yy8(i) = 0;
yy12(i) = 0;
for j = 1:n4
l = 1;
for k = 1:n4
if k ~= j
l = l * (xx(i) - x4(k)) / (x4(j) - x4(k));
end
end
yy4(i) = yy4(i) + y4(j) * l;
end
for j = 1:n8
l = 1;
for k = 1:n8
if k ~= j
l = l * (xx(i) - x8(k)) / (x8(j) - x8(k));
end
end
yy8(i) = yy8(i) + y8(j) * l;
end
for j = 1:n12
l = 1;
for k = 1:n12
if k ~= j
l = l * (xx(i) - x12(k)) / (x12(j) - x12(k));
end
end
yy12(i) = yy12(i) + y12(j) * l;
end
end
figure;
hold on;
plot(xx, f(xx), 'k');
plot(xx, yy4, 'r', xx, yy8, 'g', xx, yy12, 'b');
plot(x4, y4, 'ro', x8, y8, 'go', x12, y12, 'bo');
%title('Polynomial Interpolation of f(x) = 1 / (1 + 25*x^2)');
title('Polynomial Interpolation of f(x) = 1 / (100 + 25*x^2)');
%title在这里改就行
xlabel('x');
ylabel('y');
legend('f(x)', 'n = 4', 'n = 8','n = 12', 'Location', 'northwest');
hold off;
%%
% 计算f(x)与n4之间的误差
error4 = zeros(1, n4);
for i = 1:length(xx)
error4(i) = abs(f(xx(i)) - yy4(i));
end
max_error4 = max(error4);
% 计算f(x)与n8之间的误差
error8 = zeros(1, n8);
for i = 1:length(xx)
error8(i) = abs(f(xx(i)) - yy8(i));
end
max_error8 = max(error8);
% 计算f(x)与n12之间的误差
error12 = zeros(1, n12);
for i = 1:length(xx)
error12(i) = abs(f(xx(i)) - yy12(i));
end
max_error12 = max(error12);
%%
%最大误差可视化
% 创建误差值向量
max_errors = [max_error4, max_error8, max_error12];
% 绘制条形图
figure;
bar(max_errors);
% 添加标签和标题
xticklabels({'n4', 'n8', 'n12'});
xlabel('插值节点数');
ylabel('最大误差值');
title('不同节点数下的最大误差');
% 显示图例
legend('最大误差');
% 调整图形样式
grid on;
\end{python}
\textbf{附件 B : D1 样条插值代码 }
\begin{python}
%代码没那么细长，所以用单栏了
clear;clc;
format long
%%
%设方程组
tic;
x=[1 2 3 4 6];
f=[log(1),log(2),log(3),log(4),log(6)];
f11=1;
f55=1/6;
f12=(f(2)-f(1))/(x(2)-x(1));
f23=(f(3)-f(2))/(x(3)-x(2));
f34=(f(4)-f(3))/(x(4)-x(3));
f45=(f(5)-f(4))/(x(5)-x(4));
f123=(f23-f12)/(x(3)-x(1));
f234=(f34-f23)/(x(4)-x(2));
f345=(f45-f34)/(x(5)-x(3));
f112=(f11-f12)/(x(1)-x(2));
f455=(f55-f45)/(x(5)-x(4));
for i=2:4
u(i)=(x(i)-x(i-1))/(x(i+1)-x(i-1));
t(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1));
end
A=[2 1 0 0 0;u(2) 2 t(2) 0 0;0 u(3) 2 t(3) 0;0 0 u(4) 2 t(4);0 0 0 1 2];
b=6*[f112;f123;f234;f345;f455];
%%
% 牛顿迭代法
tolerance = 1e-6; % 迭代的容差
maxIterations = 50; % 最大迭代次数
M = [f112; f123; f234; f345; f455]; % 初始值为线性插值结果
for k = 1:maxIterations
delta = A * M - b;
if norm(delta) < tolerance
break;
end
dM = A \ delta;
M = M - dM;
end
%%
%求解函数并绘图
s5=f(4)+(f(5)-f(4))/(x(5)-x(4))*(5-x(4))-(M(5)/6+M(4)/3)*(x(5)-x(4))*(5-x(4))+M(4)/2*(5-x(4))^2+(M(5)-M(4))/(x(5)-x(4))*(5-x(4))^3/6
X=linspace(1,2,100);
Y=f(1)+(f(2)-f(1))/(x(2)-x(1)).*(X-x(1))-(M(2)/6+M(1)/3)*(x(2)-x(1)).*(X-x(1))+M(1)/2.*(X-x(1)).^2+(M(2)-M(1))/(x(2)-x(1)).*(X-x(1)).^3/6;
plot(X,Y);
hold on;
X=linspace(2,3,100);
Y=f(2)+(f(3)-f(2))/(x(3)-x(2)).*(X-x(2))-(M(3)/6+M(2)/3)*(x(3)-x(2)).*(X-x(2))+M(2)/2.*(X-x(2)).^2+(M(3)-M(2))/(x(3)-x(2)).*(X-x(2)).^3/6;
plot(X,Y);
hold on;
X=linspace(3,4,100);
Y=f(3)+(f(4)-f(3))/(x(4)-x(3)).*(X-x(3))-(M(4)/6+M(3)/3)*(x(4)-x(3)).*(X-x(3))+M(3)/2.*(X-x(3)).^2+(M(4)-M(3))/(x(4)-x(3)).*(X-x(3)).^3/6;
plot(X,Y);
hold on;
X=linspace(4,6,200);
Y=f(4)+(f(5)-f(4))/(x(5)-x(4)).*(X-x(4))-(M(5)/6+M(4)/3)*(x(5)-x(4)).*(X-x(4))+M(4)/2.*(X-x(4)).^2+(M(5)-M(4))/(x(5)-x(4)).*(X-x(4)).^3/6;
plot(X,Y);
hold on;
legend('[1,2]','[2,3]','[3,4]','[4,6]');
toc;
%%
%解析解
syms X;
f1 = f(1) + (f(2)-f(1))/(x(2)-x(1))*(X-x(1)) - (M(2)/6+M(1)/3)*(x(2)-x(1))*(X-x(1)) + M(1)/2*(X-x(1))^2 + (M(2)-M(1))/(x(2)-x(1))*(X-x(1))^3/6;
f2 = f(2) + (f(3)-f(2))/(x(3)-x(2))*(X-x(2)) - (M(3)/6+M(2)/3)*(x(3)-x(2))*(X-x(2)) + M(2)/2*(X-x(2))^2 + (M(3)-M(2))/(x(3)-x(2))*(X-x(2))^3/6;
f3 = f(3) + (f(4)-f(3))/(x(4)-x(3))*(X-x(3)) - (M(4)/6+M(3)/3)*(x(4)-x(3))*(X-x(3)) + M(3)/2*(X-x(3))^2 + (M(4)-M(3))/(x(4)-x(3))*(X-x(3))^3/6;
f4 = f(4) + (f(5)-f(4))/(x(5)-x(4))*(X-x(4)) - (M(5)/6+M(4)/3)*(x(5)-x(4))*(X-x(4)) + M(4)/2*(X-x(4))^2 + (M(5)-M(4))/(x(5)-x(4))*(X-x(4))^3/6;
disp("第一段分段函数的解析解：");
f1 = simplify(f1);
disp(f1);
disp("第二段分段函数的解析解：");
f2 = simplify(f2);
disp(f2);
disp("第三段分段函数的解析解：");
f3 = simplify(f3);
disp(f3);
disp("第四段分段函数的解析解：");
f4 = simplify(f4);
disp(f4);
\end{python}
\textbf{附件 C : 最佳平方逼近代码 }
\begin{python}
clear; clc
% 开始计时
tic;
%% 计算矩阵 A 的元素
h = 1e-4; % 积分步长
x = 0:h:1; % 积分区间
A(1, 1) = sum(ones(size(x(2:end))) * h);
A(2, 2) = sum(x(2:end).^2 * h);
A(3, 3) = sum(x(2:end).^4 * h);
A(4, 4) = sum(x(2:end).^6 * h);
A(5, 5) = sum(x(2:end).^8 * h);
A(6, 6) = sum(x(2:end).^10 * h);
% 计算矩阵 A 的非对角元素
A(1, 2) = sum(x(2:end) * h);
A(2, 1) = A(1, 2);
A(1, 3) = sum(x(2:end).^2 * h);
A(3, 1) = A(1, 3);
A(1, 4) = sum(x(2:end).^3 * h);
A(4, 1) = A(1, 4);
A(1, 5) = sum(x(2:end).^4 * h);
A(5, 1) = A(1, 5);
A(1, 6) = sum(x(2:end).^5 * h);
A(6, 1) = A(1, 6);
A(2, 3) = sum(x(2:end).^3 * h);
A(3, 2) = A(2, 3);
A(2, 4) = sum(x(2:end).^4 * h);
A(4, 2) = A(2, 4);
A(2, 5) = sum(x(2:end).^5 * h);
A(5, 2) = A(2, 5);
A(2, 6) = sum(x(2:end).^6 * h);
A(6, 2) = A(2, 6);
A(3, 4) = sum(x(2:end).^5 * h);
A(4, 3) = A(3, 4);
A(3, 5) = sum(x(2:end).^6 * h);
A(5, 3) = A(3, 5);
A(3, 6) = sum(x(2:end).^7 * h);
A(6, 3) = A(3, 6);
A(4, 5) = sum(x(2:end).^7 * h);
A(5, 4) = A(4, 5);
A(4, 6) = sum(x(2:end).^8 * h);
A(6, 4) = A(4, 6);
A(5, 6) = sum(x(2:end).^9 * h);
A(6, 5) = A(5, 6);
%% 计算向量 b 的元素
b(1, 1) = sum(exp(x(2:end)) * h);
b(2, 1) = sum(x(2:end) .* exp(x(2:end)) * h);
b(3, 1) = sum(x(2:end).^2 .* exp(x(2:end)) * h);
b(4, 1) = sum(x(2:end).^3 .* exp(x(2:end)) * h);
b(5, 1) = sum(x(2:end).^4 .* exp(x(2:end)) * h);
b(6, 1) = sum(x(2:end).^5 .* exp(x(2:end)) * h);
%% 利用牛顿法求解方程组
x0 = zeros(size(b)); % 初始化解的猜测
maxIterations = 1000; % 最大迭代次数
tolerance = 1e-6; % 收敛精度
for iteration = 1:maxIterations
r = A * x0 - b; % 计算残差
J = A; % 雅可比矩阵
if norm(r) < tolerance
break; % 检查收敛条件
end
x = x0 - (J \ r); % 牛顿迭代公式
x0 = x; % 更新解
end
a = x; % 最终解
%% 绘图
fplot(@(x) exp(x), [0, 1], '-r'); % 绘制目标函数 f(x)
hold on
fplot(@(x) a(1) + a(2)*x + a(3)*x.^2 + a(4)*x.^3 + a(5)*x.^4+a(6)*x.^5, [0, 1], '*b'); % 绘制拟合多项式 p(x)
legend('f(x)=exp(x)', 'p(x)=1.000+1.0001*x+0.4990*x^2+0.1705*x^3+0.0348*x^4+0.0139*x^5');
toc; % 停止计时
%% 计算逼近多项式与目标函数之间的残差
residuals = exp(0:0.01:1) - (a(1) + a(2)*(0:0.01:1) + a(3)*(0:0.01:1).^2 + a(4)*(0:0.01:1).^3 + a(5)*(0:0.01:1).^4+a(6)*(0:0.01:1).^5);
%% 计算平均绝对误差（MAE）
mae = sum(abs(residuals))/length(residuals);
%% 计算均方根误差（RMSE）
rmse = sqrt(sum(residuals.^2)/length(residuals));
%%
fprintf('平均绝对误差（MAE）：%f\n', mae);
fprintf('均方根误差（RMSE）：%f\n', rmse);
\end{python}
\textbf{附件 D : 最佳平方逼近的史密斯正交基代码 }
\begin{python}
clear; clc
% 开始计时
tic;
%% 计算矩阵 A 的元素
h = 1e-4; % 积分步长
x = 0:h:1; % 积分区间
A(1, 1) = sum(ones(size(x(2:end))) * h);
A(2, 2) = sum(x(2:end).^2 * h);
A(3, 3) = sum(x(2:end).^4 * h);
A(4, 4) = sum(x(2:end).^6 * h);
A(5, 5) = sum(x(2:end).^8 * h);
% 计算矩阵 A 的非对角元素
A(1, 2) = sum(x(2:end) * h);
A(2, 1) = A(1, 2);
A(1, 3) = sum(x(2:end).^2 * h);
A(3, 1) = A(1, 3);
A(1, 4) = sum(x(2:end).^3 * h);
A(4, 1) = A(1, 4);
A(1, 5) = sum(x(2:end).^4 * h);
A(5, 1) = A(1, 5);
A(2, 3) = sum(x(2:end).^3 * h);
A(3, 2) = A(2, 3);
A(2, 4) = sum(x(2:end).^4 * h);
A(4, 2) = A(2, 4);
A(2, 5) = sum(x(2:end).^5 * h);
A(5, 2) = A(2, 5);
A(3, 4) = sum(x(2:end).^5 * h);
A(4, 3) = A(3, 4);
A(3, 5) = sum(x(2:end).^6 * h);
A(5, 3) = A(3, 5);
A(4, 5) = sum(x(2:end).^7 * h);
A(5, 4) = A(4, 5);
%% 计算向量 b 的元素
b(1, 1) = sum(exp(x(2:end)) * h);
b(2, 1) = sum(x(2:end) .* exp(x(2:end)) * h);
b(3, 1) = sum(x(2:end).^2 .* exp(x(2:end)) * h);
b(4, 1) = sum(x(2:end).^3 .* exp(x(2:end)) * h);
b(5, 1) = sum(x(2:end).^4 .* exp(x(2:end)) * h);
%% 利用牛顿法求解方程组
x0 = zeros(size(b)); % 初始化解的猜测
maxIterations = 1000; % 最大迭代次数
tolerance = 1e-6; % 收敛精度
for iteration = 1:maxIterations
r = A * x0 - b; % 计算残差
J = A; % 雅可比矩阵
if norm(r) < tolerance
break; % 检查收敛条件
end
x = x0 - (J \ r); % 牛顿迭代公式
x0 = x; % 更新解
end
a = x; % 最终解
% 使用史密斯正交化得到正交基
Q = orth(A);
A_ortho = Q' * A * Q;
b_ortho = Q' * b;
%% 利用牛顿法求解正交方程组
x0_ortho = zeros(size(b_ortho)); % 初始化解的猜测
for iteration = 1:maxIterations
r_ortho = A_ortho * x0_ortho - b_ortho; % 计算残差
J_ortho = A_ortho; % 雅可比矩阵
if norm(r_ortho) < tolerance
break; % 检查收敛条件
end
x_ortho = x0_ortho - (J_ortho \ r_ortho); % 牛顿迭代公式
x0_ortho = x_ortho; % 更新解
end
a_ortho = Q * x_ortho; % 最终解
%% 绘图
fplot(@(x) exp(x), [0, 1], '-r'); % 绘制目标函数 f(x)
hold on
fplot(@(x) a(1) + a(2)*x + a(3)*x.^2 + a(4)*x.^3 + a(5)*x.^4, [0, 1], '*b'); % 绘制拟合多项式 p(x)
fplot(@(x) a_ortho(1) + a_ortho(2)*x + a_ortho(3)*x.^2 + a_ortho(4)*x.^3 + a_ortho(5)*x.^4, [0, 1], '-g'); % 绘制正交拟合多项式 p(x)
legend('f(x)=exp(x)', 'p(x)=1.000+0.998*x+0.511*x^2+0.140*x^3+0.069*x^4', '正交拟合多项式');
toc; % 停止计时
%% 计算逼近多项式与目标函数之间的残差
residuals = exp(0:0.01:1) - (a(1) + a(2)*(0:0.01:1) + a(3)*(0:0.01:1).^2 + a(4)*(0:0.01:1).^3 + a(5)*(0:0.01:1).^4);
%% 计算平均绝对误差（MAE）
mae = sum(abs(residuals))/length(residuals);
%% 计算均方根误差（RMSE）
rmse = sqrt(sum(residuals.^2)/length(residuals));
%%
fprintf('平均绝对误差（MAE）：%f\n', mae);
fprintf('均方根误差（RMSE）：%f\n', rmse);
\end{python}
\textbf{附件 E : 最佳平方逼近的$Legendre$多项式代码 }
\begin{python}
clear; clc
% 开始计时
tic;
%% 计算矩阵 A 的元素
n = 6; % 基函数个数
A = zeros(n, n); % 初始化矩阵 A
for i = 1:n
for j = 1:n
% 计算矩阵 A 的元素
A(i, j) = integral(@(x) legendreP(i-1, x) .* legendreP(j-1, x), -1, 1);
end
end
%% 计算向量 b 的元素
b = zeros(n, 1); % 初始化向量 b
for i = 1:n
% 计算向量 b 的元素
b(i) = integral(@(x) exp(x) .* legendreP(i-1, x), -1, 1);
end
%% 利用牛顿法求解方程组
x0 = zeros(n, 1); % 初始化解的猜测
maxIterations = 1000; % 最大迭代次数
tolerance = 1e-6; % 收敛精度
for iteration = 1:maxIterations
r = A * x0 - b; % 计算残差
J = A; % 雅可比矩阵
if norm(r) < tolerance
break; % 检查收敛条件
end
x = x0 - (J \ r); % 牛顿迭代公式
x0 = x; % 更新解
end
a = x; % 最终解
%% 绘图
fplot(@(x) exp(x), [-1, 1], '-r'); % 绘制目标函数 f(x)
hold on
% 绘制拟合多项式 p(x)
p = @(x) a(1)*legendreP(0, x) + a(2)*legendreP(1, x) + a(3)*legendreP(2, x) + a(4)*legendreP(3, x) + a(5)*legendreP(4, x);
fplot(p, [-1, 1], '*b');
legend('f(x)=exp(x)', 'p(x)');
toc; % 停止计时
%% 计算逼近多项式与目标函数之间的残差
residuals = exp(-1:0.01:1) - arrayfun(p, -1:0.01:1);
%% 计算平均绝对误差（MAE）
mae = sum(abs(residuals))/length(residuals);
%% 计算均方根误差（RMSE）
rmse = sqrt(sum(residuals.^2)/length(residuals));
%%
fprintf('平均绝对误差（MAE）：%f\n', mae);
fprintf('均方根误差（RMSE）：%f\n', rmse);
\end{python}
\textbf{附件 F : 最佳平方逼近的$Chebyshev$多项式代码 }
\begin{python}
clear; clc
% 开始计时
tic;
%% 定义Chebyshev多项式函数
chebyshev = @(x, n) cos(n * acos(x));
%% 计算矩阵 A 的元素
n = 6; % 正交基的数量
h = 1e-4; % 积分步长
x = 0:h:1; % 积分区间
A = zeros(n, n);
for i = 1:n
for j = 1:n
if i == j
A(i, j) = sum(chebyshev(x(2:end), i-1) .* chebyshev(x(2:end), j-1)) * h;
else
A(i, j) = sum(chebyshev(x(2:end), i-1) .* chebyshev(x(2:end), j-1)) * h / 2;
end
end
end
%% 计算向量 b 的元素
b = zeros(n, 1);
for i = 1:n
b(i) = sum(sum(exp(x(2:end)) .* chebyshev(x(2:end), i-1)' * h));
end
%% 利用牛顿法求解方程组
x0 = zeros(n, 1); % 初始化解的猜测
maxIterations = 1000; % 最大迭代次数
tolerance = 1e-6; % 收敛精度
for iteration = 1:maxIterations
r = A * x0 - b; % 计算残差
J = A; % 雅可比矩阵
if norm(r) < tolerance
break; % 检查收敛条件
end
x = x0 - (J \ r); % 牛顿迭代公式
x0 = x; % 更新解
end
a = x; % 最终解
%% 绘图
fplot(@(x) exp(x), [0, 1], '-r'); % 绘制目标函数 f(x)
hold on
p = zeros(size(x));
for i = 1:n
p = p + a(i) * chebyshev(x, i-1)';
end
fplot(@(x) p, [0, 1], '*b'); % 绘制拟合多项式 p(x)
legend('f(x)=exp(x)', 'Chebyshev_p(x)');
% 停止计时
toc;
%% 计算逼近多项式与目标函数之间的残差
residuals = exp(0:0.01:1) - arrayfun(@(x) sum(a.*chebyshev(x, 0:n-1)'), 0:0.01:1);
%% 计算平均绝对误差（MAE）
mae = sum(abs(residuals))/length(residuals);
%% 计算均方根误差（RMSE）
rmse = sqrt(sum(residuals.^2)/length(residuals));
fprintf('平均绝对误差（MAE）：%f\n', mae);
fprintf('均方根误差（RMSE）：%f\n', rmse);
\end{python}
\textbf{附件 G : 三种复化公式的实现评判代码 }
\begin{python}
clear; clc;
format long
% 定义函数
fun = @(x) exp(x);
% 精确值计算
N = 100000000; % 积分划分的子区间数目
h = 1/N; % 子区间宽度
x = 0:h:1; % 子区间节点
y = fun(x); % 子区间对应的函数值
jingque = sum(y(1:end-1)) * h;
fprintf("n 复化梯形公式 复化Simpson公式 复化Cotes公式");
fprintf("\n")
for i = 2:8
n = 2^i;
h = 1/n;
t = 0;
s = 0;
c = 0;
for k = 0:n-1
t = t + fun(k*h) + fun((k+1)*h);
end
t = h * t / 2;
for k = 0:2:n-2
s = s + fun(k*h) + 4 * fun((k+1)*h) + fun((k+2)*h);
end
s = h * s / 3;
for k = 0:4:n-4
c = c + 7 * fun(k*h) + 32 * fun((k+1)*h) + 12 * fun((k+2)*h) + 32 * fun((k+3)*h) + 7 * fun((k+4)*h);
end
c = 2 * h * c / 45;
if n < 10
fprintf("%d %5.15f %5.15f %5.15f\n", n, t, s, c)
elseif n < 100
fprintf("%d %5.15f %5.15f %5.15f\n", n, t, s, c)
else
fprintf("%d %5.15f %5.15f %5.15f\n", n, t, s, c)
end
end
fprintf("精确值：%3.15f\n", jingque);
biao_value = 1.71828182845904509;
% 计算绝对误差
abs_error = abs(jingque - biao_value);
% 计算相对误差
rel_error = abs_error / biao_value;
% 输出结果
disp(['绝对误差：', num2str(abs_error)]);
disp(['相对误差：', num2str(rel_error)]);
\end{python}
\textbf{附件 H : $Romberg$算法前两行公式的代码实现作业要求 }
\begin{python}
% 清除命令窗口和工作空间的内容
clear, clc;
% 定义被积函数
f = @(x) exp(x);
% 初始化参数
a = 0; % 积分下限
b = 1; % 积分上限
Tprev = (f(a) + f(b)) / 2; % 上一层级的近似值
j = 1; % 迭代次数
% 计算公式近似值和误差
exact = 1.718281828459;
error = abs(exact - Tprev);
fprintf('第1次迭代的近似值: %.12f\n', Tprev);
errors = []; % 存储每次迭代的误差值
while error >= 10^(-10)
h = (b - a) / (2^j);
sum_f = 0;
for i = 1:2^(j-1)
sum_f = sum_f + f(a + (2*i - 1)*h);
end
Tcurr = Tprev/2 + h*sum_f;
% 输出每次迭代的结果
fprintf('第%d次迭代的近似值：%.12f\n', j+1, Tcurr);
error = abs(exact - Tcurr);
Tprev = Tcurr;
errors = [errors error];
j = j + 1;
end
% 绘制误差值的曲线图
plot(errors);
title('每次迭代的误差值');
xlabel('迭代次数');
ylabel('误差');
\end{python}
\textbf{附件 I : $Romberg$算法前两行公式的代码实现作业要求Python实现 }
\begin{python}
import numpy as np
import matplotlib.pyplot as plt
# 定义被积函数
def f(x):
return np.exp(x)
# 初始化参数
a = 0 # 积分下限
b = 1 # 积分上限
Tprev = (f(a) + f(b)) / 2 # 上一层级的近似值
j = 1 # 迭代次数
# 计算公式近似值和误差
exact = 1.718281828459
error = abs(exact - Tprev)
print('第1次迭代的近似值:', Tprev)
errors = [] # 存储每次迭代的误差值
values = [] # 存储每次迭代的计算值
while error >= 10**(-10):
h = (b - a) / (2**j)
sum_f = 0
for i in range(1, 2**(j-1)+1):
sum_f += f(a + (2*i - 1)*h)
Tcurr = Tprev/2 + h*sum_f
# 输出每次迭代的结果
print('第{}次迭代的近似值：{}'.format(j+1, Tcurr))
error = abs(exact - Tcurr)
Tprev = Tcurr
errors.append(error)
values.append(Tcurr)
j += 1
# 绘制误差值和每次迭代的统计图
plt.figure(figsize=(12, 6))
# 误差值的曲线图
plt.subplot(1, 2, 1)
plt.plot(errors, color='blue', linewidth=2, label='Error')
plt.title('Error at Each Iteration')
plt.xlabel('Number of Iterations')
plt.ylabel('Error')
plt.grid(True) # 添加网格线
plt.legend() # 添加图例
# 每次迭代的统计图
plt.subplot(1, 2, 2)
plt.plot(values, color='red', marker='o', markersize=5, label='Values')
plt.title('Computed Values at Each Iteration')
plt.xlabel('Number of Iterations')
plt.ylabel('Computed Value')
plt.grid(True)
plt.legend()
plt.tight_layout() # 调整子图布局
plt.show()
\end{python}
\textbf{附件 J : $Romberg$算法的实现评判代码 }
\begin{python}
clear, clc;
% 设置参数
a = 0; % 积分下限
b = 1; % 积分上限
n = 5; % 迭代次数
% 初始化 Romberg 矩阵
R = zeros(n);
tic;
% 计算 h 值
h = (b-a) ./ (2.^(0:n-1));
% 计算 T(1,1)
R(1,1) = (h(1)/2) * (exp(a) + exp(b));
% 计算 T(i,1)
for i = 2:n
sum = 0;
for j = 1:2^(i-2)
sum = sum + exp(a + (2*j-1)*h(i));
end
R(i,1) = 0.5 * R(i-1,1) + h(i) * sum;
end
% 计算 T(i,j)
for j = 2:n
for i = j:n
R(i,j) = R(i,j-1) + (R(i,j-1) - R(i-1,j-1)) / ((4^j)-1);
end
end
toc;%迭代运行时间
% 表5.5给定数值得到T
T = [1.8591409142 0 0 0 0;
1.7539310925 1.7188611519 0 0 0;
1.7272219046 1.7183188419 1.7182826879 0 0;
1.7205185922 1.7182841547 1.7182818422 1.7182818288 0;
1.7188411286 1.7182819741 1.7182818287 1.7182818285 1.7182818285];
% 输出每个点的误差
for i = 1:n
for j = 1:i
error = abs(R(i,j) - T(i,j));
fprintf('第 %d 行，第 %d 列： Romberg 积分结果 %.10f，参考答案 %.10f，误差 %.10f\n', i, j, R(i,j), T(i,j), error);
end
end
fprintf('\n精确值为： ')
% 输出结果
result = R(n,n);
disp(result);
% 计算均方误差
mse = mean((R(:) - T(:)).^2);
fprintf('均方误差：%.10f\n', mse);
% 计算绝对误差
abs_error = abs(1.718701740929343 - 1.718281828459);
% 计算相对误差
rel_error = abs_error / 1.718281828459;
% 输出结果
disp(['绝对误差：', num2str(abs_error)]);
disp(['相对误差：', num2str(rel_error)]);
\end{python}
\end{appendices}
\end{document}需要注意，要把如上代码直接在overleaf跑动会报错，因为你没有导入源图片文件，这个我不提供，不过大家可以参考latex附录里的程序来复现类似的图。
具体的PDF资源还在上传中，上传好了关注我下载资源即可ww.
扫一扫
942
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
