julia 数值分析

发表于 2023-09-12 更新于 2023-11-05 Disqus：

用julia实现一些数值方法不过这些方法julia肯定有现成的库和函数解决了我只是练练手瞎写一些东西

浮点数和多项式

多项式求值

function hornor(x::Float64,n::Int64,c::Vector{Float64}=[1. for  i in 1:n+1],b::Vector{Float64}=[0. for  i in 1:n+1])
	y = c[end]
	for i in n:-1:1
		y = y*( x- b[i])+c[i]
	end

	return y

end

实现浮点数表示

\[\pm 1 . b b b \cdots b \times 2^p\]

存储时需要尾数 bbbb 和指数p

IEEE舍入最近法则
对于浮点精度，如果在二进制数右边的第 53 位是0，则向下舍去(在第 52 位后面截断)。如果第53位是1，则向上进位(在第 52位上加 1)，除了在1右边的所有已知位都是0，在这种情况下当且仅当第52位是1时在第52位上加1。

一个有趣的快速求平方根倒数算法视频链接

一个浮点数y在二进制里存贮形式Y 是 S E M（S为0和1表示符号 0为正部分 M为0和1 表示尾数这里省略了小数点前的1. E为为指数也有0和1表示不过它有正负使用实际值是补码)

$y=\left(1+\frac{M}{2^{23}}\right) \cdot 2^{E-127}$ M是尾数部分的值 E是指数部分的值

举个例子 y=1.5 Y 是 0(S表示正数) 01111111(E表示0 ) 10000000000000000000000(M本来应该是1.10000000000000000000000)

\[ log_2(y)=&log_2\left(1+\frac{M}{2^{23}}\right) + E-127\\ log_2(y) \approx &\left(\frac{M}{2^{23}}\right) + E-127\\ log_2(y) \approx &2^{-23}\left(M+2^{23}*E\right) -127\\ log_2(y) \approx &2^{-23}Y -127 \]
这里考虑到$log_2\left(1+x\right)$曲线在0和1之间接近x曲线用x 代替了$log_2\left(1+x\right)$ 当然这种近似仅在0和1处最好所以加一个修正可以是整体误差最小后面会在代码里看到修正
\[ log_2(a)=&-\frac{1}{2} log_2(y)\\ 2^{-23}A -127 \approx &-\frac{1}{2}\left(2^{-23}Y -127\right)\\ A \approx &-\frac{1}{2}Y +2^{22}\left( 381\right) \]
这里a是y的平方根倒数A是a的二进制表示

对于$\frac{1}{2}Y$右移解决虽然遇到奇数会有误差但我们只是求近似值所以可以接受

$2^{22}\left( 381\right)$十六进制为0x5400000 不过为了修正后面用了0x5f3759df

至此我们得到了平方根倒数近似值我们把它用牛顿法在近似迭代一下可以得到非常精确的值

float Q_sqrt( float number)
{
	long i;
	float x2,y;
	const float threehalfs = 1.5F;
	x2 = number * 0.5F;
	y = number;
	i = *( long * ) $y; //转换二进制类型
	i = 0x5f3759df - (i >> i );//计数平方根倒数近似值 0x5400000
	y = *( float *) &i ;//转回浮点数
	y = y *( threehalfs -(x2 * y * y));//牛顿法 求 x^(-2) -y =0 
    return y;
}

求解方程

二分法

function bisect(f::Function,a::Float64,b::Float64,tol::Float64=10^-10)
	if f(a)*f(b)>0
		error("f(a)f(b)<0 not statisfied!")
	end
		println(ceil(log(2,(b-a)/tol))+1)
	for i in 1:ceil(log(2,(b-a)/tol))+1
		global c= (a + b)*.5
		if abs(f(c)) < tol
			println("find the solution")
			return [c,f(c)]
		end
		(a,b)= f(c) < 0 ? (c,b) : (a,c)

		
	end
	return [c,f(c)]
end

不动点迭代

假设函数 $g$ 是连续可微函数, $g(r)=r, S=\left|g^{\prime}(r)\right|<1$, 则不动点迭代对于一个足够接近 $r$ 的初始估计, 以速度 $S$ 线性收玫到不动点 $r$.

function FPI(f::Function,x0::Float64,k::Int64)
	x=x0
	for i in  1:k
		x0=f(x)
		if abs(x-x0)< 10^-6	
			println("FPI convergence in ",i," step")
			return x0
		end
		x=x0
	end

	return x

end

println("FPI ",FPI(x->1-x^3.0 ,0.5,100))
println("FPI ",FPI(x->(1-x)^(1/3) ,0.5,100))
println("FPI ",FPI(x->(1+2*x^3)/(1+3*x^2) ,0.5,100))

#FPI 0.0
#FPI convergence in 39 step
#FPI 0.6823281782912355
#FPI convergence in 4 step
#FPI 0.682327803828347

牛顿法割线法

通过切线的方式不断求零点

割线法

function secant_method(f::Function,x0::Float64,x1::Float64=x0+x0*0.1,k::Int64=100,eta::Float64=10^-6)
	for i in 1:k 
		x2 = x1- f(x1)*(x1-x0)/(f(x1)- f(x0))
		if abs(x2-x1)/(x2+eta)< eta	
			println("secant_method convergence in ",i," step")
			return x2
		end
			
		

		x0 = x1
		x1 = x2
	end
	return x1
