Solve non linear equation with vector

5 visualizzazioni (ultimi 30 giorni)
Guilherme Lopes de Campos
Guilherme Lopes de Campos il 11 Ago 2023
Commentato: Sam Chak il 11 Ago 2023
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

Risposte (3)

Torsten
Torsten il 11 Ago 2023
Spostato: Torsten il 11 Ago 2023
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 Commento
Sam Chak
Sam Chak il 11 Ago 2023
@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')

Accedi per commentare.


Star Strider
Star Strider il 11 Ago 2023
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 il 11 Ago 2023
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

Categorie

Scopri di più su Introduction to Installation and Licensing in Help Center e File Exchange

Prodotti


Release

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by