计算方法 [课程学习]

最新推荐文章于 2023-04-01 23:35:57 发布

Whowell

最新推荐文章于 2023-04-01 23:35:57 发布

阅读量1.7k

点赞数 1

分类专栏：程序人生

本文链接：https://blog.csdn.net/weixin_39792252/article/details/82657731

版权

程序人生专栏收录该内容

19 篇文章 0 订阅

订阅专栏

计算方法 [课程学习] 201809

有效数字的规格化形式： $x^{*} = \pm 0.a_{1}a_{2}a_{3}...a_{n} \times 10^{m}$

计算有效数字：误差限，规格化，得到m,则有效数字是小数点后第m位到最左端第一非零数字之间的数字是有效数字，位数是有效数字位数。

一元非线性方程求根

求非线性方程 $f(x) = x^{2} - 7$ 的近似解，图像如下（找出隔根区间）：

二分法计算步骤：

(1) 输入隔根区间的端点 $a,b$ 预先给定的精度 $\xi$ ;

(2)令 $x = (a+b)/2$ ，并计算 $f(a),f(x)$

(3)若 $f(a)f(x)<0$ ，则令 $b = x$ ，否则令 $a = x$ ，转向(4)；

(4)若 $b - a < \varepsilon$ ，则输出方程满足精度的根 $x$ ，结束；否则转向(2);

double bitf(double x){
    return (x*x*x - x*x - 1.0);
}

void bits(){
    printf("二分法求方程x^3 - x^2 - 1 = 0的根\nx0 = 1.5\nx = x^3 - x^2 + x - 1\n");
    double _begin_time = clock();
    double l = 0, r = 3.0, mid;
    int cas = 0;
    while(!check(l,r)){
        mid = (l+r)/2.0;
        double f = bitf(mid), fl = bitf(l), fr = bitf(r);
        if(check(f, 0)) {
            printf("%d %.6lf\n", ++cas, mid);
            break;
        } else if(f*fl<0) r = mid;
        else l = mid;
        printf("%d %.6lf\n", ++cas, mid);
    }
    double _end_time = clock();
    printf("迭代 %d 次 用时 time = %.10lf\n\n", cas, _end_time - _begin_time);
}

迭代法：

1）方程改写： $x = \varphi (x)$ ，但是这样的形式有很多种不同形式的方程，并不是每个都能通过迭代法求得函数的近似解；

满足一下两个条件的 $\varphi (x)$ 是我们要找的：

（1）对 $\forall x \in [a,b]$ ，都有 $\varphi (x) \in [a,b]$

