|
楼主 |
发表于 2013-6-3 00:58:51
|
显示全部楼层
%得到一阶导数的值,公式如下
%2 * sin(b * x) * (c - y + a *sin(b * x))
%2 * a * x * cos(b * x) * (c - y + a * sin(b * x))
%2 * c - 2 * y + 2 * a * sin(b * x)
function Res = GetYijieDaoshu(Xy, Para)
Res = zeros(3, 1);
for i = 1 : size(Xy, 1)
Res(1, 1) = Res(1, 1) + 2 * sin(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * sin(Para(2) * Xy(i, 1)));
Res(2, 1) = Res(2, 1) + 2 * Para(1) * Xy(i, 1) * cos(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * sin(Para(2) * Xy(i, 1)));
Res(3, 1) = Res(3, 1) + 2 * Para(3) - 2 * Xy(i, 2) + 2 * Para(1) * sin(Para(2) * Xy(i, 1));
end
end
%得到二阶导数的值,公式如下
%2 * sin(b * x)^2
%2 * x * cos(b * x) * (c - y + a * sin(b * x)) + 2 * a * x * cos(b * x) * sin(b * x)
%2 * sin(b * x)
%2 * x * cos(b * x) * (c - y + a * sin(b * x)) + 2 * a * x * cos(b * x) * sin(b * x)
%2 * a^2 * x^2 * cos(b * x)^2 - 2 * a * x^2 * sin(b * x) * (c - y + a * sin(b * x))
%2 * a * x * cos(b * x)
%2 * sin(b * x)
%2 * a * x * cos(b * x)
%2
function Res = GetErjieDaoshu(Xy, Para)
Res = zeros(3, 3);
for i = 1 : size(Xy, 1)
Res(1, 1) = Res(1, 1) + 2 * sin(Para(2) * Xy(i, 1))^2;
Res(1, 2) = Res(1, 2) + 2 * Xy(i, 1) * cos(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * sin(Para(2) * Xy(i, 1))) + 2 * Para(1) * Xy(i, 1) * cos(Para(2) * Xy(i, 1)) * sin(Para(2) * Xy(i, 1));
Res(1, 3) = Res(1, 3) + 2 * sin(Para(2) * Xy(i, 1));
Res(2, 1) = Res(2, 1) + 2 * Xy(i, 1) * cos(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * sin(Para(2) * Xy(i, 1))) + 2 * Para(1) * Xy(i, 1) * cos(Para(2) * Xy(i, 1)) * sin(Para(2) * Xy(i, 1));
Res(2, 2) = Res(2, 2) + 2 * Para(1)^2 * Xy(i, 1)^2 * cos(Para(2) * Xy(i, 1))^2 - 2 * Para(1) * Xy(i, 1)^2 * sin(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * sin(Para(2) * Xy(i, 1)));
Res(2, 3) = Res(2, 3) + 2 * Para(1) * Xy(i, 1) * cos(Para(2) * Xy(i, 1));
Res(3, 1) = Res(3, 1) + 2 * sin(Para(2) * Xy(i, 1));
Res(3, 2) = Res(3, 2) + 2 * Para(1) * Xy(i, 1) * cos(Para(2) * Xy(i, 1));
Res(3, 3) = Res(3, 3) + 2;
end
end
%计算系数的值
function Res = CalPara(Xy, ParaInit, N)
Res = ParaInit;
for i = 1 : N
ParaN = GetYijieDaoshu(Xy, Res);
ParaNN = GetErjieDaoshu(Xy, Res);
Res = Res - 1 * (ParaNN \ ParaN);
end
end
clear;
clc;
a = 1.7;
b = 3.3;
c = -4.6;
X = 200 : 0.01 : 203;
Y = a * sin(b * X) + c + (rand(size(X)) - 0.5) / 1.0;
ParaABC = CalPara([X', Y'], [20.91, 3.3, -34.6]', 40);
Z = ParaABC(1) * sin(ParaABC(2) * X) + ParaABC(3);
plot(X, Y, X, Z);
|
|