首页 > 其他分享 >[数值分析]埃特肯加速法求方程解Aitken!

[数值分析]埃特肯加速法求方程解Aitken!

时间:2022-09-25 22:01:07浏览次数:52  
标签:方程解 return 迭代 Aitken theta x0 x1 def 法求

埃特肯加速法求方程解Aitken!

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return x**3 - x**2 - 1

构造迭代函数

\[\begin{aligned} x^{3} - x^{2} - 1 &= 0 \\ x = \sqrt{x^{3} - 1} \\ \end{aligned} \]

令\(\varphi(x) = \sqrt{x^{3} - 1}\)为迭代公式

def phi(x):
    return (x**3 - 1)**(1/2)
x0 = 1.5 # 初始值
while True:
    x1 = phi(x0)
    if abs(x1 - x0) < 1e-12:
        break
    x2 = phi(x1)
    theta = (x2 - x1)/(x1 - x0)
    x0 = (phi(x0) - theta * x0)/(1 - theta) 
    print(x0)

1.465365431478883
1.4655712202898679
1.4655712318767682

选取其它迭代函数

\[\begin{aligned} x^{3} - x^{2} - 1 &= 0 \\ x = {(x^{2} + 1)} ^ {\frac{1}{3}} \\ \end{aligned} \]

令\(\psi(x) = {(x^{2} + 1)} ^ {\frac{1}{3}}\)为迭代公式

# define other iteration methods
def psi(x):
    return (x**2 + 1)**(1/3)
x0 = 1.5 # 初始值
while True:
    x1 = psi(x0)
    if abs(x1 - x0) < 1e-12:
        break
    x2 = psi(x1)
    theta = (x2 - x1)/(x1 - x0)
    x0 = (psi(x0) - theta * x0)/(1 - theta) 
    print(x0)
1.4655584829667796
1.4655712318748688
1.4655712318767682

我们可以看到,\(\varphi(x)\)和 \(\psi(x)\) 的选取对迭代次数并没有明显的改善,个人认为应该是初值选取的太好了,导致迭代次数很少。下面我们来看一下不同初值的迭代次数。

选取初值

重新构造埃特肯 aitken函数

# aitken method
def aitken(x0, func = lambda x0 : x0 , eps = 1e-4 ,isprint=True):
    x0 = 1.5 # 初始值
    while True:
        x1 = func(x0)
        if abs(x1 - x0) < eps:
            break
        x2 = func(x1)
        theta = (x2 - x1)/(x1 - x0)
        x0 = (func(x0) - theta * x0)/(1 - theta) 
        print(x0)
aitken(1.5, phi, 1e-12)
1.465365431478883
1.4655712202898679
1.4655712318767682

观察不同初值的迭代次数

aitken(1, psi, 1e-12)
1.4655584829667796
1.4655712318748688
1.4655712318767682
aitken(1, phi, 1e-12)
1.465365431478883
1.4655712202898679
1.4655712318767682

更多迭代函数(free time~~~)

\(x^{3} - x^{2} - 1 = 0\)经过简单的变换可以得到下面的迭代函数

\[\begin{aligned} x &= x^{2} - \frac{1}{x} \\ x &= \frac{1}{(x - 1)^{\frac{1}{2}}} \\ x &= 1 + \frac{1}{x^{2}} \\ \end{aligned} \]

def f1(x):
    return x**2 - 1/x
def f2(x):
    return 1/((x-1)**0.5)
def f3(x):
    return 1 + 1/x**2
def f1(x):
    return x**2 - 1/x
def f2(x):
    return 1/((x-1)**0.5)
def f3(x):
    return 1 + 1/x**2
aitken(1, f1, 1e-12)
1.466725043782837
1.4655725197594736
1.465571231878372
1.465571231876768
aitken(1, f2, 1e-12)
1.4673422863283745
1.4655760852065471
1.4655712319132883
1.465571231876768
aitken(1, f3, 1e-12)
1.465858585858586
1.4655712527301703
1.4655712318767682

总结

虽然这只是数值分析的一道题但总体来说,觉得这门课如果设置成一门实践课可能也有好处,不过现在这样自己探索真的好有意思,也让我对这门课有了更深的理解。

标签:方程解,return,迭代,Aitken,theta,x0,x1,def,法求
From: https://www.cnblogs.com/Linkdom/p/16729110.html

相关文章

  • 蛮力法求矩形个数
    给定n*n的矩阵,矩阵中有0和1两个数字,现要求矩阵中只包含0的矩形的数量。枚举矩形左上角坐标......
  • 实例84 二分法求解方程
    #include<stdio.h>#include<math.h>#include<malloc.h>#include<stdlib.h>doubleFunc(double);intBisectRoot(double,double,double,double,double*,int,in......
  • 实例85 牛顿迭代法求解方程
    #include<stdio.h>#include<math.h>#include<stdlib.h>intFunction(double,double*,double*);intNewton(double*,double,int);intFunction(x,f,dy)dou......
  • 517 筛法求约数和
    视频链接: #include<iostream>usingnamespacestd;constintN=1000010;intp[N],vis[N],cnt;//g[i]表示i的最小质因子的1+p^1+...+p^kintg[N],f[N];//f[......
  • 516 筛法求约数个数
    视频链接:#include<iostream>usingnamespacestd;constintN=1000010;intp[N],vis[N],cnt;inta[N];//a[i]记录i的最小质因子的次数intd[N];//d[i]记录i......
  • 二分图最大匹配数量,匈牙利算法求解 python
    二分图最大匹配数量,匈牙利算法求解python,本质上是找增广回路"""#File:hungary.py#Time:2022/8/2821:08#Author:notomato#Description:#"""......
  • tarjan算法求强连通分量
    \(tarjan\)算法求强连通分量\(tarjan\)算法简介我在这篇博客中讲过\(tarjan\)算法的简介和求割点与桥,就不再讲述。强连通分量强连通图是指一个有向图内任意两点都能互......
  • 「学习笔记」不动点法求数列通项
    前言不动点法求数列通项是怎么回事呢?不动点法相信大家都很熟悉,但是不动点法求数列通项是怎么回事呢,下面就让小编带大家一起了解吧不动点法求数列通项,其实就是数列通项可......