注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

太阳神上的博客

青青子衿,悠悠我心,但为君故,沉吟至今。

 
 
 

日志

 
 
 
 

数值算法设计与实现笔记(简体中文版)  

2007-06-21 17:02:45|  分类: 学习 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

\documentclass{ctexart}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{times,amsfonts,algorithm,algorithmic}

\usepackage[bookmarks=true]{hyperref}
% setup hyerref for PDF navigation
\hypersetup{colorlinks, linkcolor=blue, citecolor=blue, urlcolor=black,%
    plainpages=true, bookmarksopen=false,%
    pdfhighlight=/P, %/I(inverse) /N(no effect) /O(outline) /P(inset)
    pdfauthor={Moligaloo <moligaloo@163.com>},%
    pdftitle={Notes on Numerical Algorithm and Implementation},%
    pdfsubject={NAI},%
    pdfkeywords={NAI},%
    pdfstartview=FitH, %FitBH, FitB
    pdfview=FitH,
    pdfpagemode=UseOutlines, %None, FullScreen, UseThumbs
}

\title{献给天使上--数值算法设计与实现笔记\\(简体中文版)}
\author{太阳神上}

\begin{document}
\maketitle
\tableofcontents
\newpage

\floatname{algorithm}{算法 }

\section{算法}

\subsection{算法的定义}
什么是算法?在算法的经典著作《算法导论》(由麻省理工学院出版社出版)里,有如下的定义:

\begin{quotation}
简单地讲,\textbf{算法(Algorithm)}就是以一些值(或值的集合)为\textbf{输入(input)},经过一系列定义好了的计算,
产生一些值作为\textbf{输出(output)}的计算过程,它包含一系列计算步骤来将输入转化为输出。
\end{quotation}

在算法大师高德纳(Donald Erwin Knuth)著作里,算法有着更为严格的定义,它是一个诸如$(\Omega,I,O,f)$的四元组的集合。
不过由于此定义太复杂,以我们大家的水平难以理解,所以我们还是从算法的外延而不是内涵来理解吧。

\subsection{算法的性质(外延)}
\begin{enumerate}
  \item 有穷性,算法一定要能结束
  \item 0个或以上输入
  \item 1个或以上输出
  \item 确定性,即每个步骤都是无二义的
  \item 可行性(对于计算机而言)
\end{enumerate}

\subsection{算法的描述}
\begin{description}
  \item[流程图] 不过现在已经过时了!\footnote{不要相信CCS!}
  \item[INFL语言] (INFL是INFormaL的简写) 就是曹老师上课用的那种。
  \item[高级语言] C/C++,Java,Pascal,C\#\ldots (命令式语言),或者Lisp,Scheme,ML和 Haskell(函数式语言)。
\end{description}

\subsection{算法的效率}
\textbf{数值}算法的效率我们不是用其增长函数如$O(n^2)$来描述的,而是用如下的形式:
$$ a\mathbf{A}+b\mathbf{M}+c\mathbf{D} $$
$\mathbf{A}$,$\mathbf{M}$,和 $\mathbf{D}$ 分别代表\textbf{加法(Addition)},\textbf{乘法 (Multiplication)}
和\textbf{除法(Division)},而 $a$,$b$ 和$c$则是其相应的次数,因为在计算机里,计算加法和减法是差不多的,
所以减法可以算入加法里。

\subsection{数值算法的思想}
\begin{enumerate}
  \item \textbf{离散化方法} 即将连续量转化为离散量,因为计算机只能处理离散量。典型的例子就是用梯形来计算积分。
  \item \textbf{迭代方法} 因为计算机的计算速度很快,因为重复简单计算以求精是很适用于计算机的。
\end{enumerate}

下面是迭代方法的基本步骤:
先构造出一个方程:
$$ x_{max}=\Phi(x_n) $$
使得\[ \lim_{n\rightarrow \infty} x_n  \]会收敛到问题的解。

用迭代方法来计算$\sqrt{2}$:

