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
- y1=interp1(x,y,x1,‘metdod’);
x1取值范围不能超过x，否则会产生NaN错误
- 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
- z1=interp2(x,y,z,x1,y1,‘metdod’)
xy是已知的向量，表示采样点
z是采样点对应的函数值
- 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
- z=griddata(x,y,z,cx,cy,‘metdod’)
cx行向量cy列向量
- 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)=aex+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微分基本途径
- 利用多项式直接求导
- 直接计算差分
对向量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)'内的工具实现的
符号表达式的创建
- 直接法(式子比较简单)
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)
- 函数法sym
作用:定义单个符号变量
格式:符号变量名=sym('表达式')
说明:表达式可以是字符，字符串，数字表达式或者字符表达式。
- 函数法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)
符号表达式运算
- 四则运算
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)
- 构建符号函数 compose(复合函数)
clc
clear
close all
syms x
f=sin(x)
g=x^2
compose(f,g)
compose(g,f)
- 求反函数 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)
- 因式分解与展开 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)
- 表达式简化 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)
624
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
