数值逼近高分课程设计——Runge，D1，最佳平方逼近，复化积分(matlab源代码)

最新推荐文章于 2025-03-25 14:02:59 发布

XMU_math_doctor

最新推荐文章于 2025-03-25 14:02:59 发布

阅读量1.5k

点赞数 3

分类专栏： matlab 课程设计文章标签： matlab 数学建模算法数据分析

本文链接：https://blog.youkuaiyun.com/zhulinhao/article/details/107453244

版权

课程设计同时被 2 个专栏收录

19 篇文章

订阅专栏

matlab

13 篇文章

订阅专栏

一、Runge现象
Runge1.m

clear;clc;
f=@(x)(1./(1+25.*x.^2));%被插函数
N=[4 8 12];%分段个数
a={'r','g','b'};
b={'n=4','n=8','n=12'};
X=linspace(-1,1,100);%横坐标
plot(X,f(X),'k.-','LineWidth',3);%画出原函数
hold on;
for i=1:3
    exp=0;
    n=N(i);
    x=linspace(-1,1,n+1);
    y=f(x);
    Y=interp_R(X,x,y);%进行插值
    plot(X,Y,a{i},'LineWidth',2);%画出插值函数
    xx(1:30)=linspace(-1,-0.25,30); %随机取点，来简单的估算误差
    xx(31:70)=linspace(-0.25,0.25,40);
    xx(71:100)=linspace(0.25,1,30);
    for j=1:100
        exp=exp+abs(Y(j)-f(xx(j)));
    end
    fprintf('%s的误差绝对值之和为%f\n',b{i},exp);
    hold on;
end
legend('原图','n=4','n=8','n=12','location','southeast');
title('runge现象');

interp_R.m

function [Y] = interp_R(X,x,y)
%INTERP 插值函数
n=length(x);
A=zeros(n,n);
%系数矩阵A
for ii=1:n
    temp=1;
    A(ii,1)=1;
    for jj=2:n
        temp=temp*x(ii);
        A(ii,jj)=temp;
    end
end
y=y';
%解方程求得系数
Coeff=fliplr((A\y)');
%将横坐标带入多项式求出函数值
Y=polyval(Coeff,X);
end

二、D1样条插值
D1.m

clear;clc;
x=[1 2 3 4 6];
y=log(x);
%计算差商
f1=1;%f[1,1]
f5=1/6;%f[5,5]
dy=diff(y)./diff(x);%(f[1,2],f[2,3],f[3,4],f[4,5])
d1=(f1-dy(1))/(x(1)-x(2));%f[1,1,2];
d5=(dy(end)-f5)/(x(4)-x(5));%f[4，5，5]
for i=1:3
    d2x(i)=x(i)-x(i+2);
end
figure1 = figure;
axes1 = axes('Parent',figure1);
d_mid=-diff(dy)./d2x;
%构建三弯矩方程组
Y=6*[d1;d_mid';d5];
A=zeros(5,5);
A(1,1)=2;A(1,2)=1;A(5,4)=1;A(5,5)=2;
for i=2:4
    A(i,[i-1,i,i+1])=[(x(i)-x(i-1))/(x(i+1)-x(i-1)),2,(x(i+1)-x(i))/(x(i+1)-x(i-1))];
end
M=A\Y;%求出二阶导数
%带入多项式表达
syms t
format long g
for i=1:4
    fprintf('S(%d):   ',i);
    f(i)=y(i)+(y(i+1)-y(i))/(x(i+1)-x(i))*(t-x(i))-(M(i+1)/6+M(i)/3)*(x(i+1)-x(i))*(t-x(i))+M(i)/2*(t-x(i))^2+1/6*(M(i+1)-M(i))/(x(i+1)-x(i))*(t-x(i))^3;
    p=sym2poly(f(i));
    fprintf('%dt^3+%dt^2+%dt+%d\n\n',p(1),p(2),p(3),p(4));
end
plot(x,y,'k*');
hold on;
colour={'r','g','y','m'};
for i=1:4
    X=linspace(x(i),x(i+1),100);
    plot(X,S(X,x,y,M,i),colour{i},'LineWidth',2);
    hold on;
end
legend('源数据点','S1','S2','S3','S4');
legend1 = legend(axes1,'show');
set(legend1,...
    'Position',[0.600595238095237 0.229365079365079 0.217857142857143 0.233333333333333]);
title('D1样条');

三、最佳平方逼近
pingfang.m

clear;clc
A(1, 1)=integral(@(x)ones(size(x)), 0, 1);
A(2, 2)=integral(@(x)x.^2, 0, 1);
A(1, 2)=integral(@(x)x, 0, 1);
A(2, 1)=A(1, 2);
b(1, 1)=integral(@(x)exp(x), 0, 1);
b(2, 1)=integral(@(x)x.*exp(x), 0, 1);
a=A\b;
fplot(@(x)exp(x), [0,1], '-r');
hold on
fplot(@(x)a(1)+a(2)*x, [0, 1], '--b');
legend('y=e^x', 'y=4e-10+(18-6e)x');

四、复化积分公式
fuhua.m

clear;clc
format long
for i=1:7
    n=2^i;
    h=1/n;
    T(i)=0;
    S(i)=0;
    C(i)=0;
    for j=0:(n-1)
        T(i)=T(i)+exp(j*h)+exp((j+1)*h);
        S(i)=S(i)+exp(j*h)+4*exp((j+0.5)*h)+exp((j+1)*h);
        C(i)=C(i)+7*exp(j*h)+32*exp((j+0.25)*h)+12*exp((j+0.5)*h)+32*exp((j+0.75)*h)+7*exp((j+1)*h);
    end
    T(i)=h*T(i)/2;
    S(i)=h*S(i)/6;
    C(i)=h*C(i)/90;
end