（2） $|\varphi ^{'}(x)| < 1$

2）这里选择函数 $x = x - \tfrac{1}{8.0}*(x^{2} - 7)$ ；

3）不断迭代， $|x_{k} - x_{k-1}| < \varepsilon$ ，则可取 $x_{k}$ 作为原方程的数值近似根；

//x = x - 1/8.0 * (x*x - 7);
void Iteration1(int t, double x, double esp = 0.0005) {
    double It[10000];
    It[0] = x;
    for(int i = 1; i < t; i++) {
        It[i] = It[i-1] - 1/8.0*(It[i-1]*It[i-1] - 7.0);
        printf("第%d次迭代x*的值为%.10lf %.6lf\n", i, It[i], It[i] - It[i-1]);
        if(fabs(It[i] - It[i-1])<esp) break;
    }
}

不动点迭代的加速

松弛迭代法

松弛迭代公式如下：

$\left\{\begin{matrix} \omega_{n} = \frac{1}{1-\varphi ^{'}(x_{n})} & & \\ x_{n+1} = (1-\omega_{n})x_{n} + \omega_{n}\varphi(x_{n}) & & \end{matrix}\right.$

注：这里的 $\varphi(x)$ 是上文转换后的 $x = \varphi (x)$

double Relaxation_iterationf(double x){
    return x*x*x - x*x + x - 1;
}

double Relaxation_iterationf1(double x){
    return x*x*3.0 - 2.0*x + 1;
}


void Relaxation_iteration(){
    printf("松弛迭代迭代法求方程x^3 - x^2 - 1 = 0的根\nx0 = 1.5\nx = x^3 - x^2 + x - 1\n");
    double _begin_time = clock();
    double prex = 1.5, nextx;
    int cas = 0;
    while(1){
        double w = 1.0/(1.0 - Relaxation_iterationf1(prex));
        nextx = (1.0 - w)*prex + w*Relaxation_iterationf(prex);
        printf("%d %.6lf\n", ++cas, nextx);
        if(fabs(prex - nextx) < esp) break;
        prex = nextx;
    }
    double _end_time = clock();
    printf("迭代 %d 次 用时 time = %.10lf\n\n", cas, _end_time - _begin_time);
}

埃特金加速迭代法

$\left\{\begin{matrix} iteration \ \ \ \ \ \ \ \widetilde{x_{k+1}} = \varphi(x_{k}), \overline{x_{k+1}} = \varphi(\widetilde{x_{k+1}}) & & \\ Improve \ \ \ \ \ \ \ \ x_{k+1} = \overline{x_{k+1}} - \frac{(\overline{x_{k+1}} - \widetilde{x_{k+1}})^{2}}{\overline{x_{k+1} - 2\widetilde{x_{k+1}} + x_{k}}} & & \end{matrix}\right.$

这样构造出的加速迭代公式不含导数信息, 但它需要用两次迭代值

//函数 x^3 - x^2 - 1 + x = x
double Etkinf(double x){
    return x*x*x - x*x + x - 1.0;
}


void Etkin(){
    printf("埃特金加速迭代法求方程x^3 - x^2 - 1 = 0的根\nx0 = 1.5\nx = x^3 - x^2 + x - 1\n");
    double _begin_time = clock();
    double _x[100000], _x1[100000], prex = 1.5, nextx;
    int i = 1;
    while(1) {
        _x[i] = Etkinf(prex); 
        _x1[i] = Etkinf(_x[i]);
        nextx = _x1[i] - (_x1[i] - _x[i])*(_x1[i] - _x[i])/(_x1[i] - 2.0*_x[i] + prex);
        printf("%d %.6lf\n", i++, nextx);
         if(fabs(prex - nextx) < esp) break;
        prex = nextx;
    }
    double _end_time = clock();
    printf("迭代 %d 次 用时 time = %.10lf\n\n", i-1, _end_time - _begin_time);
}

牛顿迭代法

加速迭代公式： $x_{k+1} = x_{k} - \frac{f(x_{k})}{f^{'}(x_{k})}$

注：这里的 $f(x)$ 指的是 $f(x) = 0$ 的方程；

bool check(double x, double y){
    return fabs(x-y) < esp;
}

//函数 f(x) = x^3 - x^2 - 1
double newtonf(double x) {
    return (x*x*x - x*x - 1.0);
    //return (x*x*x - 3*x - exp(x) + 2.0);
}

double newtonf1(double x) {
    return (x*x*3.0 - 2.0*x);
    //return (3*x*x - exp(x) - 3.0);
}
/*
    x[k+1] = x[k] - f(x[k])/f1(x[k]);
*/

void newton(){
    double prex = 1.5, nextx;
    //double prex = 1.0, nextx;
    printf("newton法求方程x^3 - x^2 - 1 = 0的根\nx0 = 1.5\nf(x) = x^3 - x^2 - 1\nf1(x) = 3*x^2 - 2*x\n");
    double _begin_time = clock();
    int cas = 0;
    while(1){
        nextx = prex - newtonf(prex)/newtonf1(prex);
        printf("%d %.6lf\n", ++cas, nextx);
        if(fabs(prex - nextx) <= esp) break;
        prex = nextx;
    }
    double _end_time = clock();
    printf("迭代 %d 次 用时 time = %.10lf\n\n", cas, _end_time - _begin_time);
}

弦截法

弦截迭代加速： $x_{k+1} = x_{k} - \frac{f(x_{k})}{f(x_{k}) - f(x_{k-1})}(x_{k} - x_{k-1})$

注：这里的 $f(x)$ 指的是 $f(x) = 0$ 的方程；

//函数 f(x) = x^3 - x^2 - 1
double chordcutf(double x){
    return (x*x*x - x*x - 1.0);
}

/*
    x[k+1] = x[k] - (x[k] - x[k-1])/(f(x[k]) - d(x[k-1]))*f(x[k])
*/


void chordcut(){
    printf("弦截法求方程x^3 - x^2 - 1 = 0的根\nx0 = 1.5\nf(x) = x^3 - x^2 - 1\n");
    double _begin_time = clock();
    double x0 = 1.5, x1 = chordcutf(x0), x2;
    printf("1 %.6lf\n", x1);
    int cas = 1;
    while(1) {
        x2 = x1 - (x1 - x0)/(chordcutf(x1) - chordcutf(x0))*chordcutf(x1);
        printf("%d %.6lf\n", ++cas, x2);
        if(check(x2, x1)) break;
        x0 = x1; x1 = x2;
    } 
    double _end_time = clock();
    printf("迭代 %d 次 用时 time = %.10lf\n\n", cas, _end_time - _begin_time);
}

求非线性方程 $f(x) = x^{3} - x - 1$ 的近似解，图像如下（找出隔值区间）：

这里选择函数 $x = x - \tfrac{1}{8.0}*(x^{3} - x - 1)$

迭代法：

//x = x - 1/8.0*(x*x*x - x - 1)
void Iteration2(int t, double x, double esp = 0.0005) {
    double It[10000];
    It[0] = x;
    for(int i = 1; i < t; i++) {
        It[i] = It[i-1] - 1/8.0*(It[i-1]*It[i-1]*It[i-1] - It[i-1] - 1.0);
        printf("第%d次迭代x*的值为%.10lf %.6lf\n", i, It[i], It[i] - It[i-1]);
        if(fabs(It[i] - It[i-1])<esp) break;
    }
}

所有的方法：

#include<bits/stdc++.h>
#define ll long long
using namespace std;
double esp = 1e-6;
/*
计算方法作业20180921
*/

bool check(double x, double y){
    return fabs(x-y) < esp;
}

//函数 f(x) = x^3 - x^2 - 1
double newtonf(double x) {
    return (x*x*x - x*x - 1.0);
    //return (x*x*x - 3*x - exp(x) + 2.0);
}

double newtonf1(double x) {
    return (x*x*3.0 - 2.0*x);
    //return (3*x*x - exp(x) - 3.0);
}
/*
    x[k+1] = x[k] - f(x[k])/f1(x[k]);
*/

void newton(){
    double prex = 1.5, nextx;
    //double prex = 1.0, nextx;
    printf("newton法求方程x^3 - x^2 - 1 = 0的根\nx0 = 1.5\nf(x) = x^3 - x^2 - 1\nf1(x) = 3*x^2 - 2*x\n");
    double _begin_time = clock();
    int cas = 0;
    while(1){
        nextx = prex - newtonf(prex)/newtonf1(prex);
        printf("%d %.6lf\n", ++cas, nextx);
        if(fabs(prex - nextx) <= esp) break;
        prex = nextx;
    }
    double _end_time = clock();
    printf("迭代 %d 次 用时 time = %.10lf\n\n", cas, _end_time - _begin_time);
}


//函数 f(x) = x^3 - x^2 - 1
double chordcutf(double x){
    return (x*x*x - x*x - 1.0);
}

/*
    x[k+1] = x[k] - (x[k] - x[k-1])/(f(x[k]) - d(x[k-1]))*f(x[k])
*/


void chordcut(){
    printf("弦截法求方程x^3 - x^2 - 1 = 0的根\nx0 = 1.5\nf(x) = x^3 - x^2 - 1\n");
    double _begin_time = clock();
    double x0 = 1.5, x1 = chordcutf(x0), x2;
    printf("1 %.6lf\n", x1);
    int cas = 1;
    while(1) {
        x2 = x1 - (x1 - x0)/(chordcutf(x1) - chordcutf(x0))*chordcutf(x1);
        printf("%d %.6lf\n", ++cas, x2);
        if(check(x2, x1)) break;
        x0 = x1; x1 = x2;
    } 
    double _end_time = clock();
    printf("迭代 %d 次 用时 time = %.10lf\n\n", cas, _end_time - _begin_time);
}

//函数 x^3 - x^2 - 1 + x = x
double Etkinf(double x){
    return x*x*x - x*x + x - 1.0;
}


void Etkin(){
    printf("埃特金加速迭代法求方程x^3 - x^2 - 1 = 0的根\nx0 = 1.5\nx = x^3 - x^2 + x - 1\n");
    double _begin_time = clock();
    double _x[100000], _x1[100000], prex = 1.5, nextx;
    int i = 1;
    while(1) {
        _x[i] = Etkinf(prex); 
        _x1[i] = Etkinf(_x[i]);
        nextx = _x1[i] - (_x1[i] - _x[i])*(_x1[i] - _x[i])/(_x1[i] - 2.0*_x[i] + prex);
        printf("%d %.6lf\n", i++, nextx);
         if(fabs(prex - nextx) < esp) break;
        prex = nextx;
    }
    double _end_time = clock();
    printf("迭代 %d 次 用时 time = %.10lf\n\n", i-1, _end_time - _begin_time);
} 


double Relaxation_iterationf(double x){
    return x*x*x - x*x + x - 1;
}

double Relaxation_iterationf1(double x){
    return x*x*3.0 - 2.0*x + 1;
}


void Relaxation_iteration(){
    printf("松弛迭代迭代法求方程x^3 - x^2 - 1 = 0的根\nx0 = 1.5\nx = x^3 - x^2 + x - 1\n");
    double _begin_time = clock();
    double prex = 1.5, nextx;
    int cas = 0;
    while(1){
        double w = 1.0/(1.0 - Relaxation_iterationf1(prex));
        nextx = (1.0 - w)*prex + w*Relaxation_iterationf(prex);
        printf("%d %.6lf\n", ++cas, nextx);
        if(fabs(prex - nextx) < esp) break;
        prex = nextx;
    }
    double _end_time = clock();
    printf("迭代 %d 次 用时 time = %.10lf\n\n", cas, _end_time - _begin_time);
}


double iterationf(double x){
    return (x - (x*x*x - x*x - 1.0)/4.0);
}

void iteration(){
    printf("普通迭代法求方程x^3 - x^2 - 1 = 0的根\nx0 = 1.5\nx = x^3 - x^2 + x - 1\n");
    double _begin_time = clock();
    double prex = 1.5, nextx;
    int cas = 0;
    while(1){
        nextx = iterationf(prex);
        printf("%d %.6lf\n", ++cas, nextx);
        if(fabs(prex - nextx) < esp) break;
        prex = nextx;
    }
    double _end_time = clock();
    printf("迭代 %d 次 用时 time = %.10lf\n\n", cas, _end_time - _begin_time);
}

double bitf(double x){
    return (x*x*x - x*x - 1.0);
}

void bits(){
    printf("二分法求方程x^3 - x^2 - 1 = 0的根\nx0 = 1.5\nx = x^3 - x^2 + x - 1\n");
    double _begin_time = clock();
    double l = 0, r = 3.0, mid;
    int cas = 0;
    while(!check(l,r)){
        mid = (l+r)/2.0;
        double f = bitf(mid), fl = bitf(l), fr = bitf(r);
        if(check(f, 0)) {
            printf("%d %.6lf\n", ++cas, mid);
            break;
        } else if(f*fl<0) r = mid;
        else l = mid;
        printf("%d %.6lf\n", ++cas, mid);
    }
    double _end_time = clock();
    printf("迭代 %d 次 用时 time = %.10lf\n\n", cas, _end_time - _begin_time);
}


int main()
{
    #ifndef ONLINE_JUDGE
        freopen("in.txt", "r", stdin);
        freopen("out.txt", "w", stdout);
        double _begin_time = clock();
    #endif
    newton();
    chordcut();
    Etkin();
    Relaxation_iteration();
    iteration();
    bits();
    #ifndef ONLINE_JUDGE
        double _end_time = clock();
        printf("time = %6lf\n", _end_time - _begin_time);
    #endif
    return 0;
}

Whowell

关注

1
点赞
踩
5

收藏

觉得还不错? 一键收藏
0
评论
计算方法 [课程学习]

计算方法 [课程学习] 201809 有效数字的规格化形式：计算有效数字：误差限，规格化，得到m,则有效数字是小数点后第m位到最左端第一非零数字之间的数字是有效数字，位数是有效数字位数。一元非线性方程求根求非线性方程的近似解，图像如下（找出隔根区间）： ...
复制链接

扫一扫