end

function Regula_Falsi(f::Function,a::Float64,b::Float64,k::Int64=100,eta::Float64=10^-6)
	if f(a)*f(b)>0
		return "error range"
	end
	for i in 1:k
		c= (b*f(a)-a*f(b))/(f(a)-f(b))
		if abs(f(c)) ==0
			return c

		end

		if check_tol(a,c,eta) <eta ||  check_tol(a,c,eta) <eta 
			println(" Regula_Falsi in  ", i, "step", "f(c)=",f(c)," c=",c )
			return c
		end
		(a,b) = f(a)*f(c) ? (a,c) : (c,a)
	end


	
end

试位法

function secant_method(f::Function,x0::Float64,x1::Float64=x0+x0*0.1,k::Int64=100,eta::Float64=10^-6)
	for i in 1:k 
		x2 = x1- f(x1)*(x1-x0)/(f(x1)- f(x0))
		if abs(x2-x1)/(x2+eta)< eta	
			println("secant_method convergence in ",i," step")
			return x2
		end
			
		

		x0 = x1
		x1 = x2
	end
	return x1
end

function Regula_Falsi(f::Function,a::Float64,b::Float64,k::Int64=100,eta::Float64=10^-6)
	if f(a)*f(b)>0
		return "error range"
	end
	for i in 1:k
		c= (b*f(a)-a*f(b))/(f(a)-f(b))
		if abs(f(c)) ==0
			return c

		end

		if check_tol(a,c,eta) <eta ||  check_tol(a,c,eta) <eta 
			println(" Regula_Falsi in  ", i, "step", "f(c)=",f(c)," c=",c )
			return c
		end
		(a,b) = f(a)*f(c) ? (a,c) : (c,a)
	end


	
end

求解方程组

首先提一句julia中 \函数可以求解线性方程组$A x = b$ ，其中$A$是方阵

当A的行数大于列数时，\函数可以计算超定方程的最小二乘解

一个坑要注意 ' 在julia里面是厄密不仅仅是转置

高斯消去法 LU分解


using LinearAlgebra
function back_substitution(a::(Array{Float64}),b::Array{Float64},order::Int64=1)
	n=length(b)
	x=zeros(1,n)
	l= order==1 ? (n:-1:1) : (1:n)
	for i in l
		l2= order==1 ? (i+1 : n) : (1:i-1)
		for j in l2
			b[i]-= a[i,j]*x[j]
		end
		
		x[i] = b[i]/a[i,i]

	end
	return x
end


function Gaussian_elimination(a::Array{Float64},b::Array{Float64})
	n=length(b)
	for j in  1 : n-1
		if abs(a[j,j]) < 10^-20
			error("zero pivot encountered")
		end

		for i in  j+1:n
			mult = a[i,j]/a[j,j]
			for k = j : n
				a[i,k] -=mult*a[j,k]  
			end
			b[i] -= mult*b[j]
		end

	end
	return back_substitution(a,b)

end

println(Gaussian_elimination(1.0*[1 2 -1; 2 1 -2 ; -3 1 1],[3 3 -6]*1.0))

function LU(m::Array{Float64})
n=length(m[1,:])

L=zeros(n,n)+I
	for j in  1 : n-1
		#L[j,j]=1
		if abs(m[j,j]) < 10^-20
			error("zero pivot encountered")
		end
		for i in  j+1:n
			
			L[i,j]=  m[i,j]/m[j,j]
			for k = j : n
				m[i,k] -= L[i,j]*m[j,k]  
			end
			
		end
	end

	return L,m

end
function LU_solve(L::Array{Float64},U::Array{Float64},b::Array{Float64})
	c=back_substitution(L,b,-1)
	println(c)
	back_substitution(U,c)
end
L,U=LU(1.0*[1 2 -1 ; 2 1 -2 ; -3 1 1])
println(L)
println(U)
println(LU_solve(L,U,[3 3 -6]*1.0))
println(back_substitution([1 0 ; 3 1 ]*1.0,[3 2]*1.0,-1))

误差来源

1.方阵 $A$的性质
\[ \text { 误差放大因子 }=\frac{\text { 相对前向误差 }}{\text { 相对后向误差 }}=\frac{\frac{\left\|x-x_a\right\|_{\infty}}{\|x\|_{\infty}}}{\frac{\|r\|_{\infty}}{\|b\|_{\infty}}} \]

方阵 $A$ 的条件数 $\operatorname{cond}(A)$ 为求解 $A x=b$ 时, 对于所有右侧向量 $b$, 可能出现的最大误差放大因子.
令人惊讶的是, 对于方阵有一个关于条件数的紧致的公式. 和向量范数类似, 定义 $n \times n$ 矩阵 $A$ 的矩阵范数为
\[ \|A\|_{\infty}= 每行元素绝对值之和的最大值 \]

即每行元素绝对值求和, 并把 $n$ 行求和的最大值作为矩阵 $A$ 的范数.
$ n n$ 矩阵 $A$ 的条件数是

\[ \operatorname{cond}(A)=\|A\| \cdot\left\|A^{-1}\right\| \]

