Writing function based on a vector of roots

6 vues (au cours des 30 derniers jours)
John
John le 28 Fév 2018
Commenté : Star Strider le 1 Mar 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 commentaires
John
John le 1 Mar 2018
Strangely enough, it's still not showing...
Thank you very much for your help!
Star Strider
Star Strider le 1 Mar 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.

Connectez-vous pour commenter.

Réponse acceptée

Star Strider
Star Strider le 1 Mar 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.

Plus de réponses (0)

Catégories

En savoir plus sur MATLAB dans Help Center et File Exchange

Produits

Community Treasure Hunt

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

Start Hunting!

Translated by