matla数值与符号运算

1. 多项式计算

主函数

clc
clear
close
%多项式 x^4-12*x^3+25*x+116
p=[1 -12 0 25 116]
    %多项式求值函数 polyval 与 polyvalm 
    %格式 ：Y=polyval(p,x)
    Y=polyval(p,1)
    Y=polyval(p,2)
    Y=polyval(p,3)
    Y=polyval(p,4)
    %当x为矩阵式时 polyval 结果是 A.*A.A+5*A.*A-13*A+8*eye(size(A))
    %当x为矩阵式时 polyvalm 结果是 A*AA+5*A*A-13*A+8*ones(size(A))
    %注意x一定是方阵
    Y=polyval(p,[1,2;3,4])
    Y=polyvalm(p,[1,2;3,4])
%多项式求根函数 roots
%功能 用来求多行是=0的方程的根
%格式 r=roots（p）
%说明 r是根列向量，即r(1),r(2),r(3)....r(n)分别存储多项式p的n个根
r=roots(p)
%多项式系数函数函数 poly
%功能 已知多行是=0的方程的根，反求多项式系数
%格式 p=poly(p)
%说明 r是根列向量，即r(1),r(2),r(3)....r(n)分别存储多项式p的n个根 p是多项式系数行向量
a=poly(r)
p


    %多项式相乘 conv
    a=[1 2 3]
    b=[1 2]
    p=conv(a,b)
    %表示(x^2+2*x+3)*(x+2)=(x^3+4*x^2+7*x+6)
    %多项式相除 deconv
    [q,r]=deconv(p,b)
    %q是商多项式
    %r是余多项式
    
    
%多项式求导    polyder
a=[1 2 3]
p=polyder(a)

多项式表述方法

将多项式降幂排列，并用行向量保存。

%多项式 x^4-12*x^3+25*x+116
p=[1 -12 0 25 116]

多项式求值函数 polyval 与 polyvalm

%多项式求值函数 polyval 与 polyvalm 
    %格式 ：Y=polyval(p,x)
    Y=polyval(p,1)
    Y=polyval(p,2)
    Y=polyval(p,3)
    Y=polyval(p,4)
    %当x为矩阵式时 polyval 结果是 A.*A.A+5*A.*A-13*A+8*eye(size(A))
    %当x为矩阵式时 polyvalm 结果是 A*AA+5*A*A-13*A+8*ones(size(A))
    %注意x一定是方阵
    Y=polyval(p,[1,2;3,4])
    Y=polyvalm(p,[1,2;3,4])

多项式求根函数 roots 与 多项式系数函数函数 poly

%多项式求根函数 roots
%功能 用来求多行是=0的方程的根
%格式 r=roots（p）
%说明 r是根列向量，即r(1),r(2),r(3)....r(n)分别存储多项式p的n个根
r=roots(p)
%多项式系数函数函数 poly
%功能 已知多行是=0的方程的根，反求多项式系数
%格式 p=poly(p)
%说明 r是根列向量，即r(1),r(2),r(3)....r(n)分别存储多项式p的n个根 p是多项式系数行向量
a=poly(r)
p

多项式相乘 conv 与 多项式相除 deconv

 %多项式相乘 conv
    a=[1 2 3]
    b=[1 2]
    p=conv(a,b)
    %表示(x^2+2*x+3)*(x+2)=(x^3+4*x^2+7*x+6)
    %多项式相除 deconv
    [q,r]=deconv(p,b)
    %q是商多项式
    %r是余多项式

多项式求导 polyder

%多项式求导    polyder
a=[1 2 3]
p=polyder(a)

实验

clc
clear
close

    %多项式相乘 conv
    a=[1 2 3]
    b=[1 2]
    p=conv(a,b)
    x=p+1
    %表示(x^2+2*x+3)*(x+2)=(x^3+4*x^2+7*x+6)
    %多项式相除 deconv
    [q,r]=deconv(p,b)
    [q,r]=deconv(x,b)
    %q是商多项式
    %r是余多项式

2. 数值插值

代码

