注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

highspeedlogic

FPGA/MATLAB/Simulink

 
 
 

日志

 
 

MATLAB代做|牛顿-柯特斯数值积分公式及其MATLAB的实现  

2018-04-11 22:06:49|  分类: 默认分类 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

function y=mulNewtonCotes(fun,a,b,m,n)
% 复化Newton-Cotes数值积分公式,即在每个子区间上使用Newton-Cotes公式,然后求和
% 参数说明
% fun,积分函数的句柄,必须能够接受矢量输入
% a,积分下限
% b,积分上限
% m,将区间[a,b]等分的子区间数量
% n,采用的Newton-Cotes公式的阶数,必须满足n<8,否则积分没法保证稳定性
%    (1)n=1,即复化梯形公式
%    (2)n=2,即复化辛普森公式
%    (3)n=4,即复化科特斯公式
%
% Example
% mulNewtonCotes(fun,0,2,10,4)

xk=linspace(a,b,m+1);
for i=1:m
    s(i)=NewtonCotes(fun,xk(i),xk(i+1),n);
end
y=sum(s);

牛顿-科特斯数值积分公式

function [y,Ck,Ak]=NewtonCotes(fun,a,b,n)
% y=NewtonCotes(fun,a,b,n)
% 牛顿-科特斯数值积分公式
%
% 参数说明:
% fun,积分表达式,这里有两种选择
%      (1)积分函数句柄,必须能够接受矢量输入,比如fun=@(x)sin(x).*cos(x)
%      (2)x,y坐标的离散点,第一列为x,第二列为y,必须等距,且节点的个数小于9,比如:fun=[1:8;sin(1:8)]'
% 如果fun的表采用第二种方式,那么只需要输入第一个参数即可,否则还要输入a,b,n三个参数
% a,积分下限
% b,积分上限
% n,牛顿-科特斯数公式的阶数,必须满足1<n<7,因为n>=8时不能保证公式的稳定性
%    (1)n=1,即梯形公式
%    (2)n=2,即辛普森公式
%    (3)n=4,即科特斯公式
% y,数值积分结果
% Ck,科特斯系数
% Ak,求积系数
%
% Example
fun1=@(x)sin(x);%必须可以接受矢量输入
% fun2=[0:0.1:0.5;sin(0:0.1:0.5)];%最多8个点,必须等距
% y1=NewtonCotes(fun1,0,0.5,6)
% y2==NewtonCotes(fun2)

if nargin==1
    [mm,nn]=size(fun);
    if mm>=8
        error('为了保证NewtonCotes积分的稳定性,最多只能有9个等距节点!')
    elseif nn~=2
        error('fun构成应为:第一列为x,第二列为y,并且个数为小于10的等距节点!')
    end
    xk=fun(1,:);
    fk=fun(2,:);
    a=min(xk);
    b=max(xk);
    n=mm-1;
elseif nargin==4
    % 计算积分节点xk和节点函数值fx
    xk=linspace(a,b,n+1);
    if isa(fun,'function_handle')
        fx=fun(xk);
    else
        error('fun积分函数的句柄,且必须能够接受矢量输入!')
    end
else
    error('输入参数错误,请参考函数帮助!')
end
% 计算科特斯系数
Ck=cotescoeff(n);
% 计算求积系数
Ak=(b-a)*Ck;
% 求和算积分
y=Ak*fx';

function Ck=cotescoeff(n)
% 由于科特斯系数最多7阶,为了方便我们可以直接使用,省得每次都计算
% A1=[1,1]/2
% A2=[1,4,1]/6
% A3=[1,3,3,1]/8
% A4=[7,32,12,32,1]/90
% A5=[19,75,50,50,75,19]/288
% A6=[41,216,27,272,27,216,41]/840
% A7=[751,3577,1323,2989,2989,1323,3577,751]/17280

% 当时为了体现公式,我们使用程序计算n阶科特斯系数
for i=1:n+1
    k=i-1;
    Ck(i)=(-1)^(n-k)/factorial(k)/factorial(n-k)/n*quadl(@(t)intfun(t,n,k),0,n);
end

function f=intfun(t,n,k)
% 科特斯系数中的积分表达式
f=1;
for i=[0:k-1,k+1:n]
    f=f.*(t-i);
end

  评论这张
 
阅读(32569)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018