【matlab】多元多维非线性函数的迭代求精确解（自用）

本文链接：https://blog.youkuaiyun.com/2302_79075821/article/details/144780142

待解决问题：

用六个有源RC滤波器的一阶节，实现五阶平坦的90°移相器。

电路图：

目标：

1. 一端输入正弦波，另一端输出两个正弦波，且相移为90°。

2. 实现在足够大带宽上的五阶平坦。

数学分析：

易知相移满足

$\varphi (\omega) =2\left(arctan(\frac{\omega }{\omega_1})+arctan(\frac{\omega }{\omega_2})+arctan(\frac{\omega }{\omega_3})-arctan(\frac{\omega }{\omega_4})-arctan(\frac{\omega }{\omega_5})-arctan(\frac{\omega }{\omega_6})\right) =\frac{\pi}{2}$

其中

$\omega_k=\frac{1}{R_kC_k},k=1,2,...,6$

其导数又满足

$\varphi^{(i)}=0,i=1,2,...,5$

故给定一个 $\omega_0$ ，要找到 $\omega_1,...,\omega_6$ 满足上述公式，可视为一个六元六维的向量函数求零点。

MATLAB迭代解法如下：

（注意：MATLAB自身的fsolve函数本身已经内蕴迭代法，可提供一定精度下的解。

但本题中各阶导数本身已相当小，且数量级相差太大，用一次fsolve无法得到满足要求的精确解。

故将各阶导数按不同比重扩倍，同时使函数值随迭代次数也扩倍，不让fsolve轻易判定“已解出”。

btw：不知道本题是否属于那类“条件数>1”的病态问题？（马上期末考试来不及验证力！（悲）

clc;
clear;

N = 1000000;

syms w w1 w2 w3 w4 w5 w6;
f0 = (2*(atan(w/w1)+atan(w/w2)+atan(w/w3)...
    -atan(w/w4)-atan(w/w5)-atan(w/w6))-pi/2)*N;
f1 = diff(f0,w)*1000;
f2 = diff(f0,w,2)*10000;
f3 = diff(f0,w,3)*100000;
f4 = diff(f0,w,4)*1000000;
f5 = diff(f0,w,5)*10000000;
f = [f0,f1,f2,f3,f4,f5];
f = matlabFunction(f);
f0 = matlabFunction(f0/N*180/pi+90);

function F = myfun(ww,w0,i,f)
    F = 3.9^i.*f(w0,ww(1),ww(2),ww(3),ww(4),ww(5),ww(6));
end

w0 = 800*2*pi;
ww = [10000,3000,1000,30000,5000,1000];
for i = 1:500
    [ww] = fsolve(@(ww)myfun(ww,w0,i,f),ww);
    myfun(ww,w0,i,f);
end

ww
myfun(ww,w0,0,f)
r = 1e8./ww

运行结果：

多次验证可知方程的解确实收敛（本人调参要把心态调崩了）（悲）

下面验证相移曲线确实足够平坦：

line = pi/2*0.99;
pp = plot([-2,2],[line,line],'-- k');
hold on;

fw = -2:0.0001:2;
ff0 = f0((10.^fw)*w0,ww(1),ww(2),ww(3),ww(4),ww(5),ww(6))*pi/180;
k = find(abs(ff0-line)<=1e-5);
p0 = plot(fw,ff0,'r');
plot([fw(k(1)),fw(k(1))],[0,line],': k');
plot([fw(k(2)),fw(k(2))],[0,line],': k');
bw1 = 10^(fw(k(2)))-10^(fw(k(1)));

ww = [8694,2849,468,27436,4738,1527];
ff1 = f0((10.^fw)*w0,ww(1),ww(2),ww(3),ww(4),ww(5),ww(6))*pi/180;
k = find(abs(ff1-line)<=1e-5);
p1 = plot(fw,ff1,'g');
plot([fw(k(1)),fw(k(1))],[0,line],'-. k');
plot([fw(k(2)),fw(k(2))],[0,line],'-. k');
bw2 = 10^(fw(k(2)))-10^(fw(k(1)));

xlabel('$lg\frac{w}{w_0}~~~(w_0=800*2pi)$','Interpreter','latex');
ylabel('$phy/rad$','Interpreter','latex');
legend([pp,p0,p1],{['$phy=$',num2str(line)],...
    ['$\frac{BW_1}{w_0}=$',num2str(bw1)],...
    ['$\frac{BW_2}{w_0}=$',num2str(bw2)]},'Interpreter','latex');

运行结果如下：