提高准确性和性能优化

这个例子表明,符号数学工具箱™有助于最小化错误当解决一个非线性方程组。

这个示例使用以下符号数学工具箱功能:

计算分析梯度使用梯度
计算雅可比矩阵分析使用雅可比矩阵
象征性的结果转换成数值函数使用matlabFunction
分析可视化与fplot

我们的目标是找到最低的函数。,或者是所谓的“香蕉”功能。

$f (x) = \sum_{我 = 1}^{n / 2} 100 (x_{2 我} - - - - - - x_{2 我 - - - - - - 1}^{2})^{2} + (1 - - - - - - x_{2 我 - - - - - - 1})^{2}$

使用fplot创建一个快速可视化的。两个变量的函数。

信谊xy= 1;b = 100;f (x, y) = (ax) ^ 2 + b * (x ^ 2) ^ 2

f (x, y) =
                  ${(x - - - - - - 1)}^{2} + One hundred. {(y - - - - - - x^{2})}^{2}$

fsurf (f[2 2 1 3])视图colormap (50 [-45])飞机caxis colorbar ([0, 1000])

图包含一个坐标轴对象。坐标轴functionsurface类型的对象包含一个对象。

动机

的解决非线性系统不包括雅可比矩阵(优化工具箱)在优化工具箱™解决同样的问题通过使用fsolve(优化工具箱)函数。这个例子中显示的工作流程有两个潜在的错误来源。你需要

翻译。函数的方程和雅可比矩阵的数学形式在文本数字代码。
显式计算雅可比矩阵。对于复杂的方程,这任务是容易出错。

这个例子表明,使用数学符号来描述问题的陈述和后续步骤最小化甚至消除这些错误。

重写。函数为一个非线性方程组

首先,转换函数。f一个非线性方程组F,在那里F是一个梯度f。

