解方程，解方程组，插值

最新推荐文章于 2024-09-11 07:32:36 发布

雯玥

最新推荐文章于 2024-09-11 07:32:36 发布

阅读量159

点赞数 1

文章标签： matlab

本文链接：https://blog.csdn.net/weixin_45166211/article/details/117631942

版权

这段代码展示了数值方法在求解方程和插值问题中的应用，包括二分法、不动点迭代法、牛顿法以及高斯消去法、雅可比迭代法和高斯-塞德尔迭代法。此外，还演示了拉格朗日插值法用于数据拟合。代码中涉及迭代收敛性检查和矩阵非奇异性的验证，以及实际数值的计算和结果展示。

摘要由CSDN通过智能技术生成

clc,clear;
%1_1二分法
a=-1; b=-0.6;%根据方程设置有根区间端点
e=1e-3;
%fx=6*x-sin(x)+5;
while abs(b-a)>=e
    fa=6*a-sin(a)+5;
    fb=6*b-sin(b)+5;
    mid=(a+b)/2;
    fmid=6*mid-sin(mid)+5;
    if fmid==0
        break;
    end
    if fmid*fa<0
        b=mid;
    else
        a=mid;
    end
end
x=mid;
Print=['1(1)二分法求解方程的根为x=',num2str(x)];
disp(Print);

clear;
%1_1不动点迭代
x0=-1; e=1e-3;
%fi_x=(sin(x)-5)/6;迭代函数
%fi1_x=cos(x)/6;迭代函数的导函数
if abs(cos(x0)/6)>=1
        disp('此迭代函数不收敛');
end %判断迭代函数的局部收敛性
x1=(sin(x0)-5)/6;
while abs(x1-x0)>=e
    if abs(cos(x1)/6)>=1
        disp('此迭代函数不收敛');
    end %判断迭代函数的局部收敛性
    x0=x1;
    x1=(sin(x0)-5)/6;
end
x=x1;
Print=['1(1)不动点迭代法求解方程的根为x=',num2str(x)];
disp(Print);

clear;
%1_1牛顿法
x0=-1; e=1e-3; numda=1;
%fx=6*x-sin(x)+5;
%f1x=6-cos(x);
x1=x0-numda*(6*x0-sin(x0)+5)/(6-cos(x0));
while abs(x1-x0)>=e
%     if abs(6*x0-sin(x0)+5)>=abs(6*x1-sin(x1)+5)
%         Print=['numda值变为',num2str(numda)];
%         disp(Print);
%         numda=numda/2;
%         x1=x0-numda*(6*x0-sin(x0)+5)/(6-cos(x0));
%         continue;
%     end
%     numda=1;
    x0=x1;
    x1=x0-numda*(6*x0-sin(x0)+5)/(6-cos(x0));
end
x=x1;
Print=['1(1)牛顿法求解方程的根为x=',num2str(x)];
disp(Print);

clc,clear;
%2_2高斯消去
n=100;
A=zeros(n,n);
for i=1:100
    A(i,i)=3;
end
for i=2:100
    A(i-1,i)=-1;
    A(i,i-1)=-1;
end
b=ones(n,1);
b(1)=2;b(n)=2;
for k=1:n-1
    if A(k,k)==0
        disp('A不是非奇异矩阵');
        return;
    end
    i=[k+1:n];
    m=A(i,k)/A(k,k);
    A(i,i)=A(i,i)-m*A(k,i);
    b(i)=b(i)-m*b(k);
end
x=zeros(n,1);
x(n)=b(n)/A(n,n);
for i=n-1:-1:1
    x(i)=(b(i)-sum(A(i,i+1:n)*x(i+1:n)))/A(i,i);
end
disp('2_2高斯消去法求解结果为');
disp(x);

clc,clear;
%2_2雅克比迭代法
n=100;
A=zeros(n,n);
for i=1:100
    A(i,i)=3;
end
for i=2:100
    A(i-1,i)=-1;
    A(i,i-1)=-1;
end
b=ones(n,1);
b(1)=2;b(n)=2;
e=1e-6;
max=200;m=0;
x0=zeros(n,1);x=zeros(n,1);
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
B=D\(L+U);f=D\b;
[C]=eig(B);
max_C=abs(C(1,1));
for i=1:n
    if abs(C(i,1))>max_C
        max_C=abs(C(n,i));
    end
end
if max_C>=1
    disp('雅克比迭代不收敛');
    return;
end
while m<=max
    for i=1:n
        if abs(A(i,i))<1e-10
            return;
        end
        x(i)=(b(i)- A(i,1:n)*x0(1:n) + A(i,i)*x0(i))/A(i,i);
    end
    if norm(x-x0)<e
        break;
    end
    m=m+1;
    x0=x;
end
% x=B*x0+f;
% m=1;
% while norm(x-x0)>=e
%     x0=x;
%     x=B*x0+f;
%     m=m+1;
%     if m>=max
%         disp('达到最大循环次数');
%         return;
%     end
% end
disp('2_2雅克比迭代法求解结果为');
disp(x);