clc
clear
close
%一直概率积分函数f(x)=2/squrt(pi)*积分符号 0至x e^(-x^2)dx
    %一维数据插值 interp1
    %x1取值范围不能超过x，否则会产生NaN错误
    x=0.46:0.01:0.49;
    y=[0.484,0.494,0.503,0.512];
    x1=0.472;
    % y1=interp1(x,y,x1,'metdod')
    y1(1)=interp1(x,y,x1);
    y1(2)=interp1(x,y,x1,'metdod');
    y1(3)=interp1(x,y,x1,'linear');%用线性插值方法计算
    y1(4)=interp1(x,y,x1,'nearest');%用最近点插值方法计算
    y1(5)=interp1(x,y,x1,'puchip');%三次插值方法计算
    y1(6)=interp1(x,y,x1,'spline');
    % metdod 是插值方法 为，'linear' ‘nearest’ ‘cubic’ 改成'puchip'  ‘spline’ 
    y1
%二维数据插值 interp2 与 griddata
%类似于网格节点(均匀分布) interp2
    % z1=interp2(x,y,z,x1,y1,'metdod')
    %xy是已知的向量，表示采样点
    %z是采样点对应的函数值
    % metdod 是插值方法 与一维一致
    x=0:2.5:10
    y=[0:30:60]'
    z=[95 14 0 0 0;
        88 48 32 12 6;
        67 64 54 48 41;]
    x1=0:1:10
    y1=[0:10:60]'
    z1=interp2(x,y,z,x1,y1)
    x2=7
    y2=9
    z2=interp2(x,y,z,x2,y2)
%散乱节点(不均匀分布) griddata
    % cz=griddata(x,y,z,cx,cy,'metdod')
    %cx行向量cy列向量
    % metdod 是插值方法 为，'linear'双线性插值 ‘nearest’最邻近插值‘cubic’ 改成'puchip'双三次性插值
    %‘v4’-matlab 提供插值方法省缺时，双线性插值    
    x=-3+6*rand(200,1)
    y=-2+4*rand(200,1)
    [xx,yy]=meshgrid(x,y)
    zz=(xx.^2-2*xx).*exp(-xx.^2-yy.^2-xx.*yy)
    mesh(xx,yy,zz)
    [cx,cy]=meshgrid(-3:0.1:3,-2:0.2:2)
    figure
    cz=griddata(xx,yy,zz,cx,cy)
mesh(cx,cy,cz)

一维数据插值 interp1

1. y1=interp1(x,y,x1,‘metdod’);

x1取值范围不能超过x，否则会产生NaN错误

2. metdod 是插值方法 ‘linear’ ‘nearest’ ‘cubic’改成’puchip’ ‘spline’

linear %用线性插值方法计
nearest %用最近点插值方法计算
puchip %三次插值方法计算
spline 

代码

%一维数据插值 interp1
    %x1取值范围不能超过x，否则会产生NaN错误
    x=0.46:0.01:0.49;
    y=[0.484,0.494,0.503,0.512];
    x1=0.472;
    % y1=interp1(x,y,x1,'metdod')
    y1(1)=interp1(x,y,x1);
    y1(2)=interp1(x,y,x1,'metdod');
    y1(3)=interp1(x,y,x1,'linear');%用线性插值方法计算
    y1(4)=interp1(x,y,x1,'nearest');%用最近点插值方法计算
    y1(5)=interp1(x,y,x1,'puchip');%三次插值方法计算
    y1(6)=interp1(x,y,x1,'spline');
    % metdod 是插值方法 为，'linear' ‘nearest’ ‘cubic’ 改成'puchip'  ‘spline’ 
    y1

二维数据插值 interp2 与 griddata

类似于网格节点(均匀分布) interp2

1. z1=interp2(x,y,z,x1,y1,‘metdod’)

xy是已知的向量，表示采样点
z是采样点对应的函数值

2. metdod 是插值方法 与一维一致

代码

%二维数据插值 interp2 与 griddata
%类似于网格节点(均匀分布) interp2
    % z1=interp2(x,y,z,x1,y1,'metdod')
    %xy是已知的向量，表示采样点
    %z是采样点对应的函数值
    % metdod 是插值方法 与一维一致
    x=0:2.5:10
    y=[0:30:60]'
    z=[95 14 0 0 0;
        88 48 32 12 6;
        67 64 54 48 41;]
    x1=0:1:10
    y1=[0:10:60]'
    z1=interp2(x,y,z,x1,y1)
    x2=7
    y2=9
    z2=interp2(x,y,z,x2,y2)

散乱节点(不均匀分布) griddatac

1. z=griddata(x,y,z,cx,cy,‘metdod’)

cx行向量cy列向量

2. metdod 是插值方法

