Numerical-Analysis Homework

1.Newton iteration

python

def f(x):
    return 6*x**5-5*x**4-4*x**3+3*x**2

def df(x):
    return 30*x**4-20*x**3-12*x**2+6*x

def dx(f, x):
    return abs(0-f(x))

def newtons_method(f, df, x0, e):
    delta = dx(f, x0)
    while delta > e:
        x0 = x0 - f(x0)/df(x0)
        delta = dx(f, x0)
    print('Root is at: ', x0)
    print('f(x) at root is: ', f(x0))

x0s = [0, .5, 1]
for x0 in x0s:
    newtons_method(f, df, x0, 1e-5)

2.Power method

cpp

#include <iostream.h>
#include <math.h>
#define N 3
void matrixx(double A[N][N],double x[N],double v[N])
{
    for(int i=0; i<N; i++)
    {
        v[i]=0;
        for(int j=0; j<N; j++)
            v[i]+=A[i][j]*x[j];
    }
}
double slove(double v[N])
{
    double max;
    for(int i=0; i<N-1; i++)
        max=v[i]>v[i+1]?v[i]:v[i+1];
    return max;
}
void main()
{
//data input
    double A[N][N]= {2,-1,0,0,2,-1,0,-1,2};
    double x[N]= {0,0,1};
    double v[N]= {0,0,0};
    double u[N]= {0,0,0};
    double p[N]= {0,0,0};
    double e=1e-10,delta=1;
    int k=0;
    while(delta>=e)
    {
        for(int q=0; q<N; q++)
            p[q]=v[q];
        matrixx(A,x,v);
        for(int i=0; i<N; i++)
            u[i]=v[i]/(slove(v));
        delta=fabs(slove(v)-slove(p));
        k++;
        for(int l=0; l<N; l++)
            x[l]=u[l];
    }
    cout << "迭代次数" << k << endl;
    cout << "矩阵的特征值" << slove(v) << endl;
    cout << "(" ;
    for(int i=0; i<N; i++)
        cout << u[i] << " " ;
    cout  << ")" << endl;
}

3.Gauss-Seidel iteration

cpp

#include "stdio.h"
#include "stdlib.h"
#include "math.h"


#define N 3
#define MAXITER 10000


void main()
{
    int i, j, k;
    int iter = 0;
    int count;
//系数矩阵
    double a[N][N] = {10, -1, -2, -1, 10, -2, -1, -1, 5};
//方程右边的值
    double b[N] = {72, 83, 42};
    double new_x[N] = {0};
    double old_x[N] = {0};
    double e = 0.0000001;
    bool flag = false;
//开始迭代
    while(iter<MAXITER && !flag)
    {
        count = 0;
        iter++;
        for(i=0; i<N; i++)
            old_x[i] = new_x[i];
        for(i=0; i<N; i++)
            new_x[i] = 0;
//高斯—赛德尔计算
        for(i=0; i<N; i++)
        {
            for(j=0; j<N; j++)
            {
                if(i==j)
                    continue;
                new_x[i] -= a[i][j]*old_x[j];
            }
            new_x[i] += b[i];
            new_x[i] = new_x[i]/a[i][i];
        }
//检测是否收敛
        for(i=0; i<N; i++)
            if(fabs(new_x[i]-old_x[i])>e)
            {
                count++;
                break;
            }
        if(count==0)
            flag = true;
    }
    for(i=0; i<N; i++)
    {
        for(j=0; j<N; j++)
            printf("%lf  ", a[i][j]);
        printf("%lf\n", b[i]);
    }
    for(i=0; i<N; i++)
        printf("%lf  ", new_x[i]);
    printf("\n迭代次数：%d  \n", iter);
}