MATLAB符号计算总结

1.符号计算目录

1.引导

2.基本符号对象的创建
2.1符号数

 所用函数 :
sym(num):sc=sym(100)  %非超长符号数的建立，不超过15位
sym('num'):sc1=sym('15555555555555555555') %超长
vpa(abs(b1-b2),32) %vpa设置精度值
isAlways:
e1=sym(pi)
e2=sym('pi')
isAlways(e1==e2)
d1=sym(i)
d2=sym(j)
isAlways(d1==d2)  %千万不要sym('i')

2.2符号变量和符号数组

 所用函数：
 syms pr1 pr2 prn
 sym('a',[1,n])    sym('a',[m,1])  sym('a%d%d',[m,n])
 sym syms 没有附加选项默认复数集
 x=sym('x','SET') 
 syms x y z SET
 SET可为 integer整数集 positive rational有理数集 real实数集·
建议 assume(x,'real')这种限定性假设表述 
assumeptions咨询

3.符号运算操作基础
3.1 符号算术运算符及函数
3.2 符号关系和符号逻辑表述
4.符号表达式和符号函数
4.1 创建符号表达式和符号方程

  所用函数： 
  创建符号表达式(纯数字(子)表达式中要有sym定义的数字)
  syms a
  ex3=a*sin(1+sqrt(5))+1  %wrong
  创建符号方程
  syms x
  ex21=x/(1+sin(pi/sym(2)+1))==1 %right
  ex22=x/(1+sin(sym(pi)/2+1))==1 %right
  ex23=x/(1+sin(pi/2+sym(1)))==1  %right
  ex24=x/(sym(1)+sin(pi/2+1))==1 %符号方程的wrong代码
  解方程solve（eqin,var）命令 

4.2 创建符号函数

1.创建抽象符号函数
syms f1(x,y)
2.创建具体符号函数
syms a b
f21(x,y)=[a*sin(x);b*cos(y)];
或
f22=symfun([1,a*x+b*exp(y)],[x,y]);
3.创建符号常数函数
f31(x,y)=sym(1/2)  %直接赋值法
f32=symfun(sym(1/2),[x,y])  %symfun创建函数
4.求符号函数的函数值
数值代入法办法
f21(pi,pi)
f22(sym(0.3),-5)
f22(0.3,-5)
f31(100,200)
5.函数体及函数体分量的获取
bf2=formula(f21)
bf2(2)  正确
f21(2)  错误
6.函数自变量的获取
af1=argnames(f1)
7.符号函数规模和函数体规模
sizef21=size(f21)
sizebf2=size(bf2)
symvar(expr,n) 查找符号变量除了符号常数变量

4.3符号函数和符号表达式的比较

1.抽象函数及应用示例
syms f(x) g(x)
u=f*g
Du=diff(u,x)
syms F G
w=F*G
Dw=diff(w,x)
Dw=0
2.符号具体函数与符号表达式的内涵差别
syms a b
e2=a^2*x^2+2*a*b*x+b^2
f2(x)=e2
argnames(e2) %获取自变量，没有
argnames(f2)  %有
3.简化后各自属性不变
s4=simplify(e2) %结果仍是表达式
s4=simplify(f2) %结果仍是函数
4.函数非自变量的置换只能用subs
e5=subs(e2,[a,b],[1,2])
f5=subs(f2,{a,b},{1,2})
5.符号表达式的置换只能subs函数可以subs也可以直接输入自变量
e6=subs(e5,x,1) %仍是表达式
f61=subs(f5,x,1)%仍得到符号函数，不过是符号常数函数
f62=f5(1)  %得到函数值
6.不定积分结果
e7=int(e5,x)  %仍是表达式
f7=int(f5,x)  %函数
%ne7=e5([1,2;3,4])  wrong 只能subs置换，不能直接带
ne71=subs(e5,x,[1 2;3 4])
nf7=f5([1 2;3 4])
7.定积分结果
e8=int(e5,x,0,x)  %仍是表达式
f8=int(f5,x,0,x)   %这里变成表达式的
%%函数不定积分仍是函数，定积分变成表达式
所以要求值必须借助subs函数进行x置换
ne8=subs(e8,x,[1 2 3 4 5]/10)
nf8=subs(f8,x,[1 2 3 4 5]/10)