'linear'双线性插值
‘nearest’最邻近插值
‘cubic’改成'puchip'双三次性插值
‘v4’-matlab 提供插值方法省缺时，双线性插值

代码

%散乱节点(不均匀分布) griddata
    % cz=griddata(x,y,z,cx,cy,'metdod')
    %cx行向量cy列向量
    % metdod 是插值方法 为，'linear'双线性插值 ‘nearest’最邻近插值‘cubic’ 改成'puchip'双三次性插值
    %‘v4’-matlab 提供插值方法省缺时，双线性插值    
    x=-3+6*rand(200,1)
    y=-2+4*rand(200,1)
    [xx,yy]=meshgrid(x,y)
    zz=(xx.^2-2*xx).*exp(-xx.^2-yy.^2-xx.*yy)
    mesh(xx,yy,zz)
    [cx,cy]=meshgrid(-3:0.1:3,-2:0.2:2)
    figure
    cz=griddata(xx,yy,zz,cx,cy)
mesh(cx,cy,cz)

3. 数据拟合

多项式拟合 polyfit

功能

已知一组数据(采样点和对应函数值)，计算其拟合多样式

格式 [p,s]=polyfit(X,Y,m)

自变量X
因变量Y
多项式次数m
多项式系数p
参差s(误差)

说明

xy等长向量，是采样点对应函数值

例题1 用一个三次多项式在区间[0，2*pi]内逼近函数sin(x)

x=linspace(0,2*pi,50);
y=sin(x);
[p,s]=polyfit(x,y,3)
y1=polyval(p,x,'r');
plot(x,y,'*')
hold on
plot(x,y1)

例题2分别用3次多项式和6次多项式曲线拟合实验数据

x=0:0.1:1;
y=[0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2]

x=0:0.1:1;
y=[0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2]
plot(x,y,'k.','markersize',25)
p3=polyfit(x,y,3)
p6=polyfit(x,y,6)
y3=polyval(p3,x,'r');
y6=polyval(p6,x,'b');
axis([0 1.3 -2 16])
hold on
plot(x,y3)
plot(x,y6)

预判预测

x=0:0.1:1;
y=[0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2]
plot(x,y,'k.','markersize',25)
p3=polyfit(x,y,3)
p6=polyfit(x,y,6)
y3=polyval(p3,x,'r');
y6=polyval(p6,x,'b');
axis([0 1.3 -2 16])
hold on
plot(x,y3)
plot(x,y6)
figure
t=-1:0.1:2
s3=polyval(p3,t,'r');
s6=polyval(p6,t,'b');
hold on
axis([-1 2 -10 20])
plot(x,y,'k.','markersize',25)
plot(t,s3)
plot(t,s6)

当拟合曲线相差不多是，尽量选择多项式项数少的拟合曲线

x=0:0.1:1;
y=[0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2]
plot(x,y,'k.','markersize',25)
p3=polyfit(x,y,3)
p6=polyfit(x,y,6)
y3=polyval(p3,x,'r');
y6=polyval(p6,x,'b');
axis([0 1.3 -2 16])
hold on
plot(x,y3)
plot(x,y6)
figure
t=-1:0.1:2
s3=polyval(p3,t,'r');
s6=polyval(p6,t,'b');
hold on
axis([-1 2 -10 20])
plot(x,y,'k.','markersize',25)
plot(t,s3,'r-','linewidth',3)
plot(t,s6,'b-','linewidth',3)
grid on

function my_polyfit(x,y,a,b)

function   my_polyfit(x,y,a,b)
% close all
% x=0:16;
% y=[30 29.1 28.4 28.1 28.0 27.7 27.5 27.2 27.0 26.8 26.5 26.3 26.1 25.7 25.3 24.8 24.0 ];
% a=1;
% b=3;
plot(x,y,'k.','markersize',25)
pa=polyfit(x,y,a);
pb=polyfit(x,y,b);
ya=polyval(pa,x,'r');
yb=polyval(pb,x,'b');
hold on
plot(x,ya)
plot(x,yb)
figure
m=max(x);
n=min(x);
t=n-(m+n)/2:0.1:m+(m-n)/2;
sa=polyval(pa,t,'r');
sb=polyval(pb,t,'b');
hold on
plot(x,y,'k.','markersize',25)
plot(t,sa,'r-','linewidth',1)
plot(t,sb,'b-','linewidth',1)
grid on


end

非线性拟合 lsq curve fit(最小二乘数 曲线 拟合)lsqcurvefit

