数值分析-题目3-龙贝格和高斯求解牛顿迭代节点

-Kkking

已于 2022-10-10 18:32:24 修改

阅读量782

点赞数 1

分类专栏： matlab 文章标签： matlab

于 2022-10-10 18:16:57 首次发布

本文链接：https://blog.csdn.net/Yukinowu/article/details/127249445

版权

matlab 专栏收录该内容

1 篇文章 0 订阅

订阅专栏

题目：

1.牛顿迭代求解到x=1.64485用于参考基于龙贝格和高斯三点的迭代结果

clc
clear all
syms x;

x_0=input("请输入牛顿迭代法的初始值p_0：  ");
tol=input("请输入精度E：  ");
maxK=input("请输入最大迭代次数：  ");

[p,k,Y]=NTM(x_0,tol,maxK);

DP=sprintf("使用牛顿迭代法法迭代%d次，计算以%g为迭代初始值的解为：%g",k,x_0,p);
disp(DP);
fprintf("迭代值如下：");
disp(Y);

function [p,k,Y]=NTM(x_0,tol,maxK)
%x_0表示迭代初始值
%f表示要求解的方程
%maxK表示规定的最大迭代次数
%tolr表示允许误差
%k表示最终迭代的次数
%p表示最终迭代的值
    syms x;
    x_(1)=x_0;
    k=2;
    df=((1/(2*pi)^0.5)*exp((-1/2)*(x^2))); 
    x_(k)=x_(k-1)-subs(int(df,x,0,x_(k-1))-0.45,x,x_(k-1))/subs(df,x,x_(k-1)); %第二次迭代的结果
    while k<=maxK
        err=abs(x_(k)-x_(k-1));    %err表示相邻的迭代值的差值
        if(err<tol)
            fprintf('迭代%d次即可满足允许误差值退出\n',k-1);
            break;
        end
        k=k+1;
        x_(k)=x_(k-1)-subs(int(df,x,0,x_(k-1))-0.45,x,x_(k-1))/subs(df,x,x_(k-1));  %迭代
    end      %共迭代了k-1次
    if(k-1==maxK) 
        disp("超过最大迭代次数！");
    end
    p=x_(k); 
    k=k-1;
    Y=x_;
end

2.1龙贝格求积

function R = romberge(x,accuracy)

    fx=@(x)((1/(2*pi)^0.5)*exp((-1/2)*(x^2)));
    
    h = x;
    T = zeros(16 +1,16 +1);%生成一个16 * 16的0矩阵来存储T
    T(0 +1,0 +1) = (h/2) * (fx(0) + fx(x));%计算T(0,0)
    
    for k = 1:16 %T(1,0) --- T(16, 0)
        Sum = 0;
        deltaX = (x - 0) / 2^(k - 1); %Xk,X(k+1)之间的距离
        delta = deltaX / 2; %Xk,X(k+1/2)之间的距离
        for i = 0:(2^(k-1)-1)
            Sum = Sum + fx(0 + i * deltaX + delta);
        end
        T(k +1, 0 +1) = 0.5 * T(k-1 +1, 0 +1) + h/2 * Sum; 
        %求加速值
        for m = 1 : k
            T(k +1, m +1) = (4^m / (4^m -1)) * T(k +1, m-1 +1) - (1/(4^m - 1)) * T(k-1 +1, m-1 +1);
        end
        %判断是否到达精度
        if abs(T(k +1, m +1) - T(k-1 +1, m-1 +1)) <= accuracy
           disp('龙贝格算法积分结果:')
           result =  T(k +1, m +1) - 0.45;
           break;
        end
        h = h / 2;
    end
    
    R = result;

    disp('T表如下:')
    disp(T(1:k+1, 1:m+1));
end

2.2运用龙贝格求积计算牛顿迭代的 $gif.latex?x_%7Bk%7D$

clc
clear all
syms x;

x_0=input("请输入牛顿迭代法的初始值p_0：  ");
tol=input("请输入精度E：  ");
maxK=input("请输入最大迭代次数：  ");

[p,k,Y]=NTM(x_0,tol,maxK);

DP=sprintf("使用牛顿迭代法法迭代%d次，计算以%g为迭代初始值的解为：%g",k,x_0,p);
disp(DP);
fprintf("迭代值如下：");
disp(Y);

function [p,k,Y]=NTM(x_0,tol,maxK)
%x_0表示迭代初始值
%f表示要求解的方程
%maxK表示规定的最大迭代次数
%tolr表示允许误差
%k表示最终迭代的次数
%p表示最终迭代的值
    syms x;
    x_(1)=x_0;
    k=2;
    df=((1/(2*pi)^0.5)*exp((-1/2)*(x^2))); 
    x_(k)=x_(k-1)-romberge(x_(k-1),tol)/subs(df,x,x_(k-1)); %第二次迭代的结果
    while k<=maxK
        err=abs(x_(k)-x_(k-1));    %err表示相邻的迭代值的差值
        if(err<tol)
            fprintf('迭代%d次即可满足允许误差值退出\n',k-1);
            break;
        end
        k=k+1;
        x_(k)=x_(k-1)-romberge(x_(k-1),tol)/subs(df,x,x_(k-1));  %迭代
    end      %共迭代了k-1次
    if(k-1==maxK) 
        disp("超过最大迭代次数！");
    end
    p=x_(k); 
    k=k-1;
    Y=x_;
end

3.1三点高斯积分

function I = legender(a,b)
%Gauss-Lagendre求积
%a为下限，b为上限
syms t
syms x
f=((1/(2*pi)^0.5)*exp((-1/2)*(t^2)));
f1=subs(f,t,(b+a)/2+(b-a)*t/2)*(b-a)/2;

    I=5/9*subs(f1,t,0.6^0.5)+...
      5/9*subs(f1,t,-0.6^0.5)+...
      8/9*subs(f1,t,0)-0.45;
end

3.2运用三点高斯勒让德求积计算牛顿迭代的 $gif.latex?x_%7Bk%7D$

clc
clear
syms x;

x_0=input("请输入牛顿迭代法的初始值p_0：  ");
tol=input("请输入精度E：  ");
maxK=input("请输入最大迭代次数：  ");

[p,k,Y]=NTM(x_0,tol,maxK);

DP=sprintf("使用牛顿迭代法法迭代%d次，计算以%g为迭代初始值的解为：%g",k,x_0,p);
disp(DP);
fprintf("迭代值如下：");
disp(Y);

function [p,k,Y]=NTM(x_0,tol,maxK)
%x_0表示迭代初始值
%f表示要求解的方程
%maxK表示规定的最大迭代次数
%tolr表示允许误差
%k表示最终迭代的次数
%p表示最终迭代的值
    syms x;
    x_(1)=x_0;
    k=2;
    df=((1/(2*pi)^0.5)*exp((-1/2)*(x^2))); 
    x_(k)=x_(k-1)-legender(0,x_(k-1))/subs(df,x,x_(k-1)); %第二次迭代的结果
    while k<=maxK
        err=abs(x_(k)-x_(k-1));    %err表示相邻的迭代值的差值
        if(err<tol)
            fprintf('迭代%d次即可满足允许误差值退出\n',k-1);
            break;
        end
        k=k+1;
        x_(k)=x_(k-1)-legender(0,x_(k-1))/subs(df,x,x_(k-1));  %迭代
    end      %共迭代了k-1次
    if(k-1==maxK) 
        disp("超过最大迭代次数！");
    end
    p=x_(k); 
    k=k-1;
    Y=x_;
end

4.1龙贝格结果