Solve non linear equation with vector

3 ビュー (過去 30 日間)
Guilherme Lopes de Campos
Guilherme Lopes de Campos 2023 年 8 月 11 日
コメント済み: Sam Chak 2023 年 8 月 11 日
Dear,
I would like to solve a non linear equation, on the form:
yv = 1:1:19;
for k1 = 1:length(yv)
y = yv(k1);
f = @(x) exp(-0.5*((x-45)/dados.mass(k1))^2) == dados.par(k1) ;
x(k1) = fzero(f, 0.5);
end
Send : x the variable
dados.mass(k1) vector with value, contain 19 lines
dados.par (k1) equal aspect to top.
I would like to solve each equation related a value of dados.mass(1) until (19) concomitantly at dados.par(1) until (19), obtaining the vector x with respectivly values.
Can help me, please?
Yours faithfully

回答 (3 件)

Torsten
Torsten 2023 年 8 月 11 日
移動済み: Torsten 2023 年 8 月 11 日
x(k1) = 45+dados.mass(k1)*sqrt(-2*log(dados.par(k1)))
or
x(k1) = 45-dados.mass(k1)*sqrt(-2*log(dados.par(k1)))
  1 件のコメント
Sam Chak
Sam Chak 2023 年 8 月 11 日
@Torsten, 👍 That's a good analytical solution without using fzero()! Did I plot the visualization of the solutions correctly?
m = linspace(0, 300, 31);
p = linspace(0, 1, 31); % assuming that Gaussian function has a range of 0 to 1
[M, P] = meshgrid(m, p);
X1 = 45 - sqrt(2)*M.*sqrt(log(1./P));
X2 = 45 + sqrt(2)*M.*sqrt(log(1./P));
tiledlayout(2,2);
nexttile
surf(M, P, X1), xlabel('mass'), ylabel('par'),
title({'$x^{-}$'}, 'interpreter', 'latex', 'fontsize', 16)
nexttile
surf(M, P, X2), xlabel('mass'), ylabel('par'),
title({'$x^{+}$'}, 'interpreter', 'latex', 'fontsize', 16)
nexttile
contourf(M, P, X1), xlabel('mass'), ylabel('par')
nexttile
contourf(M, P, X2), xlabel('mass'), ylabel('par')

サインインしてコメントする。


Star Strider
Star Strider 2023 年 8 月 11 日
It looks like you want to do a nonlinear regression.
Perhaps this —
dados = array2table(sortrows([10*randn(12,1)+35,rand(12,1)],1), 'VariableNames',{'mass','par'})
dados = 12×2 table
mass par ______ _______ 20.408 0.66176 21.47 0.37043 24.976 0.14994 28.527 0.72292 29.961 0.91367 33.376 0.67683 33.456 0.64286 34.064 0.2926 34.807 0.49665 37.383 0.91716 39.025 0.984 58.526 0.60612
f = @(x,m) exp(-0.5*((x-45)./m).^2) ;
x0 = rand;
mdl = fitnlm(dados, f, x0)
mdl =
Nonlinear regression model: par ~ F(x,m) Estimated Coefficients: Estimate SE tStat pValue ________ ______ _____ _________ b1 15.791 3.7571 4.203 0.0014781 Number of observations: 12, Error degrees of freedom: 11 Root Mean Squared Error: 0.25 R-Squared: 0.0592, Adjusted R-Squared 0.0592 F-statistic vs. zero model: 74.6, p-value = 3.13e-06
[y,yci] = predict(mdl, dados.mass);
figure
hp1 = plot(dados.mass, dados.par, 'b.', 'DisplayName','Data');
hold on
hp2 = plot(dados.mass, y, 'DisplayName','Regression');
hp3 = plot(dados.mass, yci, '--r', 'DisplayName','95% Confidence Intervals');
hold off
grid
xlabel('mass')
ylabel('par')
legend([hp1,hp2,hp3(1)], 'Location','best')
If you do not have the Statistics and Machine Learning Toolbox (for fitnlm), an alternative would be the fminsearch function —
xest = fminsearch(@(x) norm(dados.par - f(x,dados.mass)), x0)
xest = 15.7910
figure
hp1 = plot(dados.mass, dados.par, 'b.', 'DisplayName','Data');
hold on
hp2 = plot(dados.mass, f(xest,dados.mass), 'DisplayName','Regression');
hold off
grid
xlabel('mass')
ylabel('par')
legend([hp1,hp2], 'Location','best')
.

Sam Chak
Sam Chak 2023 年 8 月 11 日
I have fixed the function f in the code. Now it should be working correctly. The dados.mass and dados.par data are not provided. Thus I made up some value to test the code.
mass = linspace(110, 190, 19);
par = linspace(0.1, 0.19, 19);
yv = 1:1:19;
for k1 = 1:length(yv)
% y = yv(k1); % unused
f = @(x) exp(-0.5*((x - 45)/mass(k1)).^2) - par(k1); % <-- fix it here
x(k1) = fzero(f, 0.5);
end
% Solutions
x
x = 1×19
-191.0563 -197.9780 -204.7954 -211.5110 -218.1269 -224.6453 -231.0681 -237.3972 -243.6343 -249.7812 -255.8394 -261.8103 -267.6954 -273.4959 -279.2132 -284.8485 -290.4028 -295.8772 -301.2727
subplot(2,1,1)
plot(mass, x), xlabel('mass'), ylabel('x'), grid on
subplot(2,1,2)
plot(par, x), xlabel('par'), ylabel('x'), grid on

カテゴリ

Help Center および File ExchangeSupport Vector Machine Regression についてさらに検索

製品


リリース

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by