遇到一个很神奇的问题,这两天逛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
输入变量含义:
比如对于:
f
(
x
)
=
{
−
x
−
1
x
≤
−
1
−
x
2
+
1
−
1
<
x
<
1
(
x
−
1
)
3
x
≥
1
f(x)=\left\{−x−1x≤−1−x2+1−1<x<1(x−1)3x≥1
可这样获取分段函数:
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