说明

根据给定的数据xdata，ydata(自变量，应变量)
按函数文件fun给定函数定义，以x0为初值作为最小二乘曲线拟合
返回函数fun中的系数向变量x和参差平方和resnorm
用残差平方和来评价拟合效果

例题

求参数abc使得曲线在f(x)=ae^x+b*x2+cx^3与已知数据点在最小二乘意义上充分接近

x=0:0.1:1;
y=[3.1 3.27 3.81 4.5 5.18 6 7.05 8.56 9.69 11.25 13.17];

主函数my_lsqcurvefit

clc
clear
close all
x=0:0.1:1;
y=[3.1 3.27 3.81 4.5 5.18 6 7.05 8.56 9.69 11.25 13.17];
a=[0 0 0]
[a,resnorm]=lsqcurvefit(@nihehanshu,a,x,y)%非线性拟合
plot(x,y,'k.','markersize',25)
z=nihehanshu(a,x)
hold on
plot(x,z,'r')

功能函数

function   f=nihehanshu(a,x)
f=a(1)*exp(x)+a(2)*x.^2+a(3)*x.^3;
end

4. 数据微积分

数值微分

基本原理

差分法

前向差分
后向差分
中心差分

matlab微分基本途径

1. 利用多项式直接求导
2. 直接计算差分

对向量X求前向差分：DX=diff(X)
对向量X求n阶前向差分：DX=diff(X,n)
对矩阵X求n阶前向差分：DX=diff(X,n，dim)

代码示例

clc
clear
close all
x=0:0.1:2*pi;
p=polyfit(x,sin(x),5)
dp=polyder(p)
dpx=polyval(dp,x)
dx=diff(sin([x,6.4]))/0.1
plot(x,dpx+01,'x',x,dx+2,'*',x,cos(x),'o')

数值积分

数值积分基本公式

基本梯形求积公式
    s1=(b-a)/2*[f(a)+f(b)]
基本辛普森求积公式
    s2=(b-a)/6*{f(a)+4*f[(a+b)/2]+f(b)}
复合梯形求积公式
    s3=
复合辛普森求积公式
    s4=

matlab积分函数 quad 与 quadl

quad(filename,a,b,tol,trace)
quadl(filename,a,b,tol,trace)

说明

filename 被积函数
a,b 积分上下限
tol 积分精度(默认10^(-6))
trace 0 不显示积分过程 非0 显示积分过程(默认值是零)
x1=quad('sin(x)',0,pi)
x2=quadl('sin(x)',0,pi)

trapz

trapz(x,y)%x,y是等长向量
trapz的被积函数是用表格形式来定义的.也就是x表示自变量y表示因变量，直接求积分

一重积分 quad 与 quadl 与 trapz

clc
clear
close all
x1=quad('sin(x)',0,pi)
x2=quadl('sin(x)',0,pi)
x=0:0.01:pi;
y=sin(x);
trapz(x,y)%x,y是等长向量

二重积分 dblquad 与 与(10分钟)

例题

求二重积分

主函数

clc
clear
close all
global ki;
ki=0;
z1=dblquad('fxy',-2,2,-1,1)
f=inline('exp(-x.^2/2).*sin(x.^2+y)','x','y')
z2=dblquad(f,-2,2,-1,1)

功能函数 f= fxy(x,y)

function  f= fxy(x,y)
global ki;
ki=ki+1;
f=exp(-x.^2/2).*sin(x.^2+y);
end

5. 符号运算基础

定义符号变量
定义符号矩阵
定义符号运算规则
符号计算与数值计算的转换

实验比较

clc
clear
close all
a=sym('a')
b=sym('b')
c=sym('c')
d=sym('d')
w=[a b;c d;]
det(w)

clc
clear
close all
w=10;
x=5;
y=-8;
z=11;
b=[w x;y z;]
det(b)

符号运算与数值运算区别

数值运算必须先对变量赋值，然后才能参与运算
符号运算无需对独立变量赋值，运算结果以符号形式表达
matlab的符号运算功能是通过调用 '符号运算工具箱(Symbolic Math Toolbox)'内的工具实现的

符号表达式的创建

1. 直接法(式子比较简单)

f1='a*x^2+b*x+c'%多项式 字符串构建
f2='a*x^2+b*x+c=0'%方程
f3='Dy+y^2=1'%微分方程

实例

x=sym('1+sqrt(6)')
y=cos(x)

