matlab：多项式微积分、数值微积分

最新推荐文章于 2024-04-19 04:48:30 发布

INawaken

最新推荐文章于 2024-04-19 04:48:30 发布

阅读量1.6k

点赞数 1

分类专栏： matlab笔记文章标签： matlab 开发语言算法

本文链接：https://blog.csdn.net/weixin_67370287/article/details/128050674

版权

matlab笔记专栏收录该内容

7 篇文章 8 订阅

订阅专栏

一、多项式微分

1.呈现多项式

2.polyval() (多项式计算

polyval(a,x);%a：多项式，x：x取值范围

>> a=[9,-5,3,7];
>> x=-2:0.01:5;
>> f=polyval(a,x);
>> plot(x,f,'LineWidth',2);
>> xlabel('x');
>> ylabel('y');
>> set(gca,'FontSize',14);

3.polyder() （多项式微分

>> p=[5 0 -2 0 1];
>> polyder(p)
ans =
20     0    -4     0
>>  polyval(polyder(p),7)%当x=7时
ans =
        6832

>> p=[5 -7 5 10];
q=[0 4 12 -3];
x=-2:0.01:1;
>> a=polyval(p,q);
>> a=conv(p,q);%卷积和多项式乘法
>> y=polyval(a,x);%多项式计算
>> d=polyder(a);%求微分
>> f=polyval(d,x);%多项式计算，微分方程
>> hold on
>> plot(x,y,'--b',x,f,'r');
>> legend('f(x)','f\prime(x)');%f’(x)

二、多项式的积分

polyint(d,k); %d：多项式，k：自定义常数项

>> p=[5 0 -2 0 1];
d= polyint(p,3)
d =
  列 1 至 5
    1.0000         0   -0.6667         0    1.0000
  列 6
    3.0000
>> polyval(d,6) %当x等于6时
ans =
        7641

三、数值的微分

1.diff()计算差分

>> x=[1 2 5 2 1];
>> diff(x)
ans =
     1     3    -3    -1

2.计算斜率

>> x=[1,2];
>> y=[5,7];
>> slope=diff(x)./diff(y);%1/2
>> slope
slope =
0.5000

例1：

>> x0=pi/2;
>> h=0.1;
>> x=[x0 x0+h];
>> y=[sin(x0) sin(x0+h)];
>> slope=diff(y)./diff(x)
slope =
   -0.0500

例2：

>> x0=pi/2;
for i=1:1:7
h(i)=0.1^i;
dx(i)=diff([x0 x0+h(i)]);
dy(i)=diff([sin(x0) sin(x0+h(i))]);
m(i)=dy(i)./dx(i);
end
>> h
h =
  列 1 至 5
    0.1000    0.0100    0.0010    0.0001    0.0000
  列 6 至 7
0.0000    0.0000
>> m
m =
  列 1 至 5
   -0.0500   -0.0050   -0.0005   -0.0000   -0.0000
  列 6 至 7
   -0.0000   -0.0000

3.某个范围的微分

>> h=0.05;
>> x=0:h:2*pi;
>> y=sin(x);
>> m=diff(y)./diff(x);
>> hold on
>> plot(x,y);
>> plot(0:h:2*pi-h,m,'-o');%-h:差分一定比原来少一位x(1:end-1)
>> hold off
>> set(gca,'YTick',-1:1:2);
>> set(gca,'XTick',0:pi/2:2*pi);
>> set(gca,'XLim',[0,2*pi]);
>> set(gca,'YLim',[-1,2]);
>> set(gca,'XTickLabel',{'0','\pi/2','\pi','3\pi/2','2\pi'});
>> h=legend('sin(x)',strcat('sin','''','(x)'));%sin’(x)同\prime
>> set(h,'FontName','Time New Roman');%设置字体
>> box on;

>> g=colormap(lines);
>> hold on;
>> for i=1:4
x=0:power(10,-i):pi;
y=sin(x);
m=diff(y)./diff(x);
plot(x(1:end-1),m,'Color',g(i,:));
end
>> hold off
>> set(gca,'Xlim',[0,pi/2]);
>> set(gca,'YLim',[0,1.2]);
>> set(gca,'YTick',0:0.2:1.2);
>> set(gca,'XTick',0:pi/4:pi/2);
>> set(gca,'XTickLabel',{'0','\pi/4','\pi/2'});
>> set(gca,'FontName','symbol');
>> set(gca,'FontName','null');
>> h=legend('h=0.1','h=0.01','h=0.001','h=0.0001');
>> set(h,'FontName','Times New Roman');
>> box on;

>> g=colormap(lines);
hold on
for i=1:3
x=0:power(10,-i):2*pi;
y=exp(-x).*sin(x.^2./2);
m=diff(y)./diff(x);
plot(x(1:end-1),m,'Color',g(i,:));
end
>> set(gca,'Xlim',[0,2*pi]);
>> set(gca,'Ylim',[0,0.2]);
>> set(gca,'Ylim',[-0.2,0.2]);
>> set(gca,'Ylim',[-0.3,0.3]);
>> set(gca,'YTick',-0.3:0.1:0.3);
>> set(gca,'XTick',0:pi/2:2*pi);
>> set(gca,'XTickLabel',{'0','\pi/2','\pi','3*\pi/2','2*\pi'});
>> legend('h=0.1','h=0.01','h=0.001');

4.二重和三重微分

二重：(dx/dy)/dx

m=diff(f(x)),/diff(x);

m2=diff(m)./diff(x);

>> x=-2:0.005:2;
y=x.^3;
m=diff(y)./diff(x);
m2=diff(m)./diff(x(1:end-1));
plot(x,y,x(1:end-1),m,x(1:end-2),m2);
>> xlabel('x','FontSize',18);
>> ylabel('y','FontSize',18);
legend('f(x)=x^3','f''(x)','f''''(x)','Location','southeast');
>> set(gca,'FontSize',18);

四、数值积分

1.数值积分三种解法

解法一、中值法则

midpoint=(x(1:end-1)+x(2:end))./2;推导

>> h=0.05;
>> x=0:h:2;
>> midpoint=(x(1:end-1)+x(2:end))./2;
>> y=4.*midpoint.^3;
>> s=sum(h*y);
>> s
s =
   15.9950
缩小步长会使答案更加精准
>> h=0.0005;
x=0:h:2;
midpoint=(x(1:end-1)+x(2:end))./2;
y=4.*midpoint.^3;
s=sum(h*y);
s
s =
   16.0000

解法二、梯形法则

>> h=0.05;
>> x=0:h:2;
>> y=4.*x.^3;
>> s=h*trapz(y);%trapz梯形数值积分(fi+fj)/2（0<=i<=n-1,1<=j<=n）
>> s
s =
   16.0100
这里的trapz相当于
trapezoid=(y(1:end-1)+y(2:end))/2;
>> s=h*sum(trapezoid);
>> s
s =
   16.0100

解法三、辛普森公式

>> h=0.05;
>> x=0:h:2;
>> y=4.*x.^3;
>> s=h/3*(y(1)+2*sum(y(3:2:end-2))+...
4*sum(y(2:2:end))+y(end));
>> s
s =
16

三种解法比较：

2.integral()求解积分

(1).辅助知识：handle句柄

function [y] =xy_plot(input,x)
y=input(x);
plot(x,y,'r--');
xlabel('x');
ylabel('function(x)');
end

>> xy_plot(@sin,0:0.01:2*pi);

（2）.integral()

>> y=@(x) 1./(x.^3-2*x-5);
>> integral(y,0,2)
ans =
   -0.4605

(2).二重积分integral12()、三重积分integral3()

>> f=@(x,y) y.*sin(x)+x.*cos(y);
>> integral2(f,pi,2*pi,0,pi)
ans =
   -9.8696

>> f=@(x,y,z) y.*sin(x)+z.*cos(y);
>> integral3(f,0,pi,0,1,-1,1)
ans =
    2.0000