2.淹没

高斯消去法会因为第一列主元过小而导致大数小数相减而忽略小数引入误差

PA=LU分解

可以解决0主元和淹没的问题


function PA_LU(m::Array{Float64})
	println("PA_LU")
	n=length(m[1,:])
	L=zeros(n,n)+I
	p=zeros(n,n)+I
	order=collect(1:n)
	for j in  1 : n-1
		#L[j,j]=1
		change_row=findmax(m[j:n,j])[2]
		if change_row !=j
			m[j,:],m[change_row,:]=m[change_row,:],m[j,:]
			
			order[j],order[change_row]=order[change_row],order[j]
			
		end



		if abs(m[j,j]) < 10^-20
			error("zero pivot encountered")
		end
		for i in  j+1:n
			
			L[i,j]=  m[i,j]/m[j,j]
			for k = j : n
				m[i,k] -= L[i,j]*m[j,k]  
			end
			
		end
	end
    if order!=collect(1:n)
		p=p*Permutation_matrix(order)
	end
	return p,L,m # P L U

end

P,L,U=PA_LU(1.0*[1 2 -1 ; 2 1 -2 ; -3 1 1])
println("P: "),display(P)
println("L: "),display(L)
println("U: "),display(U)
Pb=[i for i in P*[3;3;-6]]
println("Pb: "),display(Pb)
println(LU_solve(L,U,Pb*1.0))
println(LU_solve(L,U,[3.0  -6.0  3.0]))

P,L,U=PA_LU(1.0*[2 1 5 ; 4 4 -4 ; 1 3 1])
println("P: "),display(P)
println("L: "),display(L)
println("U: "),display(U)
Pb=[i for i in P*[5;0;6]]
println("Pb: "),display(Pb)
println(LU_solve(L,U,Pb*1.0))

雅可比方法

\[ \begin{split} &x_0= 初始向量 \\ &x_{k+1}=D^{-1}\left(b-(L+U) x_k\right), k=0,1,2, \cdots \end{split} \]

对于严格对角占优矩阵($\left|a_{i j}\right|>\sum\limits_{i \neq j}\left|a_{i j}\right|$) 是非奇异的且雅可比方法有唯一解

高斯-塞徳方法（Gauss-Seidel)

\[ \begin{aligned} x_0 & =\text { 初始向量 } \\ x_{k+1} & =D^{-1}\left(b-U x_k-L x_{k+1}\right), k=0,1,2, \cdots \end{aligned} \]

连续过松驰 (SOR)

\[ \begin{aligned} x_0 & =\text { 初始向量 } \\ x_{k+1} & =(\omega L+D)^{-1}\left[(1-\omega) D x_k-\omega U x_k\right]+\omega(D+\omega L)^{-1} b, k=0,1,2, \cdots \end{aligned} \]

正如雅可比和高斯-塞德尔方法, 另一种导出 SOR 的方法是将该系统看做不动点迭代. 问题 $A x=b$ 可以写成 $(L+D+U) x=b$, 乘上 $\omega$ 并重新组织方程,
\[ \begin{aligned} (\omega L+\omega D+\omega U) x & =\omega b \\ (\omega L+D) x & =\omega b-\omega U x+(1-\omega) D x \\ x & =(\omega L+D)^{-1}[(1-\omega) D x-\omega U x]+\omega(D+\omega L)^{-1} b \end{aligned} \]

function L_inv(m::Matrix{Float64},n::Int64=length(m[1,:]))

	for i in 1:n
		m[i,i]= m[i,i]!=0 ?  1/m[i,i] : 0
		for j in 1:i-1
			
				m[i,j]=-sum([ m[i,k]*m[k,j] for k in j:i-1])*m[i,i] #(1/m_ii)*sum(m_i,k*inv m_kj)
			
		end
	end
	
	return m