5.符号表达式的简化重组和子对象置换
5.1 符号表达式的简化

nexpr=simplify(expr)
nexpr=simplify(expr,name,value)
expr:输入量：符号表达式 符号函数 或他们组成的数组
nexpr 输入量形式简洁
(1)'steps'选项的应用
syms x
f1=((exp(-x*i)*i)/2-(exp(x*i)*i)/2)/(exp(-x*i)/2+exp(x*i)/2)
s11=simplify(f1)  %默认简化模式执行一次
s12=simplify(f1,'Steps',40)  %不同简化模式执行40步
（2）critertion的应用
f2=sym(1i)^(1i+1)
s21=simplify(f2)
s22=simplify(f2,'Criterion','preferReal','Steps',100) %preferReal使i尽可能在外层
（3）ignoreanalyticconstrainys的应用
syms  x y a b
f3=sqrt(x^-2)+log(x*y)+(y^a)^b+asin(sin(x))
s31=simplify(f3,'Steps',100)  %100步简化也没变化
s32=simplify(f3,'IgnoreAnalyticConstraints',true)%当x=-2就是错的了，非严格解析等价的简化

5.2符号表达式的重写

  rewrites
  syms a x t
f=a*cos(t)
rf=rewrite(f,'exp')
f1=rectangularPulse(0,0.5,t) %创建0:0.5区间高度为1的方波脉冲
rf1=rewrite(f1,'heaviside')

5.3 符号表达式的子对象置换

rs=subs(s,old,new) old 可以是符号单变量，多变量，符号表达式，符号方程、函数中的子对象，数组，混合的元胞数组

用数组置换表达式中的标量
f6=cos(a*t)*exp(-b*t)  %创建符号表达式
tt=0:0.01:2;%定义双精度数组
s61=subs(f6,{a,b,t},{6,1.2,tt}); %f6中的符号变量被双精度数值置换
%要花括号？中括号不行 既有符号变量又有双精度数值
s62=subs(f6,[a,b],[6,1.2]);  %不置换数值
ctt=class(tt)  %double对象
cs61=class(s61) %可知仍是符号对象sym
subplot(211)
fplot(s62,[0 2]); 
subplot(212)
plot(tt,s61) %可接受双精度和符号数据绘图

[nexpr,sisgma]=subexpr(expr) %对expr 自动识别将其冗子式用默认的sigma代替

[nexpr,var]=subexpr(expr,var) %，指定变量更换
pretty(expr)  

6.变精度计算及数字类型转化
6.1有限精度符号数和变精度符号数

数值计算一定存在截断误差，累计误差
符号计算不产生累积误差，但降低计算速度，所以将纯符号运算变为所谓的变精度算法，
环境精度设置命令 digits digits(n)
变精度命令vpa(x)  vpa(x,n)
(1) 两种符号数值
a=sym(1)/3  %第一种精准符号数
digits %获取符号环境默认近似精度
b=vpa(a)    %第二种有限精度符号数
whos a b
（2） 精准符号数与有限近似数的算术运算在符号环境有限精度上进行
isAlways(a==b) %32位符号环境下，结果为逻辑1
s21=a-b  %之差为0
s22=b^2
isAlways(s22==a^2)  %不相同了
（3）精准符号数与有限精度近似数的函数运算
s31=cos(a)-cos(b)
isAlways(s31==0)
s32=cos(a)/cos(b)
s33=vpa(s32)  %32为环境精度下，s33为1
isAlways(s32==sym(1)) %s32不是

6.2符号数字转换成双精度数字

1.符号数字转换为双精度数字double
数字的类别可用class判别
nd=double(num_sym)  把符号数num_sym转化为双精度数字nd
nd=double('num') 将num中的数字字符转化为ascii码
format long
nd=double(sym(1/3))
%vpa是默认转化为环境精度一般是32位的有限精度近似数