因为$\sqrt{2}$是$f(x)=x^2-2=0$的正根
$$
\begin{array}{lll}
x_{n+1} & = & x_n-\displaystyle{\frac{f(x_n)}{f'(x_n)}} \\
& = & x_n-\displaystyle{\frac{x_n^2-2}{2x_n}} \\
& = & \displaystyle{\frac{1}{2}x_n+\frac{1}{x_n}} \\
\end{array}
$$

于是我们得到以下公式:
$$
\left\{
\begin{array}{lll}
x_0 & = & 2 \\
x_{n+1} & = & \displaystyle{\frac{x_n}{2}+\frac{1}{x_n}}
\end{array}
\right.
$$
根据上面的递推公式,我们可以算得:
$$
\begin{array}{lll}
x_0 & = & 2 \\
x_1 & = & \displaystyle{\frac{x_0}{2}+\frac{1}{x_0} = \frac{2}{2}+\frac{1}{2}}=1.5 \\
x_2 & = & \displaystyle{\frac{x_1}{2}+\frac{1}{x_1} = \frac{1.5}{2}+\frac{1}{1.5}}=0.75+0.666=1.416 \\
x_3 & = & \displaystyle{\frac{x_2}{2}+\frac{1}{x_2} =
\frac{1.416}{2}+\frac{1}{1.416}}=0.708+0.706=1.414
\end{array}
$$
计算到$x_3$时,我们得到了$\sqrt{2}$的一个近似值。

\subsection{误差}
\subsubsection*{绝对误差}
$$ e=|x-x^*| $$
$x$是测量值, $x^*$是准确值。而在实际中,绝对误差是无意义的。

\subsubsection*{相对误差}
$$ e=\frac{|x-x^*|}{|x^*|}  $$
因为$x$和$x^*$相差不大,所以下式和上式近似
$$ e=\frac{|x-x^*|}{|x|} $$
相对误差也就是\textbf{误差限}。

\subsubsection*{切断误差}
这是由于计算机表示数据的有限性引起的,因为计算机只能表示有理数。

\subsubsection*{舍入误差}
这是由于四舍五入引起的误差。

\subsubsection*{精度(误差范围)}
“事实上,我们平时生活中所提到的误差就是精度”(曹老师语)。

\subsection{有效数字}
由四舍五入产生的数字属于有效数字,如:
$$ \pi=3.1415926535 8979323846 2643383279 5028841971 6939937510\ldots $$
\begin{center}
\begin{tabular}{|l|c|}
  \hline
  近似值 & 有效数字位数 \\
  \hline
  3.14 & 三位 \\
  3.141 & 三位 \\
  3.142 & 四位 \\
  3.1415 & 四位 \\
  3.14159 & 六位 \\
  \hline
\end{tabular}
\end{center}

\subsection{适定性}
如果微小的输入改变会引起结果的巨大变化(蝴蝶效应),那么它就是不适定的。

适定性有两个方面,问题本身的适定性与算法的适定性。

\textbf{混沌问题(Chaos Problem)}就是一个典型的非适定性问题。

下面是一个非适定性的例子:
\begin{equation*}
\left\{
 \begin{aligned}
x+2y &= 3 \\
0.499x+y&=1.5
\end{aligned}
\right.
\end{equation*}

下面也是:
$$ f(x)=\frac{1}{x-\sqrt{x^2-1}} $$
当$x$很大时,$x$和$\sqrt{x^2-1}$的差别就很小了,结果一除,问题就来了。
解决方法是将分母有理化如下:
$$ f(x)=\frac{1(x+\sqrt{x^2-1})}{(x-\sqrt{x^2-1})(x+\sqrt{x^2-1})}=\frac{x+\sqrt{x^2-1}}{x^2-(\sqrt{x^2-1})^2}= x+\sqrt{x^2-1} $$
这样就没问题了。

下面是一个更复杂的计算积分的方法:
\newcommand{\ }{\textrm{d}}
$$ I=\int_0^1e^{x-1}x^{20}\ x $$

首先我们先弄出一个数学模型如下:
\begin{eqnarray*}
\textrm{设} I_n & = & \int_0^1 x^n e^{x-1} \ x \\
& = & \int_0^1 x^n (e^{x-1})^1 \ x \\
& = & x^n e^{x-1} \Big|_0^1 - \int_0^1 nx^{n-1}e^{x-1} \ x \\
& = & 1-nI_{n-1}
\end{eqnarray*}
得到以下递推公式:
$$
\left\{
\begin{array}{lll}
I_n & = & 1-nI_{n-1} \\
I_0 & = & \displaystyle{1-\frac{1}{e}}
\end{array}
\right.
$$
尽管上面的公式是一个完美的数学模型,但是,如果我们按照上面的方法来计算。
误差就会严重扩大,因为每一步都会将上一次的误差扩大,
结果算到$n=20$时,它就被放大了$20!$次,这可是一个巨大的天文数字。

因此,上面是一个非适定的算法,适定的算法如下
$$
\left\{
\begin{array}{lll}
I_{40} & = & 1 \\
I_{n-1} & = & \displaystyle{\frac{1-I_n}{n},n=40,39,\ldots 21}
\end{array}
\right.
$$
这在数学上一个错误的模型,可是在本题里却能获得满足要求的结果,可是数学理论和实践是有很大区别的。
下面是一个总结:
\begin{quotation}
数值算法的不稳定性的来源是\textbf{误差的扩大}。
\end{quotation}


\section{搜索算法}
\subsection{黄金分割搜索算法}
\textbf{黄金分割搜索算法}就是迭代方法的一种具体体现。黄金分割率是由下面的方程解出来的:
$$ x^2+x-1=0 $$
于是得到:
$$ \varphi=\frac{\sqrt{5}-1}{2}\approx 0.618 $$
黄金分割搜索算法用于求最大值。

下面是相关的算法伪代码:

\begin{algorithm}
\caption{黄金分割搜索算法}
\begin{algorithmic}[1]
\STATE Input$\{a,b,e\}$
\STATE $x_1 \leftarrow a+(b-a)*0.382$
\STATE $x_2 \leftarrow a+(b-a)*0.618$
\STATE $f_1 \leftarrow \cos (x_1)$
\STATE $f_2 \leftarrow \cos (x_2)$
\WHILE{$b-a>e$}
\IF{$f_2>f_1$}
\STATE $a \leftarrow x_1$
\STATE $x_1 \leftarrow x_2 $
\STATE $f_1 \leftarrow f_2 $
\STATE $x_2 \leftarrow a+0.618*(b-a)$
\STATE $f_2 \leftarrow \cos(x_2)$
\ELSE
\STATE $b \leftarrow x_2$
\STATE $x_2 \leftarrow x_1 $
\STATE $f_2 \leftarrow f_1 $
\STATE $x_1 \leftarrow a+0.382*(b-a)$
\STATE $f_1 \leftarrow \cos(x_1) $
\ENDIF
\ENDWHILE
\STATE Output$\{\displaystyle{\frac{a+b}{2}}\}$
\end{algorithmic}
\end{algorithm}

\subsection{二分搜索算法}
二分搜索用于求解在一段区间上\textbf{单调}且\textbf{连续}的函数的解,
在上学期的数据结构课里,我们也可以用这种方法来在一个已排好了序的数组里找元素。

下面是该算法的伪代码,设$f(x)$在区间$[a,b]$单调递增,输出值为当$f(x)=0$时的解。
\begin{algorithm}
\caption{二分搜索算法}
\begin{algorithmic}[1]
\STATE Input$\{a,b,e\}$
\WHILE{$b-a>e$}
\STATE $mid \leftarrow \displaystyle{\frac{a+b}{2}} $
\IF{$f(mid)>0$}
\STATE $b \leftarrow mid$
\ELSE
\STATE $a \leftarrow mid$
\ENDIF
\ENDWHILE
\STATE Output $\{mid\}$

\end{algorithmic}
\end{algorithm}


\subsection{爬山算法}

\subsubsection{问题}
设$f(x) \in \Omega,x \in \mathbb{R}^n$,($f(x)$只有单峰)求:
$$ x_{max}=\{x|\forall y \in \mathbb{R}^n,f(y)\leq f(x) \} $$
也就是多元函数求最大值。
$x_{k+1}=x_k+\lambda \nabla f(x_k) $.
$\lambda$为步长,最大(上山,+),最小(下山,-)。
它的基本思想是着梯度方向爬,此算法是迭代思想的一种体现。

\subsubsection{方法}
\begin{enumerate}
  \item 给定初值$x^{(0)}$及精度$\epsilon$,若$||\nabla f(x^{(0)})||^2 \leq \epsilon||$,则$x^{(0)}$即为近似极小值。
  $\nabla f(x)=\displaystyle{(\frac{\partial f}{\partial x_1},\cdots,\frac{\partial f}{\partial x_n})^T},||x||^2=x_1^2+x_2^2+\cdots+x_n^2$
  \item 若$||\nabla f(x^{(0)})||^2>\epsilon$,用适当步长$\lambda_0$代入下式计算:
  $$
  x^{(1)}=x^{(0)}-\lambda_k \nabla f(x^{(0)})
  $$
  \item 一般,若$|| \nabla x^{(k)} || ^2 \leq \epsilon$,则$x^{(k)}$为近似极小值,否则用适当步长$\lambda_k$确定下一个近似值,
  直到满足精度为止:
  $$
  x^{(k+1)}=x^{(k)}-\lambda_k \nabla f (x^{(k)})
  $$
\end{enumerate}

\section{插值逼近算法}

\subsection{问题描述}
\textbf{插值逼近}就是将\textbf{离散}的函数转化为\textbf{连续}的函数,这种做法广泛运用于服装设计和造船工业上。

顺便说一句,曹老师还补充了数学家感兴趣的三个方面:
\begin{enumerate}
\item 问题是否有解?
\item 问题有一个还是多个解?
\item 算法是不是适定的?
\end{enumerate}

\subsection{拉格朗日插值算法}
有$n+1$个离散点,它们是
$$ (x_0,y_0),(x_1,y_1)\ldots (x_n,y_n)$$
并且$x_0<x_1<\ldots<x_n$,怎样获得一个$n$次多项式$P_n(x)$使得:
$$
y_k=P_n(x_k),k=0,1,\ldots,n
$$

方法是:
$$
f(x)=\sum_{k=0}^n \Big(y_k \prod_{i=0,i\neq k}^n \frac{x-x_i}{x_k-x_i}\Big)
$$
上面的括号是我自己加的。

该问题存在解且唯一,范德蒙德行列式$\pi(x_i-x_j) \Rightarrow$
由于$x_0,\ldots,x_n$互异,因此行列式不为0$\Rightarrow$解存在且唯一。

\subsubsection{实现}
\begin{algorithm}
 \caption{拉格朗日插值算法的实现}
\begin{algorithmic}[1]
 \STATE Input$\{N,x_0,y_0,\ldots,x_n,y_n\}$
 \REPEAT
 \STATE Input $\{x\}$
 \STATE $S \leftarrow 0$
 \FOR{$k=0,1,2,\ldots,N$}
 \STATE $p=y_k$
 \FOR{$i=0,1,2,\ldots,N$}
 \IF{$i\neq k$}
 \STATE $p\leftarrow \displaystyle{\frac{p(x-x_i)}{x_k-x_i}}$
 \ENDIF
 \ENDFOR
 \STATE $S\leftarrow S+p$
 \ENDFOR
 \STATE Output $\{y\}$
 \UNTIL{No more $x$}
\end{algorithmic}
\end{algorithm}

\subsubsection{效率}
$$
(2n+1)(n+1)\mathbf{A}+n(n+1)\mathbf{M}+n(n+1)\mathbf{D}
$$

\subsubsection{提高算法效率}
一种提高拉格朗日插值算法的效率的方法是使用秦九韶算法(中国人管它叫秦九韶算法,西方人叫\emph{Hornor's Rule})。


\begin{equation*}
\begin{aligned}
P_n(x) &= P_{n-1}(x)+Q_n(x) \\
Q_n(x) &= P_n(x)-P_{n-1}(x) \\
Q_n(x_k) &= 0,k=0,1,2,\ldots,n \\
\end{aligned}
\end{equation*}
$\Rightarrow Q_n(x)=a_n(x-x0)\cdots(x-x_{n-1}) a_n$是$P_n(x)$的系数
$$
\therefore a_n=\sum_{k=0}^{n}\big(y_k\prod_{i=0,i\neq k}^{n}\frac{1}{x_k-x_i}\big)
$$

\subsection{牛顿插值算法}
牛顿插值算法的公式是:
$$
P_n(x)=a_0+a_1(x-x_0)+\cdots+a_n(x-x_0)(x-x_1)\cdots(x-x_{n-1})
$$
\begin{algorithm}
 \caption{牛顿插值算法}
\begin{algorithmic}[1]
\STATE Input$\{x_0,x_1,\ldots,x_n,a_1,a_2,\ldots,a_n \}$
\REPEAT
\STATE Input $\{x\}$
\STATE $p \leftarrow a_n$
\FOR{$i=n-1,n-2,\ldots,a_n$}
\STATE $p\leftarrow p(x-x_i)+a_i$
\ENDFOR
\STATE Output$\{p\}$
\UNTIL{No more $x$}
\end{algorithmic}
\end{algorithm}

下面是得到系数$a_n$的算法。
\begin{algorithm}
 \caption{得到系数$a_n$的算法}
\begin{algorithmic}[1]
 \STATE Input $\{x_0,y_0,\ldots,x_n,y_n\}$
 \FOR{$j=0$ to $n$}
 \STATE $a_j \leftarrow y_j$
 \ENDFOR
 \FOR{$k=1$ to $n$}
 \FOR{$j=n$ to $k$ \textbf{step} $-1$}
 \STATE $a_j \leftarrow \displaystyle{\frac{a_j-a_{j-1}}{x_j-x_{j-k}}}$
 \ENDFOR
 \ENDFOR
 \STATE Output$\{a_0,a_1,\ldots,a_n\}$
\end{algorithmic}
\end{algorithm}
算法的效率:
$$
2n\mathbf{A}+n\mathbf{M}
$$
系数计算次数
$$
\sum_{k=1}^{n}2(n-k+1)\mathbf{A}+\sum_{k=1}^{n}(n-k+1)\mathbf{D}
$$

牛顿插值算法和拉格朗日插值算法所画出来的曲线是一样的。

用拉格朗日插值算法和牛顿插值算法都有\textbf{龙格现象(Runger Phenomenon)},
多项式中的非适定性,因此要让阶数尽量低。从这个思想出发,于是有了我们的三次自然样条插值算法。

\subsection{三次自然样条插值算法}
\subsubsection{定义}
前面两个算法只具有教学价值,而工业上主要使用的是三次自然样条插值算法。
下面是这个算法的四个要点:
\begin{enumerate}
 \item 通过已知的离散点(\textbf{B样条插值算法}会更平滑,但是不一定过已知离散点)。
 \item 具有二阶连续函数(弯距,无集中扭距,不用应力)。
 \item 分段三次多项式。
 \item 具有最小弹性势能。
\end{enumerate}

给定值表函数
$$
{x_0,x_1,\cdots,x_n} \choose {y_0,y_1,\cdots,y_n}
$$
其中
$$
a=x_0 < x_1 < x_2 < \cdots < x_n =b
$$
则满足以下四个条件的函数$S(x)$被称为三次自然样条:
\begin{enumerate}
  \item $S(x_k)=y_k,k=0,1,\ldots,n$
  \item $S(x),S'(x),S''(x) \in \mathbf{C}[a,b]$
  \item 在每个子区间$[x_k,x_{k+1}]$上$S(x)$是一个三次多项式,$k=0,1,\ldots,n$
  \item 对于任意满足1,2的函数$g(x)$有
  $$
  \int_a^b[S''(x)]^2 \ x \neq \int_a^b[g''(x)]^2 \ x \Rightarrow S''(a)=S''(b)=0
  $$
\end{enumerate}

\subsubsection{三次自然样条插值算法的导出}
(略)

\begin{algorithm}
 \caption{三次自然样条插值算法}
\begin{algorithmic}[1]
 \STATE Input$\{ x_0,y_0,\ldots,x_n,y_n\}$
 \FOR{$k=0,1,2,\ldots,n-1$}
 \STATE $h_k \leftarrow x_{k+1}-x_k$
 \ENDFOR
 \STATE $a_1 \leftarrow 2*(h_0+h_1)$
 \FOR{$k=2,3,\ldots,n-1$}
 \STATE $a_k\leftarrow 2*(h_{k-1}+h_k)-h_{k-1}^2/a_{k-1}$
 \ENDFOR
 \FOR{$k=1,2,\ldots,n$}
 \STATE $c_k\leftarrow(y_k-y_{k-1})/h_{k-1}$
 \ENDFOR
 \FOR{$k=1,2,\ldots,n-1$}
 \STATE $d_k=6*(c_{k+1}-c_k)$
 \ENDFOR
 \STATE $b_1\leftarrow d_1$
 \FOR{$k=2,3,\ldots,n-1$}
 \STATE $b_k=d_k-b_{k-1}*h_{k-1}/a_k$
 \ENDFOR
 \STATE $S''_{n-1}\leftarrow b_{n-1}/a_{n-1}$
 \FOR{$k=n-2,n-3,\ldots,2,1$}
 \STATE $S''_k=(b_k-h_kS''_{k+1})/a_k$
 \ENDFOR
 \STATE $S''_0\leftarrow 0$
 \STATE $S''_n\leftarrow 0$
 \FOR{$k=0,1,2,\ldots,n-1$}
 \FOR{$x=x_k$ to $x_k+1$}
 \STATE $S'\leftarrow c_{k+1}-S''_{k+1}h_k/6-S''_k h_k/3$
 \STATE $S\leftarrow y_k+S'(x-x_k)+S''_k(x-x_k)^2/2+(S''_{k+1}-S''_k)*(x-x_k)^3/(6h_k)$
 \STATE Output $\{ S \}$
 \ENDFOR
 \ENDFOR
\end{algorithmic}
\end{algorithm}

\section{人工神经网络}
\subsection{人工神经元的基本思想}
单神经元
$$
\sigma = \sum w_k x_k
$$
$w$为权值,当$\sigma$大于一个值时,神经元兴奋,反之则抑制。

$$
y=
\left\{
\begin{array}{lll}
1  , \sigma > \theta \\
0  , \sigma \leq \theta
\end{array}
\right.
$$

其中
$$
\sigma(x)=\frac{1}{1+e^{-x}}
$$
$\sigma(x)$学名叫sigmoid函数。

\subsection{BP神经网络模拟黑箱子方法}
典型的BP网络是三层网络,包括输入层,隐含层和输出层。
$$
\sigma(\sum_k\bar{w}_{ki}\sigma(\sum w_{ki}-\theta_k)-\bar{\theta_i})=y_j
$$

\subsection{如何学习权值}
$$
x_1^{(k)}\cdots x_m^{(k)} \rightarrow y_1^{(k)} \cdots y_m^{(k)},k=1\cdots P
$$
误差$\epsilon=\sum_{q=1}^{m} (\bar{y}_q^{(k)}-y_q)^2$

$\sigma(x)$的导数:
$$
\sigma'(x)=\frac{e^{-x}+1-1}{(1+e^x)^2}=\frac{1}{1+e^{-x}}-\frac{1}{(1+e^{-x})}=\sigma(x)(1+\sigma(x))
$$


\section{随机数产生算法}
随机数在科学实验中运用得很广,一般的科学实验都得先在计算机模拟,没问题后才能真正到实验室去上做。

\subsection{线性同余法}

\subsubsection{数学模型}
线性同余法是一种很简单地产生均匀分布随机数的方法。

给定三个质数$a,b,M$,其中$M$很大。和一个种子$n_0$,可以根据下式得到一个随机数列。
$$
n_{k+1}=(a*n_k+b) \bmod M
$$

\subsubsection{C语言的实现}
在C语言中,使用\textbf{srand(unsigned int)}来设定随机种子,用\textbf{rand()} 来返回一个从0到0x7FFF(32767)之间的随机数。
一般使用当前的时间来作为随机种子,在C语言中,time(NULL)返回从1970年1月1日到程序运行到这个函数时所经过的秒数。

\subsubsection{0-1均匀分布随机数的数值特征}
数学期望:
$$
E(\xi)=\int_{-\infty}^{\infty}xf(x)\ x=\int_0^1 x\ x=\frac{1}{2}x^2|_0^1=\frac{1}{2}
$$
方差:
$$
\begin{array}{lll}
D(\xi) & = & \displaystyle{ \int_{-\infty}^{\infty}(x-E(\xi))^2f(x)\ x }\\
& = & \displaystyle{ \int_0^1(x-\frac{1}{2})^2 \ x }\\
& = & \displaystyle{ \frac{1}{3}(x-\frac{1}{2})^3|_0^1 }\\
& = & \displaystyle{ \frac{1}{3}((\frac{1}{2})^3-(-\frac{1}{2})^3) }\\
& = & \displaystyle{ \frac{1}{12} }
\end{array}
$$
平均误差:
$$
\sqrt{\frac{1}{12}}
$$

\subsection{正态分布随机数}

\subsubsection{背景}
一般的测量误差都服从正态分布。

\subsubsection{计算模型}
\textbf{中心极限定理}:\\
设 $\{\xi_k\}$ 是0-1均匀分布随机数列,
$$
\eta=\lim_{n\rightarrow \infty}\frac{\xi_1+\xi_2+\cdots+\xi_n}{n}
$$
则$\eta$服从正态分布。

定义rnd()为rand()/32767.0,即rnd()返回0-1之间均匀分布的随机数。
$$
E(\xi)=\frac{1}{n}\sum_{k=1}^n E(\textrm{rnd}())=\frac{1}{2}
$$

$$
D(\xi)=(\frac{1}{n})^2 \sum_{k=1}^n D(\textrm{rnd}())=\frac{1}{12n}
$$

\section{积分算法}
\subsection{数值积分算法的计算模型}
\subsubsection{梯形模型}
计算
$$
I=\int_a^b f(x)\ x,\textrm{设}
h=\frac{b-a}{n}
$$
$a$和 $b$分别是积分下限和积分上限,$n$是分的段数, $h$是步长。
通过$I$下面的公式可以算得其积分。
$$
I \approx \sum_{k=1}^{h} \frac{h[f(a+kh)+f(a+(k-1)h)]}{2}
$$
我们并不能通过增加$n$的值来使得结果更精确,因为当$n$足够大,分成的梯形越小时,
计算同舍入误差也就越大。

\subsubsection{复化递推法}

$$
I_n=\sum_{k=1}^{n} \frac{h[f(a+(k-1)h)+f(a+kh)]}{2}
$$
于是
$$
I_{2n}=\frac{1}{2}I_n+\frac{n}{2}\sum_{k=1}^{n}f(x_{k-\frac{1}{2}})
$$
$$
I_{\frac{h}{2}}=\frac{1}{2}I_h+\frac{h}{2}\sum_{k=1}^{n}f(x_{k-\frac{1}{2}})
$$
引理(lemma):
$$
|I_{2n}-I_n|<\epsilon \Rightarrow |I-I_{2n}|<\frac{1}{3}\epsilon
$$

\subsection{变步长数值积分算法}
\begin{algorithm}
 \caption{变步长数值积分算法}
\begin{algorithmic}[1]
 \STATE Input$\{a,b,e\}$
 \STATE $h\leftarrow b-a$
 \STATE $T_0 \leftarrow \displaystyle{\frac{h(f(a)+f(b))}{2}}$
 \STATE $T_1 \leftarrow \displaystyle{\frac{T_0}{2}+h*f(a+\frac{h}{2})}$
 \STATE $n \leftarrow 1$
 \WHILE{$|T_1-T_0| \geq e$}
   \STATE $h \leftarrow h/2$
   \STATE $n \leftarrow n*2$
   \STATE $T_0 \leftarrow T_1$
   \STATE $T_1 \leftarrow 0 $
   \FOR{$k=1$ to $n$}
      \STATE $T_1 \leftarrow T_1+f(a+(k-\frac{1}{2})*h)$
   \ENDFOR
 \ENDWHILE
 \STATE Output$\{T_1\}$
\end{algorithmic}
\end{algorithm}

\subsection{自适应数值积分算法}
这个算法使用了堆栈的思想。
$$
T_1=\frac{(f(a)+f(b))h}{2}
$$
$$
T_2=\frac{T_1}{2}+\frac{f(\frac{a+b}{2})h}{2}
$$

下面是其伪代码
\begin{algorithm}
\caption{自适应数值积分算法}
\begin{algorithmic}[1]
\STATE Input$\{a,b,e\}$
\STATE $E \leftarrow \displaystyle{\frac{e}{b-a}}$
\STATE $I \leftarrow 0$
\STATE $S \leftarrow 0$
\STATE $U \leftarrow a$
\STATE $V \leftarrow b$
\STATE \textbf{start:}
\STATE $h \leftarrow V-U$
\STATE $T_1 \leftarrow \displaystyle{\frac{(f(U)+f(V))*h}{2}}$
\STATE $W \leftarrow \displaystyle{\frac{U+V}{2}}$
\STATE $T_2 \leftarrow \displaystyle{\frac{T_1}{2}}+\displaystyle{\frac{h*f(W)}{2}}$;
\IF{$|T_2-T_1|>=E*h$}
\STATE $I \leftarrow I+1$
\STATE $ST[I] \leftarrow V$
\STATE $V \leftarrow W$
\STATE \textbf{goto start}
\ENDIF
\STATE $S \leftarrow S+T_2$
\IF{$I \neq 0$}
\STATE $U \leftarrow V$
\STATE $V \leftarrow ST[I]$
\STATE $I \leftarrow I-1$
\STATE \textbf{goto start}
\ENDIF
\STATE Output$\{S\}$
\end{algorithmic}
\end{algorithm}
\end{document}

  评论这张
 
阅读(827)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017