清楚n = 64;x =符号(“x”[n, 1]);f =总和(100 * (x (2:2: n) - x (1:2: n)。^ 2)。^ 2 + (1 - x (1:2: n)) ^ 2);F =梯度(F (x);

显示第一个10方程:

F (1:10)

ans =
                  
                    $(\begin{array}{c} 2 x_{1} - - - - - - 400年 x_{1} (x_{2} - - - - - - {x_{1}}^{2}) - - - - - - 2 \\ 200年 x_{2} - - - - - - 200年 {x_{1}}^{2} \\ 2 x_{3} - - - - - - 400年 x_{3} (x_{4} - - - - - - {x_{3}}^{2}) - - - - - - 2 \\ 200年 x_{4} - - - - - - 200年 {x_{3}}^{2} \\ 2 x_{5} - - - - - - 400年 x_{5} (x_{6} - - - - - - {x_{5}}^{2}) - - - - - - 2 \\ 200年 x_{6} - - - - - - 200年 {x_{5}}^{2} \\ 2 x_{7} - - - - - - 400年 x_{7} (x_{8} - - - - - - {x_{7}}^{2}) - - - - - - 2 \\ 200年 x_{8} - - - - - - 200年 {x_{7}}^{2} \\ 2 x_{9} - - - - - - 400年 x_{9} (x_{10} - - - - - - {x_{9}}^{2}) - - - - - - 2 \\ 200年 x_{10} - - - - - - 200年 {x_{9}}^{2} \end{array})$

中所示的中间结果解决非线性系统不包括雅可比矩阵(优化工具箱),使用相同的形式F。

F (1:2: n) =简化(F (1:2: n) + 2 * x (1:2: n)。* F (2:2: n));F (1:2: n) = - F (1:2: n) / 2;(2:2:n) = F (2:2: n) / 20;

现在的系统方程是相同的:

F (1:10)

ans =
                  
                    $(\begin{array}{c} 1 - - - - - - x_{1} \\ 10 x_{2} - - - - - - 10 {x_{1}}^{2} \\ 1 - - - - - - x_{3} \\ 10 x_{4} - - - - - - 10 {x_{3}}^{2} \\ 1 - - - - - - x_{5} \\ 10 x_{6} - - - - - - 10 {x_{5}}^{2} \\ 1 - - - - - - x_{7} \\ 10 x_{8} - - - - - - 10 {x_{7}}^{2} \\ 1 - - - - - - x_{9} \\ 10 x_{10} - - - - - - 10 {x_{9}}^{2} \end{array})$

计算雅可比矩阵表示的方程组

使用雅可比矩阵计算雅可比矩阵F。这个函数计算雅可比矩阵象征意义,从而避免错误与数值近似相关的衍生品。

摩根富林明=雅可比矩阵(F (x);

显示前10雅可比矩阵的行和列。

摩根富林明(1:10,1:10)

ans =
                  
                    $(\begin{array}{c} - - - - - - 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ - - - - - - 20. x_{1} & 10 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & - - - - - - 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & - - - - - - 20. x_{3} & 10 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & - - - - - - 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & - - - - - - 20. x_{5} & 10 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & - - - - - - 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & - - - - - - 20. x_{7} & 10 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - - - - - - 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - - - - - - 20. x_{9} & 10 \end{array})$

大部分的雅可比矩阵的元素摩根富林明是0。不过,当你将这个矩阵转换成一个matlabFunction,结果是一个密集的数字矩阵。通常,稀疏矩阵的操作的效率比相同的密度矩阵的操作。

因此,创建高效的代码,把唯一的非零元素摩根富林明在调用前一个符号向量matlabFunction。是和js表示的稀疏模式摩根富林明。

(是,js) =找到(摩根富林明);摩根富林明=摩根富林明(摩根富林明~ = 0);

方程组和雅可比矩阵转换为一个MATLAB函数

系统的方程F代表。函数,是一个象征性的矩阵,由符号表达式。能够解决的fsolve功能,该系统转换为一个matlabFunction。在这一步中,将很方便F及其雅可比矩阵,摩根富林明一个基于文件的MATLAB函数,FJFfun。原则上,这允许F和摩根富林明重用变量。或者,您可以将它们转换为单独的MATLAB函数。

matlabFunction ([F;摩根富林明),“var”,x,“文件”,“FJFfun”);

FJFfun接受变量列表,但是fsolve预计一个向量的变量。这个函数genericObj将向量转换为一个列表。匿名函数bananaObj定义修复参数值genericObj。注意,使用稀疏模式的功能genericObj。进一步指出,这种方法的使用FJFfun和genericObj在一起不是绑定到特定的表达式为代表F可以作为任何F。

bananaobj = @ (y) genericObj (y, @FJFfun, js)

bananaobj =function_handle与价值:@ (y) genericObj (y, @FJFfun, js)

使用`fsolve`数值求解非线性方程组

使用fsolve对方程组转化为MATLAB函数。使用的起始点x (i) = -1.9奇数项,x (i) = 2偶数项。集“显示”来“通路”看到解决的进展。集的雅可比矩阵来“上”利用雅可比矩阵中定义bananaobj。

x0 (1: n, 1) = -1.9;x0 (2:2: n, 1) = 2;选择= optimoptions (@fsolve,“显示”,“通路”,的雅可比矩阵,“上”);[sol1, Fsol exitflag、输出江淮]= fsolve (x0, bananaobj选项);

一阶规范信赖域迭代Func-count | | f (x) | | ^ 2步最优半径0 329 8563.84 615 3093.71 1 1 2 1 1 2 3 3 4 2.5 34.8 2.5 212.48 6.25 34.1 225.104 6.25 - 4 5 212.48 6.25 34.1 6.25 5 6 6 7 1.5625 34.1 1.56 116.793 0.390625 5.79 212.48 99.4696 0.391 7 8 109.947 0.390625 0.753 0.391 8 9 10 83.6416 2.44141 7.13 2.44 0.976563 1.2 0.977 9 11 11 12 77.7663 2.44141 9.94 2.44 2.44141 9.94 2.44 77.7663 12 13 43.013 0.610352 1.38 0.61 13 14 36.4334 0.610352 1.58 0.61 14 15 15 16 18.0108 1.52588 4.91 1.53 34.1448 1.52588 6.71 1.53 16 17 8.48336 1.52588 3.74 1.53 17 18 18 19 1.52588 3.58 1.53 1.46166 1.52588 3.32 3.74566 1.53 19 20 20 21 0 0 0.0547695 1.53 0.29997 1.24265 1.94 1.53方程解决。fsolve完成因为函数值接近于零的向量的值函数的宽容,问题出现普通的梯度。

使用vpasolve数值求解非线性方程组

非线性方程组,F,或者将其解决vpasolve函数,它是一个数值解算器中可用符号数学工具箱。一个关键的区别vpasolve和fsolve是,fsolve使用双精度运算来解决方程,而vpasolve使用可变精度运算,因此,允许您控制计算的准确性。

sol2 = vpasolve (F);

如果你分配的解决方案的一个变量的方程组,sol2,然后vpasolve返回一个结构数组的形式的万博尤文图斯解决方案。您可以访问每个解决方案(结构数组的每个字段):

sol2.x1

ans =
                   $1.0$

你也可以把sol2一个象征性的向量,然后访问一系列的解决方案。万博尤文图斯例如,显示5日25日解决方案:万博尤文图斯

为k = (1:64) sol2_array (k) = sol2。(char (x (k)));结束sol2_array (25)

ans =
                   $(\begin{array}{c} 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 \end{array})$

辅助函数

函数F [J] = genericObj (x, FJFfun, js)% genericObj采取作为输入,向量x,和雅可比矩阵评价函数%函数处理FJFfun稀疏模式,js并返回作为输出,%的数值函数和雅可比矩阵:F和Jfunction% fj(1:长度(x))是数值函数的向量% fj(长度(x) + 1:结束)是数值向量对雅可比矩阵的非零元素伊势亚= num2cell (x);fj = FJFfun(伊势亚{:});F = fj(1:长度(x));J =稀疏(js, fj(长度(x) + 1:结束));结束