6.3 不同类型数字的相互转换

1.定义三种类型的数字
Nd=1/2+3^(2/3)  %创建双精度数，输出是浮点数
Nsym=sym(1)/2+3^(sym(2)/3) %符号数
Nstr='1/2+3^(2/3)'  %字符串
whos Nd Nsym Nstr
2.符号和字符数字转化为双精度数字double str2num str2double
N21=double(Nsym)  %符号数值转化为双精度数值
N22=str2num(Nstr) %字符串数值转化为双精度数值
N23=double(Nstr)  %ASCII码
N24=str2double(Nstr)  %该命令只能转化浮点型字符串数值，不能正确转化字符串广义有理数
3.符号或双精度数转化为字符串数char num2str int2str 
N31=char(Nsym)  %符号噢数值转化为字符串数字
N32=num2str(Nd)  %把双精度数转化为默认短格式字符串数字
N33=num2str(Nd,'%0.15f\n')  %把双精度数转化为指定格式字符串数字
N34=sprintf('%0.15f',Nd) %同上
%fprintf把输出送到显示屏或文件中，而sprintf把输出返回到一个字符串中
campare:
t=sprintf(' A circle of radius %.4g has an area of %.4g. ', rad, area);
disp(t)
 A circle of radius 2.5 has an area of 19.63.
fprintf(' A circle of radius %.4g has an area of %.4g.\n ' , rad, area)
A circle of radius 2.5 has an area of 19.63.
这里%.4g是用在函数num2str中的数据格式。%.4g就是用指数或定点标记，不管哪一种更短些，只显示至4位数字。除了g格式，还可用e (指数)和f (定点)转换。
%int2str() 双精度整数转化  num2str() 双精度浮点数转化

7.符号微积分
7.1 序列、级数的符号求和

F=symsum(f,k)  k取[0，k-1]
F=symsum(f,k,[a,b])
F=cumsum(A,dim,directon)  在dim维度，沿direction方向
f可以是符号函数，符号表达式，符号数组，符号常数

7.2 极限的求取

limit(exf)  默认变量趋于0的极限
limit（exf,x,a）双向极限
limit(exf,x,a,'right')  右极限
limit(exf,x,a,'left')  左极限

（1）单输入为符号表达式
clear all
syms w t
f1=sin(w*t)/t
F11=limit(f1)
F12=limit(f1,t,0)

(2)单输入为符号函数
f2(t)=sin(w*t)/t
F21=limit(f2)
F22=limit(f2,t,0)

(3)左右双向极限
f3=1/t
F31=limit(f3,t,0)
F32=limit(f3,t,0,'right')
F33=limit(f3,t,0,'left')

7.3 符号导数和级数展开

df=diff(f,x,n)  求f关于x的n阶导
df=diff(f,x1,x2,x3……xn) 混合导
df=jacobian(f,v)  求多元函数的jacobian矩阵
在数值计算中diff是求差分

对隐函数求导
clear
syms x y(x)  %创建基本符号对象
g=cos(x+sin(y))==sin(y)  %创建隐函数符号方程
dgdx=diff(g,x)
syms dydx
dgdx1=subs(dgdx,diff(y(x),x),dydx)
dydx=solve(dgdx1,dydx)

泰勒级数和帕德近似
st=taylor(f,var,a)  f在var=a处的5阶Taylor展开,默认3阶
st=taylor(f,var,name,value)
sp=pade(f,var,a)  三阶帕德
sp=pade(f,var,name,value)

7.4 符号积分

Fx=int(f,var)
Fx=int(f,var,[a,b])
Fx=int(f,var,[a,b],name,value)

（1）求一般定积分
syms x
f1=log(1+sqrt(x))
F1=int(f1,x,[0 1])

(2)求正态分布函数的积分
syms m a t1 t2
assume(a>0),assume(m,'real') %m>=0为实数，限定性假设
f2(x)=exp(-(x-m)^2/(2*a^2))/(a*sqrt(sym(2)*pi)) %均值为m，标准差为a，的正太概率密度函数，的符号函数表达形式
F21=int(f2,x)  %f2的原函数
F22=int(f2,x,[t1 t2]) %求f2的定积分

