MATLAB | solve函数求解析解时不支持分段函数的解决方案

Posted 2023-01-11 slandarer

MATLAB符号求解功能居然不能求分段函数？？这么离谱的事情你敢信？

离谱的问题

遇到一个很神奇的问题，这两天逛CSDN的时候发现了一个提问：

这个人在求解多元方程组的时候，遇到了以下问题，即求解时遇到了一些特殊的函数，例如分段函数就无法求解：

syms x y 
a=x+y;
if x>0
   b=1;
else
    b=2;
end
eqns = [a + b*x == 1, a - b == 2];
S=solve(eqns,[x y]);

无法从 sym 转换为 logical。

出错 demo2 (第 3 行)
if x>0

如果不是求解析解只要数值解的话vpasolve可以轻松解决，但如果非得要解析解呢？

看到这个问题我就想直接使用逻辑表达式来表示分段函数，试了一下不行：

syms x y 
a=x+y;
b=1.*(x>0)+2.*(x<=0);
eqns = [a + b*x == 1, a - b == 2];
S=solve(eqns,[x y]);

错误使用 mupadengine/feval_internal
System contains an equation of an unknown type.

出错 sym/solve (第 293 行)
sol = eng.feval_internal(‘solve’, eqns, vars, solveOptions);

出错 demo3 (第 5 行)
S=solve(eqns,[x y]);

MATLAB符号求解功能居然不能求分段函数？？这么离谱的事情你敢信？

离谱的解决方案

但是我发现，solve函数是可以在里面包含heaviside这样指示正负的函数的(就是变量取值为正数、负数和零时返回值吧不同)，于是我灵机一动把代码改为了一个很离谱的格式:

syms x y 
a=x+y;
b=floor(heaviside(x)) - 2*abs(2*heaviside(x) - 1) + 2*floor(-heaviside(x)) + 4;
b=piecewiseSym(x,0,[2,1],2);
eqns = [a + b*x == 1, a - b == 2];
S=solve(eqns,[x y]) 

S =
包含以下字段的 struct:

x: -3/2
y: 11/2

分段函数表达式生成器

但是只分成两段都得写出

b=floor(heaviside(x)) - 2abs(2heaviside(x) - 1) + 2*floor(-heaviside(x)) + 4;

这么长的表达式，那多分几段，复杂度不是要上天？这么麻烦的工作当然要封装成函数来做，于是有了下面这玩意：

工具函数

function pwFunc=piecewiseSym(x,waypoint,func,pfunc)
% @author : slandarer

gSign=[1,heaviside(x-waypoint)*2-1];
lSign=[heaviside(waypoint-x)*2-1,1];


inSign=floor((gSign+lSign)/2);
onSign=1-abs(gSign(2:end));

inFunc=inSign.*func;
onFunc=onSign.*pfunc;

pwFunc=simplify(sum(inFunc)+sum(onFunc));
end

函数介绍

输入变量含义：

• x : 自变量
• waypoint : 分段函数分段点
• func : 每一段上的函数
• pfunc : 分段点上的函数值

例子

比如对于：
$$f(x)=\{\begin{array}{cc}-x-1& x\le -1\\ -x^2+1& -1<x<1\\ (x-1{)}^3& x\ge 1\end{array}$$

可这样获取分段函数：

syms x
%            变量 分段点     每段上的函数         端点处的函数
f=piecewiseSym(x,[-1,1],[-x-1,-x^2+1,(x-1)^3],[-x-1,(x-1)^3]);

比如求该分段函数与f=0.4的交点解析解并绘图：

% 定义分段函数
syms x
%            变量 分段点     每段上的函数         端点处的函数
f=piecewiseSym(x,[-1,1],[-x-1,-x^2+1,(x-1)^3],[-x-1,(x-1)^3]);
% 求解
S=solve(f==.4,x)
% 绘图
xx=linspace(-2,2,500);
f=matlabFunction(f);
yy=f(xx);
plot(xx,yy,'LineWidth',2);

hold on
scatter(double(S),.4.*ones(length(S),1),50,'filled')

S =
-7/5
(2 ^ (1/3)*5^(2/3))/5 + 1
-15^(1/2)/5
15^(1/2)/5

前例

对于之前的例子只需要将代码改写如下形式即可：

syms x y 
a=x+y;
b=piecewiseSym(x,0,[2,1],2);
eqns = [a + b*x == 1, a - b == 2];
S=solve(eqns,[x y]) 

S =
包含以下字段的 struct:

x: -3/2
y: 11/2

---