2. 函数法sym

作用:定义单个符号变量
格式:符号变量名=sym('表达式')
说明:表达式可以是字符，字符串，数字表达式或者字符表达式。

3. 函数法syms

作用:定义单个或者多个符号变量
格式:syms var1  var2  var3  var4 ...  varn
说明:表达式可以是字符，字符串，数字表达式或者字符表达式。

例题：计算3阶范德蒙行列式的值

clc
clear
close all
syms a b c;
u=[a,b,c]
%建立范德蒙符号矩阵
t=[1,1,1;
    u;
    u.^2]
%计算矩阵t行列式值
det(t)

符号表达式运算

1. 四则运算

clc
clear
close all
syms x y z;
f=2*x+x^2*x-5*x+x^3
f=2*x/(5*x)
f=(x+y)*(x-y)

2. 构建符号函数 compose(复合函数)

clc
clear
close all
syms x
f=sin(x)
g=x^2
compose(f,g)
compose(g,f)

3. 求反函数 finverse

clc
clear
close all
syms x
f=sin(x)
g=x^2
f=finverse(f)
g=finverse(g)
m=2*x+3
n=finverse(m)
m=finverse(n)

4. 因式分解与展开 factor 与 expand

factor(s) 把多项式转为乘积模式
expand(s) 将乘积展开为多项式

clc
clear
close all
syms x a b
t=(2*a^2*b^3*x^2-4*a*b^4+10*a*b^6*x^4)/4
m=factor(t)
n=expand(m)

collect(s) 将s合并同类项
collect(s,v) 将s按变量v合并同类项(指定变量名)
    s是符号表达式或者是符号矩阵

clc
clear
close all
syms x y
s=(-7*x^2-8*y^2)*(-x^2+3*y^2)
expand(s)
collect(s,x)
collect(s,y)
factor(s)

5. 表达式简化 simplify 与 simple

simplify 用函数规则对s函数进行化简。即利用各种恒等式化简代数式
simple 调用matlab的其他函数对表达式进行综合化简，并显示化简过程

clc
clear
close all
syms x y
s=(-7*x^2-8*y^2)*(-x^2+3*y^2)
simplify(s)

符号表达时转换(符号表达是与数值表达式之间) eval

clc
clear
close all
syms x y
y=sin(x)
x=pi/4
y
f=eval(y)

clc
clear
close all
p=1.414
q=sym(p)

eval应用

clc
clear
close all
syms x y
y=x^4+x^3-x^0.1
x=2
y
eval(y)

符号矩阵 用 sym 创建

clc
clear
close all
syms a b c d
x=sym([1 2 3 4])
% x=sym('[a b c d]')%不成功
y=sym([a b c d])
z=sym([a b c 2])
w=sym([a b c 90])
x*2
y*2
z*2
w+z

符号矩阵 用 字符串 直接 创建 符号矩阵

6. 符号运算应用

符号函数 求 导 diff(f,v,n)

diff 符号函数 数值函数 两种情况下都可以运用
f 符号函数
v 符号变量
n 缺省值1 n阶导数

clc
clear
close all
syms a t
x=a*sin(t)
y=a*cos(t)
m1=diff(y)/diff(x)
m2=( diff(x)*diff(y,2)-diff(y)*diff(x,2) )/ (diff(x))^3

符号函数 求 不定积分

符号函数 求 定积分

clc
clear
close all
syms t x
f=5*x*t/(1+x^2);
m1=int(f)
m11=int(f,x)%对x求不定积分
m2=int(f,t)%对t求不定积分
m3=int(f,1,2)%对x求定积分
m4=int(f,t,1,2)%对t求定积分

符号 代数方程求解(等式)表达式(=0) solve(eq,v)

eq 是方程 而且是 等式方程
v 是求解变量
solve(eq1,eq2,eq3,eqn,v1,v2,v3,vn)
    代数方程组的求解

代数方程求解 一元解方程 solve(eq,v)

clc
clear
close all
syms y x
% y=x^2+2*x=3
y=x^2+2*x-3
solve(y)
solve(x^2+2*x-3)

代数方程求解 多元解方程 solve(eq1,eq2,eq3,eqn,v1,v2,v3,vn)

clc
clear
close all
syms y x
% 2*x+y+5=0
% x+y-4=0
m=2*x+y+5
n=x+y-4
[x,y]=solve(m,n,x,y)
solve(m,n,x,y)