P21=F21(inf)-F21(-inf) % 利用原函数求无限区间定积分
P22=subs(F22,[t1,t2],[-inf inf])  %置换法求无穷区间定积分

P24=int(f2,[m-a,m+a])  %求均值两侧一个标准差区间的定积分
P25=vpa(P24)  %求有限精度数值

（3）求瑕积分
syms a
f3=1/sqrt(abs(x));
F3_0L=limit(int(f3,x,[-1 a]),a,0,'left')  %瑕点左积分存在
F3_0R=limit(int(f3,x,[a 2]),a,0,'right')   %瑕点有积分
F3=int(f3,x,[-1 2])

(4)瑕积分不存在，计算积分的柯西主值
f4=1/x
F4_0L=int(f4,x,[-1 0])
F4_0R=int(f4,x,[0 2])
F4=int(f4,x,[-1 2])
F41=int(f4,x,[-1 2],'PrincipalValue',true)

(5)不能产生封闭形式的定积分也许存在有限精度的数值定积分
f5=sin(x)/(x^2+x+1);
F51=int(f5,x,[0 5])
F52=vpa(F51)

8.符号变换及应用

Gw=fourier(gt,t,w)
gt=ifourier(Gw,w,t)

syms A t w tao
yt=A*(heaviside(t+tao/2)-heaviside(t-tao/2));
Yw0=fourier(yt,t,w)
Yw=simplify(Yw0)
Yt=ifourier(Yw,w,t)
Yts=simplify(Yt)
yt12=subs(yt,[A tao],[1 2]);
fplot(yt12,'r','LineWidth',2);
hold on
plot([-1 1;-1 1],[0 0;1 1],'w','LineWidth',1.5);%间断处曲线白化 注意矩阵数据线的对应
plot([-1 1;-1 1],[0 0;1 1],'ro','LineWidth',2); %注意矩阵数据点的对应
h0=heaviside(0);
plot([-1 1],[h0 h0],'r.','MarkerSize',25) %注意矩阵数据点的对应
axis equal
axis([-3 3 -0.5 1.5])
xlabel('t')
title(char(yt12))

Yw12=subs(Yw,[A tao],[1 2]);
figure
fplot(Yw12,'b','LineWidth',2)
axis equal
axis([-9 9 -1 3])
xlabel('w')
title(char(Yw12))

laplace变换、反变换
Gs=laplace(gt,t,s)
gt=ilaplace(Gs,s,t)

(1)对函数标量的变换反变换
clear all
syms t a s
gt1=heaviside(t-a)  %对a没有限制，
gs11=laplace(gt1,t,s) %变换失败
assume(a>0)
gs12=laplace(gt1,t,s)
gt1i=ilaplace(gs12,s,t)

（2）对函数数组的变换，反变换
syms b n x
assume(n>0);
g1(t)=dirac(t)
g2(x)=sym(1)  %常数函数
g3(t)=exp(-a*t)*sin(b*t)
g4(x)=x^n
gt2=[g1(t),g3(t);g2(x),g4(x)]  %由符号函数构建的2*2数组
tx=[t t;x x]  %根据gt2数组，构造自变量数组
gs2=laplace(gt2,tx,s)
gt2i=ilaplace(gs2,s,tx)

z变换和差分方程求解
gz=ztrans(gn,n,z)
gn=iztrans(gz,z,n)

(1)离散时域函数标量的z变换反变换
clear all
syms n m z T
gn1(n)=sin(n*T)
gz1=ztrans(gn1,n,z)
gn1i=iztrans(gz1,z,n)
(2)离散时域表达式数组的z变换反变换，不能函数
g1=kroneckerDelta(n)
g2=heaviside(m)
g3=sym(1)
gn2=[g1 g2 g3]
gz2=ztrans(gn2,[n m m],z)
gn2i=iztrans(gz2,z,[n m m])
gn2i=rewrite(gn2i,'heaviside') %重写失败