clc,clear;
%2_2高斯赛德尔迭代法
n=100;
A=zeros(n,n);
for i=1:100
    A(i,i)=3;
end
for i=2:100
    A(i-1,i)=-1;
    A(i,i-1)=-1;
end
b=ones(n,1);
b(1)=2;b(n)=2;
e=1e-3;
max=200;m=0;
x0=zeros(n,1);x=zeros(n,1);
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
B=(D-L)\U;f=(D-L)\b;
[X,C]=eig(B);
max_C=abs(C(1,1));
for i=1:n
    if abs(C(i,i))>max_C
        max_C=abs(C(i,i));
    end
end
if max_C>=1
    disp('G-S迭代不收敛');
    return;
end
while m<=max
    for i=1:n
        if abs(A(i,i))<1e-10
            return;
        end
        if i==1
            x(i)=  (b(i)-A(i,i+1:n)* x0(i+1:n) )/A(i,i) ;
        elseif i==n
            x(i)=  (b(i)-A(i,1:i-1)* x(1:i-1))/A(i,i);
        else
            x(i)=  (b(i)-A(i,1:i-1)* x(1:i-1) - A(i,i+1:n)*x0(i+1:n))/A(i,i);
        end
    end
        if norm(x-x0)<e
            break;
        end
        m=m+1;
        x0=x;
end
    
% x=B*x0+f;
% n=1;
% while norm(x-x0)>=e
%     x0=x;
%     x=B*x0+f;
%     n=n+1;
%     if n>=max
%         disp('达到最大循环次数');
%         break;
%     end
% end
disp('2_2G-S迭代法求解结果为');
disp(x);

资料来源：雅可比迭代法的MATLAB程序 - 百度文库 (baidu.com)

x=[25,40,50,60];
y=[95,75,63,54];
xx=70;
n=length(x);
f=0;
for i=1:n
    t=1;
    for j=1:n
        if j~=i
            t=t*(xx-x(j))/(x(i)-x(j));
        end
    end
    f=f+t*y(i);
end
Print=['拉格朗日插值法求温度70摄氏度时，预期寿命为',num2str(sprintf('%.1f',f)),'千时'];
disp(Print);

% clc, clear;
% A = [1994 67.052 0 0 0 0 0 0 0 0 0
%      1995 68.008 0 0 0 0 0 0 0 0 0
%      1996 69.803 0 0 0 0 0 0 0 0 0
%      1997 72.024 0 0 0 0 0 0 0 0 0
%      1998 73.400 0 0 0 0 0 0 0 0 0
%      1999 72.063 0 0 0 0 0 0 0 0 0
%      2000 74.669 0 0 0 0 0 0 0 0 0
%      2001 74.487 0 0 0 0 0 0 0 0 0
%      2002 74.065 0 0 0 0 0 0 0 0 0
%      2003 76.777 0 0 0 0 0 0 0 0 0]
%  n=9;
% for j=3:n+2
%     for i=j-1:n+1
%         A(i,j)=( A(i,j-1)-A(i-1,j-1) )/( A(i,1)-A(i+2-j,1) );
%     end
% end
% format long g
% A
% 
% syms x;y=0;
% for i=1:n+1
%     S=1;
%     for j=i-1:-1:1
%         S=S*(x-A(j,1));
%     end
%     y=y+A(i,i+1)*S;
% end
% y
% 
% xx=2010;yy=0;
% for i=1:n+1
%     S=1;
%     for j=i-1:-1:1
%         S=S*(xx-A(j,1));
%     end
%     yy=yy+A(i,i+1)*S;
% end
% Print=['2010年油产量预计为',num2str(yy)];
% disp(Print);
% 
% x=1994:0.01:2003;
% y=vectorize(y);
% y=eval(y);
% plot(x,y);


clc,clear;
year=[1994:2003];
bbl=[67.05,68.01,69.80,72.02,73.40,72.06,74.67,74.49,74.07,76.78];
A=zeros(10,11);
n=9;
for i=1:10
    A(i,1)=year(i);
    A(i,2)=bbl(i);
end
for j=3:n+2
    for i=j-1:n+1
        A(i,j)=(A(i,j-1)-A(i-1,j-1))/(A(i,1)-A(i+2-j,1));
    end
end
format long g
A
yy=zeros(7,1);
for xx=2004:2010
    for i=1:n+1
        s=1;
        for j=i-1:-1:1
            s=s*(xx-A(j,1));
        end
        yy(xx-2003)=yy(xx-2003)+A(i,i+1)*s;
    end
end
yy
i=[2004:2010];
plot(i,yy)

牛顿前插：等距点插值法向牛顿前插值

科学网—等距点插值法向牛顿前插值matlab程序 - 殷春武的博文 (sciencenet.cn)

三次样条插值

(2条消息) 数值分析（二续）三次样条插值二类边界完整matlab代码_cugautozp的博客-CSDN博客_matlab三次样条插值自然边界

埃尔米特插值

分段三次Hermite插值Matlab实现 - 灰信网（软件开发博客聚合） (freesion.com)