Writing function based on a vector of roots

1 visualización (últimos 30 días)
John
John el 28 de Feb. de 2018
Comentada: Star Strider el 1 de Mzo. de 2018
In the following code, one of lines plotted is a vector of solutions to a nonlinear function (plotted against a parameter). I want to find the intersection point, but I don't know how to do it. Does anyone have any suggestion?
clc;
close all;
clear all;
lambda = 1;
Fprime = @(z) lambda.*exp(-lambda.*z);
F = @(z) 1-lambda.*exp(-lambda.*z);
F2prime = @(z) -lambda^2.*exp(-lambda.*z);
intz = @(z) (1-(lambda.*z+1).*exp(-lambda.*z))./lambda;
funz = @(z) z.*(1-F(z))+intz(z);
ph = 0.8;
pl = 0.4;
t = 0.6;
theta =1.3;
sigma = 0.6;
alpha=0.6;
epsilon = 0.1;
y = @(z,psi) (((ph-pl)/(1-t))*psi*(z.*(1-F(z))+intz(z))).^(1/theta);
u = @(z,psi) ph.*((1+epsilon).*y(z,psi)).^alpha + (1-ph).*y(z,psi).^alpha;
N2 = 200;
psigrid = linspace(.5,5,N2);
zbar2 = zeros(1,N2);
Gfun = @(z,psi) psi.*z;
rhs = @(z,psi) (sigma/.8).*(u(z,psi)-Gfun(z,psi)+(1-ph).*psi.*(1-F(z)).*z);
for k=1:N2
rhstest = @(z) rhs(z,psigrid(1,k));
Gfuntest = @(z) Gfun(z,psigrid(1,k));
zbarfun = @(z) rhstest(z) - Gfuntest(z);
zbar2(1,k) = fsolve(zbarfun,.2);
end
figure
plot(psigrid,zbar2);
hold on;
plot(psigrid,.3.*ones(1,N2));
  3 comentarios
John
John el 1 de Mzo. de 2018
Strangely enough, it's still not showing...
Thank you very much for your help!
Star Strider
Star Strider el 1 de Mzo. de 2018
As always, my pleasure!
I’ve not heard back yet from Rena Berman. I’ll let you know. She might post back here as well, since I included the URL to your Question in my email to her.

Iniciar sesión para comentar.

Respuesta aceptada

Star Strider
Star Strider el 1 de Mzo. de 2018
This will calculate the x-coordinate of the intersection, and plot a green pentagram there:
intx = interp1(zbar2, psigrid, 0.3, 'linear');
figure
plot(psigrid,zbar2)
hold on;
plot(psigrid,.3.*ones(1,N2))
plot(intx, 0.3, 'pg', 'MarkerFaceColor','g', 'MarkerSize',10)
hold off
The rest of your code (before and including the for loop) is unchanged, so I didn’t post it.

Más respuestas (0)

Categorías

Más información sobre MATLAB en Help Center y File Exchange.

Productos

Community Treasure Hunt

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

Start Hunting!

Translated by