9.常微分方程的符号解法
9.1符号解法和数值解法的互补作用
9.2常微分方程ODES概述
9.3 常微分方程的符号解算命令

S=dsolve(eqs)  单输出量S承接eqs常微分方程（组）的解
S=dsolve(eqs,cond) cond条件
S=dsolve(eqs,cond,name,value)
[Y1 Y2 ……Yn]=dsolve(eqs,cond,name,value) 多输出量
syms x(t); 先定义抽象函数，
用DX(t)=diff(x,t)    D2x(t)=diff(x,t,2)
然后写成符号方程或符号表达式用作eqs输入量
x(0) Dx（0） D2x(0)当做条件

1.求解dx/dt=y,dy/dt=-x的解

syms x(t) y(t)  %创建抽象符号函数，同时指定自变量t
eq11=diff(x,t)-y  %默认等于0的表达式描述被解微分方程
eq22=diff(y,t)==-x;
cond=x(0)==sym(0.5); %只一个，所以有待定系数
S1=dsolve(eq11,eq22,cond) %单解S1是构架，它的域想，y是根据被解抽象函数名自动生成
disp(' ')
disp('微分方程组的解 S1')
disp([blanks(10),'x(t)',blanks(18),'y(t)'])
disp([S1.x,S1.y])

2.求解xy''-3y'=x^2,y(1)=0 y'(1)=-0.75
clear all
syms y(x)
Dy=diff(y,x);D2y=diff(y,x,2);
eq1=x*D2y-3*Dy==x^2    %方程注意关系符==
cond=[y(1)==0 Dy(1)==sym(-0.75)];  %注意关系符==
y1=dsolve(eq1,cond)
fplot(y1,[-1 6],'linewidth',2)
hold on
plot(1,0,'r.','MarkerSize',25) %画边值点【1 0】
hd=quiver(1,0,1,-0.75,'Color','r'); %在【1 0】点画横长1，中长0.75的红色矢量线
hd.MaxHeadSize=1;%矢量头长度为1
hold off
grid on
text(1,1,{'y(1)=0';'y{''}(1)=-0.75'})  %{''}代表一次导，{''''}两次
ylabel('y')
title(['y(x)=',char(y1)])

10.符号矩阵分析和代数方程解
10.1 符号矩阵分析

1:求解A=[a11 a12
     a21   a22]的行列式，逆，特征根，符号矩阵的输入，运算结果的准确性和繁琐，冗子式替换
     
clear all
A=sym('A%d%d',[2,2]) %自动生成带下标的符号矩阵
nA=ndims(A)  %符号矩阵A的维数
sA=size(A)  %规模
dA=det(A)  %行列式
iA=inv(A)  %逆
siA=subexpr(iA,'ID') %用冗子式ID简化表达逆矩阵
eA=eig(A)  %特征值
seA=subexpr(eA,'D')

2:G=[cost -sint  A=[3^(1/2)/2  1/2
     sint  cost]    1/2    3^(1/2)/2]的旋转作用
 %1 产生符号矩阵，和旋转变换
clear all 
 syms t
 A=sym([sqrt(3)/2 1/2;1/2 sqrt(3)/2])
 G=[cos(t) -sin(t);sin(t) cos(t)]
 GA=G*A
 
%2:旋转的动态演示
 figure
 Op=[0;0]; %原点坐标
 v1=[Op,A(:,1)]';%第一列绘制矢量用
 v2=[Op,A(:,2)]';
 Lh=plot(v1(:,1),v1(:,2),'--k',v2(:,1),v2(:,2),'b',A(:,1),A(:,2),'.r')
%黑虚表示A第一列，实蓝表示A第二列
axis([-1 1 -1 1]); axis square
Lh(1).LineWidth=4;Lh(2).LineWidth=4;Lh(3).MarkerSize=30;
Th=title(['旋转0度']);
 Gh=legend([Lh(1),Lh(2)],{'A(:,1)','A(:,2)'},'Location','southoutside','Orientation','horizontal');
 aa=(1:30)*12/sym(180)*pi;  %符号数组
 for k=1:30
   An=subs(GA,t,aa(k));  %旋转后的两个列向量
   u1=[Op,An(:,1)]'; u2=[Op,An(:,2)]';
   set(Lh(1),'XData',u1(:,1),'YData',u1(:,2))  %重置向量位置
    set(Lh(2),'XData',u2(:,1),'YData',u2(:,2))
    set(Lh(3),'XData',An(1,:),'YData',An(2,:))
    Th.String=['旋转',int2str(12*k),'度']; %动态轴图名
    if k==1;Gh.String={'An(:,1)','An(:,2)'}; end  %改图例
   pause(0.1)  %强迫更新图形，控制速度
 end

10.2 线性方程组的符号解

AX=B  %B为列向量，A行数等于B行数，列数小于等于B行数
1:X=A\B  
2:[X,r]=linsolve(A,B) 若A是方阵，r是A的倒条件数，不是方阵，r是A的秩
说明双精度除法或求逆，对矩阵元素的截断误差是铭感的

10.3 各类等式/不等式方程的符号解

S=solve(eqIn,var)
S=solve(eqIn,var,name,value)
S=vpasolve(eqIn,var,init_guess,'random',RV) 有限精度符号=解S
eqIn输入量，可以是等式不等式符号方程，也可以是表达式默认值为0
syms x
eq=x^3-3*x-6==0;
S11=solve(eq,x)

S2=vpa(S11) %隐式解的显性化处理

S3=solve(eq,x,'MaxDegree',3);
pretty(S3)

S41=solve(eq,x,'PrincipalValue',true,'MaxDegree',3) %只返回方程的一个显式主值

S42=solve(eq,x,'Real',true,'MaxDegree',3)  %只返回方程的实数解

syms y
eq1=sin(x)-cos(y)^2;  %符号表达式==0
eq2=x/2-y;
S5=solve([eq1,eq2],[x,y]) %S是架构

S5.x
S5.y

10.4 代数状态方程求符号传递函数

代数状态方程法，是为了区别求解微分状态方程法S-传函，和差分状态方程法Z-传函 三者相似
1：在结构框图上标识状态变量xi
2：建立代数状态方程，x=Ax+bU;Y=cx  根据框图得出A b c
代入G=Y/U=c(I-A)(-1)b

syms G1 G2 G3 G4 H1 H2 H3
A=[0 0 0 0 0 0 -G1;
   G2 0 0 0 0 -G2 0;
   0 G3 0 0 G3 0 0;
   0 0 G4 0 0 0 0;
   0 0 0 H2 0 0 0;
   0 0 0 H1 0 0 0;
   0 0 0 H3 0 0 0];
b=[G1;0;0;0;0;0;0];
c=[0 0 0 1 0 0 0];
Y2Ua=c*((eye(size(A))-A)\b)  %左除取代求逆inv，计算传函
disp([blanks(5),'传递函数Y2Ua为'])
pretty(Y2Ua)

% 方块参数具体化时的传递函数
syms s
Sblock={100/(s+10),1/(s+1),(s+1)/(s^2+4*s+4)，(s+1)/(s+6),(2*s+12)/(s+1),(s+1)/(s+2),1};
ww=subs(Y2Ua,{G1,G2,G3,G4,H1,H2,H3},Sblock);
Y2Ub=simplify(ww)
[NN,DD]=numden(Y2Ub);%分离出分子分母多项式
NN=expand(NN); %分子多项式展开
disp('参数具体化的传递函数Y2Ub为')
pretty(NN/DD)

11.符号函数的可视化
11.1功能绘图命令汇集

1：直接利用获得的符号表达式绘图
2：据获得的符号表达式或符号数值结果，进而转化得到数值数据，再利用数值绘图命令
f ez新老命令比较
11.2 线图绘制及修饰

fplot(f) 绘制平面函数最简调用格式
fp=fplot(f,[xmin,xmax],LineSpec,Name,Value)
生成平面函数曲线对象的详尽调用格式
fplot(xt,yt)  绘制平面参数曲线最简调用格式
fp=fplot(xt,yt,[tmin,tmax],LineSpec,Name,Value)
fplot3(xt,yt,zt) 绘制空间参数曲线的最简调用格式
fp=fplot3(xt,yt,zt,[tmin,tmax],LineSpec,Name,Value)
不管什么曲线只有一个自由度，一个独立的自由变量
%符号表达式
clear
figure
syms x
f1=exp(-x)*sin(x);
Lh1=fplot(f1,'r:','LineWidth',8)
%符号函数
f2(x)=exp(-x)*sin(x);
hold on
fplot(f2,'b-','LineWidth',4)
%参数表达式
syms t
u=t;
v=exp(-t)*sin(t);
fplot(u,v,'y-','LineWidth',1)
hold off
legend('符号表达式f1','符号函数f2(x)','参数表达式 u v')
title('输入f可用三种不同形式编码')
xlabel('x');ylabel('y')

11.3画图绘制及修饰

不管什么曲面只有2个自由度，2个独立的自由变量
fsurf(fxy) 
fm=fsurf(fxy,[xmin,xmax,ymin,ymax],LinSpec,Name,Value)
fsurf(Xuv,Yuv,Zuv)
fs=fsurf(Xuv,Yuv,Zuv,[umin,umax,vmin,vmax],LinSpec,Name,Value)
可多个符号表达式构成数组，绘制多个曲面，但符号函数只能一个

使用球坐标系参变表达式绘制第三卦缺失的上半球壳
曲面的参变量函数表述，感受所绘曲面形态与该曲面对象句柄属性间的关联，借助对象属性的点调用格式重置属性，使图形改变
clear
syms s t
x(s,t)=cos(s)*cos(t);  %用参量st表述x
y(s,t)=cos(s)*sin(t);
z(s,t)=sin(s);

figure
Sh=fsurf(x,y,z,[0,pi/2,-pi/2,pi])
xlabel('x'),ylabel('y'),zlabel('z')
str=['x(s,t)=',char(x),',y(s,t)=',char(y),',z(s,t)',char(z)]
title(str)

Sh.EdgeColor='none';
camlight

11.4 符号数字和可视化应用

1.y=f(x)=1-2/(1+exp(x))  sf(x)其不定积分
clear all
syms x y real %限定xy 在实数范围
f(x)=1-2/(1+exp(x));
disp('f(x)=')
pretty(f)
sf=int(f,x)

figure
fh=fplot([f,sf],'LineWidth',3); %功能函数绘制两条曲线
fh(1).Color='r';fh(2).Color='g';
fh(2).LineStyle=':';
grid on
title('函数f（x）和其不定积分')
xlabel('x'),ylabel('y')
legend('\it f(x)','\it\int f(x) dx','Location','North')

(2)求其反函数和其不定积分
g(y)=subs(finverse(f),x,y) %求反函数g(y)
sg=int(g,y)

gf=g(f(x))
fg=f(g(y))
%函数积分和反函数积分的互补性
syms a b real
ya=f(a),yb=f(b);

sgyb_ya=sg(yb)-sg(ya);
sfb_a=sf(b)-sf(a);
ss=sfb_a+sgyb_ya-(b*f(b)-a*f(a))
sss=simplify(ss,'IgnoreAnalyticConstraints',true,'Steps',200)
%弱化约束进行200次简化
isAlways(sss==0)

%函数积分和反函数积分的互补性的图形可视化
figure
a=1;b=3;
af=a+(b-a)/20*(0:20); %[a,b]区间数值数组
Sfxx=[a,af,b,a];Sfxy=[0,f(af),0,0];
fill(Sfxx,Sfxy,'g'); %f(x)的积分区间填充绿色
hold on
ag=f(a)+(f(b)-f(a))/20*(0:20);
Sgyx=[0 g(ag) 0 0];Sgyy=[f(a) ag f(b) f(a)];
fill(Sgyx,Sgyy,'y')
fh=fplot(f,'Color','r','LineWidth',3);
plot([-5 5;0 0]',[0 0;-1 1]','k')  %绘制过原点的xy轴线，列对应一组点连成线
hold off
xlabel('x'),ylabel('y');
title(['验证：\it\int_{a}^{b} f(x)dx+\int_{Y_{a}}^{Y_{b}} g(y)dy', ...
    '=bf(b)-af(a)'])
text(0.1,-0.06,'O')
text(a+0.1,f(a),'C'),text(b+1,f(b),'D');
text(-0.2,f(a),'Y_{a}','Rotation',-90)
text(-0.2,f(b),'Y_{b}','Rotation',-90)
text(1.4,0.3,'\it\int_{a}^{b} f(x)dx')
text(0.7,0.87,'\int\_{Y_{a}}^{Y_{b}} g(y)dy','Rotation',-90)
%混合了符号绘图命令fplot和数值绘图命令fill
%反函数积分的互补法是解决隐式反函数积分的有效途径

展开点领域内二元泰勒近似误差的可视化
（1）生成8阶截断taylor级数
clear all
syms x y
F(x,y)=sin(x^2+y);
F7=taylor(F,[x,y],'Order',8)  %关于[x,y]的8阶截断泰勒展开
（2）在较大范围内原函数和8阶泰勒展开的图形比较
subplot(121)
fsurf(F,[-2 2],'MeshDensity',30) %控制区域及网格密度绘曲面
line(0,0,F(0,0),'Marker','.','MarkerSize',20,'Color','r');
%在展开点红点标记
daspect([1 1 1]) %控制三个轴的单位长度相同
view(-16,26) %控制视角，突出误差变化
camlight,camlight left %采用光照，增强立体感
xlabel('x'),ylabel('y');
title('\it F(x,y)=sin(y+x^2)');
subplot(122)
xx=linspace(-2,2,30);
yy=xx;
[X,Y]=meshgrid(xx,yy);  %生成xy平面网格点
Z=double(F7(X,Y));  %由符号函数计算Z数据，转化成双精度
IZ=abs(Z)>1;  %标注超范围数据位置
C=Z; C(IZ)=0;
surf(X,Y,Z,C)  %用C确定曲面着色,
%fsurf(F7,[-2 2],'MeshDensity',30) 
line(0,0,F7(0,0),'Marker','.','MarkerSize',20,'Color','r')
daspect([1 1 1]);
axis([-2 2 -2 2 -1 1])
view(-16,26);
camlight,camlight left %采用光照，增强立体感
xlabel('x'),ylabel('y');
title('F的8阶截断泰勒近似');
%左边fsurf,右边surf为了曲面的着色相匹配
%%在-0.2<=x<=0.2,-0.2<=y<=0.2范围内求误差曲面
DF=abs(F-F7);
figure
xx=linspace(-0.2,0.2,50);yy=xx;
[X,Y]=meshgrid(xx,yy);
Z=double(DF(X,Y));
subplot(121)
surf(X,Y,Z);
line(0,0,DF(0,0),'Marker','.','MarkerSize',20,'Color','r');
shading interp  %插值着色
colorbar  %显示色条
caxis([0 2e-9]);  %设置色轴，加强0处色变
camlight;camlight left
view(-35,31);
axis([-0.2 0.2 -0.2 0.2 -0.2e-7 1e-7]);
pbaspect([1 1 1]) %三轴等长
xlabel('x'),ylabel('y');
title('abs(f(x,y)-F7)的数值绘图')
%%下面不必了
subplot(122)
fsurf(DF);
line(0,0,DF(0,0),'Marker','.','MarkerSize',20,'Color','r');
shading interp  %插值着色
colorbar  %显示色条
caxis([0 2e-9]);  %设置色轴，加强0处色变
camlight;camlight left
view(-35,31);
axis([-0.2 0.2 -0.2 0.2 -0.2e-7 1e-7]);
pbaspect([1 1 1]) %三轴等长
xlabel('x'),ylabel('y');
title('abs(f(x,y)-F7)的符号绘图');
%fsurf太卡了，还奇怪