end
function SOR(a::Array{Float64},b::Array{Float64},w::Int64=1,k::Int64=100,eta::Float64=10^-6)
	if w<1 
		error("w must be greater than 1 ")
	end
	n=length(b)

	x=rand(n,1)

	D=zeros(n,n)
	L=zeros(n,n)
	U=zeros(n,n)
	for i in 1:n
		for j in 1:n
			if i==j
				D[i,j]=a[i,j]
			elseif i<j
				U[i,j]=a[i,j]
			else
				L[i,j]=a[i,j]
			end
		end
	end
	println(n)
	wLD_inv=
    L_inv((w*L+D),n)
	
	w_wLD_inv_b=w*wLD_inv*(b')
	wd=(1-w)*D-w*U
	display(wd)
	display(wLD_inv*(w*L+D))
	display(w_wLD_inv_b)
	display(wLD_inv)
	
	for i in 1:k
		x=wLD_inv*(wd*x)+w_wLD_inv_b
		
		r=sum([i^2 for i in (a*x-b')])^0.5
		println("r ",r)
		if  r <eta
			println("SOR convergence in ",i," step")
			return x'			
		end
	end

	return x'

end
println(SOR(1.0*[3 -1 0 0 0 0.5; -1 3 -1 0 0.5 0;0 -1 3 -1 0 0 ;0 0 -1 3 -1 0 ;0 0.5 0 -1 3 -1 ; 0.5 0 0 0 -1 3],[5/2.0 3/2.0 1.0 1.0 3/2.0 5/2.0]*1.0))

迭代法可以用来求稀疏矩阵比较有效

非线性方程组

\[ \begin{aligned} f_1(x_1, x_2,..., x_n)&=0 \\ f_2(x_1, x_2,..., x_n)&=0 \\ \vdots\\ f_m(x_1, x_2,..., x_n)&=0 \end{aligned} \]

多变量牛顿法

类似牛顿法直接用雅可比矩阵迭代求解

Broyden方法

\[ \begin{aligned} & x_0=\text { 初始向量 } \\ & A_0=\text { 初始矩阵 } \\ & \text { for } i=0,1,2, \cdots \\ & \qquad x_{i+1}=x_i-A_i^{-1} F\left(x_i\right) \\ & \qquad A_{i+1}=A_i+\frac{\left(\Delta_{i+1}-A_i \delta_{i+1}\right) \delta_{i+1}^{\mathrm{T}}}{\delta_{i+1}^{\mathrm{T}} \delta_{i+1}} \end{aligned}\notag \]
其中 $\delta_{i+1}=x_{i+1}-x_i, \Delta_{i+1}=F\left(x_{i+1}\right)-F\left(x_i\right)$

Broyden 方法II

\[ \begin{aligned} & x_0=\text { 初始向量 } \\ & B_0=\text { 初始矩阵 } \\ & \text { for } i=0,1,2, \cdots \\ & \qquad x_{i+1}=x_i-B_i F\left(x_i\right) \\ & \qquad B_{i+1}=B_i+\frac{\left(\delta_{i+1}-B_i \Delta_{i+1}\right) \delta_{i+1}^{\mathrm{T}} B_i}{\delta_{i+1}^{\mathrm{T}} B_i \Delta_{i+1}} \end{aligned}\notag \]
其中 $\delta_{i+1}=x_{i+1}-x_i, \Delta_{i+1}=F\left(x_{i+1}\right)-F\left(x_i\right)$

插值

定理 3.4 假设 $P(x)$ 是 $n-1$ 或者更低阶的插值多项式, 其拟合 $n$ 个点 $\left(x_1, y_1\right), \cdots$, $\left(x_n, y_n\right)$. 插值误差是
\[ f(x)-P(x)=\frac{\left(x-x_1\right)\left(x-x_2\right) \cdots\left(x-x_n\right)}{n !} f^{(n)}(c) \]
其中 $c$ 在最小和最大的 $n+1$ 个数字 $x, x_1, \cdots, x_n$ 之间.

拉格朗日插值

一般来说, 如果我们有 $n$ 个点 $\left(x_1, y_1\right), \cdots,\left(x_n, y_n\right)$. 对于 $1 \sim n$ 之间的每个 $k$, 定义 $n-1$ 次多项式
\[ L_k(x)=\frac{\left(x-x_1\right) \cdots\left(x-x_{k-1}\right)\left(x-x_{k+1}\right) \cdots\left(x-x_n\right)}{\left(x_k-x_1\right) \cdots\left(x_k-x_{k-1}\right)\left(x_k-x_{k+1}\right) \cdots\left(x_k-x_n\right)} \]
$L_k$ 具有一个有趣的性质, 即 $L_k\left(x_k\right)=1$, 而 $L_k\left(x_j\right)=0$, 其中 $x_j$ 是任何一个其他的数据点. 然后定义了 $n-1$ 次多项式
\[ P_{n-1}(x)=y_1 L_1(x)+\cdots+y_n L_n(x) \]
这是对 (3.1) 中多项式直接的推广, 并以相同的方式工作. 用 $x_k$ 替代 $x$ 得到
\[ P_{n-1}\left(x_k\right)=y_1 L_1\left(x_k\right)+\cdots+y_n L_n\left(x_k\right)=0+\cdots+0+y_k L_k\left(x_k\right)+0+\cdots+0=y_k \]
可以看出该多项式满足设计要求.

function Lagrange_inter(x0::Array{Float64},y0::Array{Float64})
	return x-> sum([ y0[j]*mapreduce(m->m,*,[(x-x0[i])/(x0[j]-x0[i]) for i in 1:length(y0) if i!= j ])  for j in 1:length(y0)])
end

f_Lagrange=Lagrange_inter([0.0 2.0 3],[1.0 2 4])
println(ff(1.0))

牛顿差商

\[ \begin{aligned} f\left[x_k\right] & =f\left(x_k\right) \\ f\left[x_k x_{k+1}\right] & =\frac{f\left[x_{k+1}\right]-f\left[x_k\right]}{x_{k+1}-x_k} \\ f\left[x_k x_{k+1} x_{k+2}\right] & =\frac{f\left[x_{k+1} x_{k+2}\right]-f\left[x_k x_{k+1}\right]}{x_{k+2}-x_k} \\ f\left[\begin{array}{lll} x_k x_{k+1} & x_{k+2} & x_{k+3} \end{array}\right] & =\frac{f\left[x_{k+1} x_{k+2} x_{k+3}\right]-f\left[x_k x_{k+1} x_{k+2}\right.}{x_{k+3}-x_k} \end{aligned} \]

\[ \begin{array}{l|lll} x_1 & f\left[x_1\right] & \\ & & f\left[\begin{array}{ll} x_1 x_2 \end{array}\right] & \\ x_2 & f\left[x_2\right] & & f\left[\begin{array}{lll} x_1 & x_2 & x_3 \end{array}\right] \\ & & f\left[x_2 x_3\right] & \\ x_3 & f\left[x_3\right] & \end{array} \]

function newtdd(x::Array{Float64},y::Array{Float64})
	n = length(x)
	v = zeros(n,n)
	for j in 1:n
		v[j,1]=y[j]
	end

	for i in 2:n

		for j in 1: n+1-i
			v[j,i]= (v[j+1,i-1]-v[j,i-1])/(x[j+i-1]-x[j])
		end
	end


	return(v[1,:])

切比雪夫插值节点

在区间 $[a, b]$,
\[ x_i=\frac{b+a}{2}+\frac{b-a}{2} \cos \frac{(2 i-1) \pi}{2 n} \]
$i=1, \cdots, n$ ，不等式
\[ \left|\left(x-x_1\right) \cdots\left(x-x_n\right)\right| \leqslant \frac{\left(\frac{b-a}{2}\right)^n}{2^{n-1}} \]
在区间 $[a, b]$ 上成立.

function mysin(x::Float64,n::Int64=10)
	xsin= [ pi/4 + (pi/4)* cos(i*pi/(2*n)) for i in 1:2:2*n-1]
	ysin= map(sin,xsin)
	c= newtdd(xsin,ysin)
	s=1
	x1=mod(x,2*pi)
	if x1>pi
		x1= 2*pi-x1
		s=-1
	end
	if x1 > pi/2
		x1 = pi - x1
	end

	return s*hornor(x1,n-1,c,xsin)  #hornor(x::Float64,n::Int64,c::Vector{Float64}=[1. for  i in 1:n+1],b::Vector{Float64}=[0. for  i in 1:n+1])
	end


println(sin(10)," ",mysin(10.0))

三次样条插值

为了更加精确地定义三次样条的性质, 我们做出如下的定义: 假设给定 $n$ 个点 $\left(x_1\right.$, $\left.y_1\right), \cdots,\left(x_n, y_n\right)$, 其中 $x_i$ 不同, 并且升序. 通过点 $\left(x_1, y_1\right), \cdots,\left(x_n, y_n\right)$ 的三次样条 $S(x)$ 是一组三次多项式
\[ \begin{aligned} S_1(x) & =y_1+b_1\left(x-x_1\right)+c_1\left(x-x_1\right)^2+d_1\left(x-x_1\right)^3 \text { 在区间 }\left[x_1, x_2\right] \text { 上 } \\ S_2(x) & =y_2+b_2\left(x-x_2\right)+c_2\left(x-x_2\right)^2+d_2\left(x-x_2\right)^3 \text { 在区间 }\left[x_2, x_3\right] \text { 上 } \\ & \vdots \\ S_{n-1}(x) & =y_{n-1}+b_{n-1}\left(x-x_{n-1}\right)+c_{n-1}\left(x-x_{n-1}\right)^2+d_{n-1}\left(x-x_{n-1}\right)^3 \text { 在区间 }\left[x_{n-1}, x_n\right] \text { 上 } \end{aligned} \]
并具有如下性质:
性质 1 $ S_i(x_i)=y_i, S_i(x_{i+1})=y_{i+1}$, 其中 $i=1, \cdots, n-1$.
性质 2 $ S_{i-1}^{}(x_i)=S_i{}(x_i)$, 其中 $i=2, \cdots, n-1$.
性质 3 $ S_{i-1}^{}(x_i)=S_i{}(x_i)$, 其中 $i=2, \cdots, n-1$.
性质 1 保证样条 $S(x)$ 插值数据点. 性质 2 使得相邻的样条段在它们相遇的地方斜率相同, 性质 3 则保证在两条样条段相邻的地方曲率相同, 该曲率由二阶导数表示.

这三个性质只能提供3n-5个独立方程所以要给出一个确定的样条插值需要添加额外的性质从而构成了不同的样条

自然样条

性质 4a $ S_{1}^{{}(x_1)=0,S_{n-1}}{}(x_n)=0$

使用自然样条条件 (性质 $4 a$ ) 可以得到另外的两个方程:
\[ \begin{aligned} S_1^{\prime \prime}\left(x_1\right) & =0 \rightarrow 2 c_1=0 \\ S_{n-1}^{\prime \prime}\left(x_n\right) & =0 \rightarrow 2 c_n=0 \end{aligned} \]
引人额外的未知变量 $c_n=S_{n-1}^{\prime \prime}\left(x_n\right) / 2$, 计算会更简单. 同时, 我们引人速记表示法 $\delta_i=x_{i+1}-x_i, \Delta_i=y_{i+1}-y_i$

对于 $n$ 个未知变量 $c_i$, 一共给出了 $n$ 个方程, 可以写成矩阵形式
\[ \left[\begin{array}{cccccc} 1 & 0 & 0 & & & \\ \delta_1 & 2 \delta_1+2 \delta_2 & \delta_2 & \ddots & & \\ 0 & \delta_2 & 2 \delta_2+2 \delta_3 & \delta_3 & & \\ & \ddots & \ddots & \ddots & \ddots & \\ & & & \delta_{n-2} & 2 \delta_{n-2}+2 \delta_{n-1} & \delta_{n-1} \\ & & & 0 & 0 & 1 \end{array}\right]\left[\begin{array}{c} c_1 \\ \vdots \\ c_n \end{array}\right]=\left[\begin{array}{c} 0 \\ 3\left(\frac{\Delta_2}{\delta_2}-\frac{\Delta_1}{\delta_1}\right) \\ 3\left(\frac{\Delta_{n-1}}{\delta_{n-1}}-\frac{\Delta_{n-2}}{\delta_{n-2}}\right) \\ 0 \end{array}\right] \]

function Spline_Interpolation(x::Array{Float64}, y::Array{Float64})
	n = length(x)
	
	delta=[(x[i+1] - x[i], y[i+1] -y[i]) for i in 1:n-1]
	
	m=zeros(n,n)
	m[1,1]=1
	m[n,n]=1
	b=zeros(n,1)
	for  i in 2:n-1
	
		m[i,i-1], m[i,i] ,m[i,i+1] = delta[i-1][1] , 2*delta[i-1][1] + 2*delta[i][1] ,delta[i][1]
		
		b[i]= 3* (delta[i][2]/delta[i][1] -delta[i-1][2]/delta[i-1][1])
	end
	
	
	c=m\b
	
	d=zeros(1,n)
	b=zeros(1,n)
	for i in 1:n-1
		d[i]=(c[i+1]-c[i])/(3*delta[i][1])#[((c[i+1]-c[i])/(3*delta[i][1]), delta[i][2]/delta[i][1] - delta[i][1](2* c[i] + c[i+1])/3) for i in 1:n-1]
		b[i]=delta[i][2]/delta[i][1] - delta[i][1]*(2* c[i] + c[i+1])/3
	end
return y , b , c ,d
end




println(Spline_Interpolation([0 1 2]*1.0,[3 -2 1]*1.0))

贝塞尔曲线

最小二乘法

首先提一句julia中 \函数可以求解线性方程组$A x = b$ ，其中$A$是方阵

当A的行数大于列数时，\函数可以计算超定方程的最小二乘解

法线方程

最小二乘的法线方程
对于不一致系统
\[ A x=b \]
求解 $\bar{x}$ 的计算公式如下. 我们已经知道
\[ (b-A \bar{x}) \perp\left\{A x \mid x \in R^n\right\} \]
把垂直性表示为矩阵的乘法, 我们发现对于 $R^n$ 上所有的 $x$,
\[ (A x)^{\mathrm{T}}(b-A \bar{x})=0 \]
使用前面关于转置的事实, 我们可以重写该式子, 即对于 $R^n$ 上所有的 $x$,
\[ x^{\mathrm{T}} A^{\mathrm{T}}(b-A \bar{x})=0 \]
求解
\[ A^{\mathrm{T}} A \bar{x}=A^{\mathrm{T}} b \]
得到的 $\bar{x}$, 就是最小二乘解, 它可以最小化余项 $r=b-A x$ 的欧氏长度.

function least_sqaure(a::(Array{Float64}),b::Array{Float64})

	A=a'*a
	ab=a'*b'
	P,L,U=PA_LU(A)
	Pb=[i for i in P*ab]
	return LU_solve(L,U,Pb*1.0)

end

QR 分解

经典格拉姆-施密特正交

令 $A_j(j=1, \cdots, n)$ 为线性无关向量.
$$
\[\begin{gathered} \begin{array}{l} \text { for } j=1,2, \cdots, n \\ \quad \quad y=A_j \\ \quad \quad \text { for } i=1,2, \cdots, j-1 \\ \quad \quad \quad r_{i j}=q_i^{\mathrm{T}} A_j \\ \quad \quad \quad y=y-r_{i j} q_i\\ \quad \quad \text { end } \\ %\begin{array}{l} \quad \quad y=\|y \|_2 \\ \quad \quad r_{j j}=y / r_{j j}\\ \text { end } %\end{array} \end{array} \\ \end{gathered} \notag\]

改进格拉姆-施密特正交

\[ \begin{gathered} \begin{array}{l} \text { for } j=1,2, \cdots, n \\ \quad \quad y=A_j \\ \quad \quad \text { for } i=1,2, \cdots, j-1 \\ \quad \quad \quad r_{i j}=q_i^{\mathrm{T}} y \\ \quad \quad \quad y=y-r_{i j} q_i\\ \quad \quad \text { end } \\%\begin{array}{l} \quad \quad y=\|y \|_2 \\ \quad \quad r_{j j}=y / r_{j j}\\ \text { end }%\end{array} \end{array} \\ \end{gathered} \notag \]

豪斯霍尔徳方法

豪斯霍尔德反射子令 $x$ 和 $w$ 是向量, $\|x\|_2=\|w\|_2$, 并定义 $v=w-$ $x$. 则 $H=I-2 v v^{\mathrm{T}} / v^{\mathrm{T}}$ 是对称正交矩阵, 并且 $H x=w$.

求出A矩阵每列向量的豪斯霍尔徳反射子H 得到 Q= H1 H2 H3 R= H3 H2 H1( 如果矩阵是列数为3的话)


function Householder(x::Vector{Float64},w::Vector{Float64})
	v=w-x
	P= v* (v')/(v'*v)
	return I- 2*P
end

function QR(m::Matrix{Float64},method::String="QR class")
	n=size(m)
	Q=zeros(n[1],n[2])
	Q[:,1]=m[:,1]
	r=zeros(n[1],n[2])
	println(method)
	if method =="Householder method"
		H=I+zeros(n[1],n[1])
	end

	for j in 1:n[2]
		if method == "Householder method"
			H2=I +zeros(n[1],n[1])
	 		H2[j:n[1],j:n[1]] =Householder(m[j:n[1],j],[k== j ? norm(m[j:n[1],j]) : 0.0 for k in j:n[1]])
	 		H= H*H2
	 		m=H2*m
			
		else 
			y=m[:,j]
			for i in 1:j-1

				r[i,j] =  method =="QR class" ? Q[:,i]' *m[:,j] :  Q[:,i]' *y
				y= y- r[i,j]* Q[:,i]
			end
			r[j,j]=  sqrt(sum([ i^2 for i in y ]))#norm(y)
			Q[:,j]= y/r[j,j]

		end

	end

	if method =="Householder method"
		Q,r =H , m
	end

	return Q ,r
end


Q,r = QR([1 -4 ; 2 3 ;2 2 ]*1.0,"QR class")
display(Q)
display(r)
Q,r = QR([1 -4 ; 2 3 ;2 2 ]*1.0,"QR")
display(Q)
display(r)

Q,r = QR([1 -4 ; 2 3 ;2 2 ]*1.0,"Householder method")
display(Q)
display(r)

通过 $Q R$ 分解实现最小二乘

给定 $m \times n$ 不一致系统
\[ A x=b \]
找出完全 $Q R$ 分解 $A=Q R$, 令$\hat{R}=R$ 的上 $n \times n$ 子矩阵$\hat{d}=d=Q^{\mathrm{T}} b$ 的上面的 $n$ 个元素求解 $\hat{R} \bar{x}=\hat{d}$ 得到最小二乘解 $\bar{x}$.

非线性最小二乘法

高斯-牛顿方法的重要应用是拟合具有非线性系数的模型. 令 $\left(t_1, y_1\right), \cdots,\left(t_m, y_m\right)$ 是数据点, $y=f_c(x)$ 是要进行拟合函数, 其中 $c=\left[c_1, \cdots, c_p\right]$ 是一组选择的参数, 用以最小化余项的平方和
\[ \begin{gathered} r_1(c)=f_c\left(t_1\right)-y_1 \\ \vdots \\ r_m(c)=f_c\left(t_m\right)-y_m \end{gathered} \]
的特定情况被更一般地看待, 以保证这里的特殊处理.

高斯一牛顿方法

忽略了二阶导（海森矩阵）

为了最小化
\[ r_1(x)^2+\cdots+r_m(x)^2 \notag \]
令 $x^0=$ 初始向量,

\[ \begin{aligned} \text { for } k=0,1 ， 2 \text {, } \cdots \\ A & =\operatorname{Dr}\left(x^k\right) \\ A^{\mathrm{T}} A v^k & =-A^{\mathrm{T}} r\left(x^k\right) \\ x^{k+1} & =x^k+v^k \\ \text{end} \end{aligned} \notag \]

Levenberg-Marquardt 方法

为了最小化
\[ r_1(x)^2+\cdots+r_m(x)^2 \notag \]
令 $x^0=$ 初始向量,$\lambda =$ 常数
\[ \begin{aligned} \text { for } k=0,1 ， 2 \text {, } \cdots \\ A & =\operatorname{Dr}\left(x^k\right) \\ (A^{\mathrm{T}} A + \lambda diag(A^{\mathrm{T}} A))v^k & =-A^{\mathrm{T}} r\left(x^k\right) \\ x^{k+1} & =x^k+v^k \\ \text{end} \end{aligned} \notag \]

数值微分和积分

积分

复合梯形法则

\[ \int_a^b f(x) \mathrm{d} x=\frac{h}{2}\left(y_0+y_m+2 \sum_{i=1}^{n-1} y_i\right)-\frac{(b-a) h^2}{12} f^{\prime \prime}(c) \]
其中 $h=(b-a) / m, c$ 在 $a$ 和 $b$ 之间.

推导思路对二个点进行拉格朗日插值再积分得到

复合辛普森公式

\[ \int_a^b f(x) \mathrm{d} x=\frac{h}{3}\left[y_0+y_{2 m}+4 \sum_{i=1}^m y_{2 i-1}+2 \sum_{i=1}^{m-1} y_{2 i}\right]-\frac{(b-a) h^4}{180} f^{(4)}(c) \]
其中$h=(b-a) / m $ c 在 $a$ 和 $b$ 之间.

辛普森公式推导思路对三个点进行拉格朗日插值在积分得到

复合中点法则

\[ \int_a^b f(x) \mathrm{d} x=h \sum_{i=1}^m w_{i}+\frac{(b-a) h^2}{24} f^{\prime \prime}(c) \]

其中$h=(b-a) / m $ $w_{i}$ 是在 a-b 中 m个相等子区间的中点 c 在 $a$ 和 $b$ 之间.

龙贝格积分Romberg quadrature formula

function romberg(f::Function,a::Float64,b::Float64,n::Int64,tol::Float64=10^-10)
	l=b-a
	r=zeros(Float64,n,n)
	r[1,1]=(l)*(f(a)+f(b))/2
	#r=[ (l)*(f(a)+f(b))/2]
	h= l
	@simd for j in 2:n
		h = h/2
		r[j,1]=0.5* r[j-1,1] + h * sum([ f(a+(2*i-1)*h)   for i in 1:2^(j-2)] )
		@simd for k in 2:j

			r[j,k]=(4^(k-1) * r[j,k-1] -r[j-1,k-1])/(4^(k-1)-1)
			
		  end

			if abs(r[j,j] - r[j-1,j-1]) <= tol
				println("finsh in ", j ,"step")
				return r[j,j]
			end

	end
	println("finsh in ", j ,"step" ," error :", abs(r[n,n] - r[n-1,n-1]))
	return r[n,n]
end

println(romberg(x-> log(x),1.0,2.0,100))

自适应积分

function trapezoidal_rule(f::Function,a::Float64,b::Float64)
	return  (f(a)+f(b))*(b-a)/2.0
	
end

#println(trapezoidal_rule(x->1+ sin(exp(3*x)),-1.0,1.0))
function  Simpson_rule(f::Function,a::Float64,b::Float64)
		h=(b-a)/2.0
		return  h*(f(a)+f(b)+4*f(a+h))/3.0

	end


function adapquad(f::Function,a::Float64,b::Float64,S::Function,Tol::Float64=10^-5,n::Int64=100,a_orig::Float64=a,b_orig::Float64=b)
	
	c = (a+b)/2.0
	Sab=S(f,a,b)
	Sac=S(f,a,c)
	Scb=S(f,c,b)
	get_Tol=abs(Sab-Sac-Scb)
	if  get_Tol < Tol*(b-a)/(b_orig-a_orig)
		#println("convergence  in " , n, "step")
		return  Sac + Scb
	elseif n<0
		println("no convergence get_Tol:" ,get_Tol)
		return  Sac + Scb
	else
		 return adapquad(f,a,c,S,Tol,n-1,a_orig,b_orig) + adapquad(f,c,b,S,Tol,n-1,a_orig,b_orig)

	end
end

println(adapquad(x->1.0+ sin(exp(3.0*x)),-1.0,1.0,trapezoidal_rule))

println(adapquad(x->1.0+ sin(exp(3.0*x)),-1.0,1.0,Simpson_rule))

随机数和应用

随机数生成器的目的是使输出的数字满足独立同分布。

线性同余生成器

$x_i=a x_{i-1}+b(\bmod m)$, a为乘子 b为偏移 m为模数（如 a=13,b=0,m=31，以30为周期)

Park 和 Miller 提出一中随机数生成的方案使用梅森（Mersennse) 素数 $2^p -1$ p为整数（如 a=$7^5$,b=0,m=$2^{31}-1$，重复周期理论最大值为$2^{31}-2$)

function LCG(x::Int64,n::Int64=1,a::Int64=7^5,m::Int64=2^31-1,b::Int64=0)
	result=[]

	for i in 1:n
		x= (a*x +b )%m
		push!(result,1.0*x/m)

	end 
	return result
end
#求 x^2 在0-1的积分值
println(sum([ i^2 for i in LCG(1,3000)])/3000)

生成正态分布

蒙卡1型问题

求函数的均值

蒙卡2型问题

问题无法转换为计算函数平均值的问题如求一个区域的面积

两类问题用伪随机数得到的误差

$\text { Error } \propto n^{-\frac{1}{2}}$

扩展

延迟斐波那契生成器重复周期超过$2^{1400}$

拟随机数

在条件允许下放弃对随机数的独立性要求

halton 选定一种进制如 2 进制第i随机数为如 6 在二进制下为110 倒转这个表示为 011 写成小数 0.011 然后变回原来十进制 0.375

function halton(p::Int64,n::Int64=1)
	num=Int(ceil(log(p,n)))+1
	b=zeros(num,1)
	eps=10^-6
	u=zeros(n,1)
	for j in 1:n
		
		i=1
		b[1] += 1
		while (b[i] > (p-1 +eps))
			
			b[i] = 0

			i+=1
			b[i] += 1
		end
		
		u[j] = 0
		for k in 1:num
			
			u[j] = u[j] + b[k]*(1.0*p)^(-k)
			
		end


	end
	return u
end

两类蒙卡问题用拟随机数得到的误差

1型蒙卡 $\text { Error } \propto(\ln n)^d n^{-1}$

2型蒙卡 $\text { Error } \propto(n)^{\frac{1}{2}} n^{-\frac{1}{2 d}}$

d是维度

奇异值分解

参考书籍

数值分析作者: 萨奥尔 (Timothy Sauer) 译者: 裴玉茹 / 马赓宇

浮点数 和多项式