MATLAB Answers

Problems with integral in fsolve

3 views (last 30 days)
Hey
I have three functions with three unknown variables, but I can't find a solution for them. the features look like this:
functions.PNG
I've tried to do the following code:
m = 0.31;
r = 0.05;
delta = 0.02;
y_bar = 2;
w_1 = 5;
w_2_bar = 5;
sigma = 1;
b_bar = 3;
z_0=[0.5,3,4];
sol=fsolve(@(z)([
w_1*(1-z(1))-z(2)-((1+delta)/((1-m)*(1+r)))*int((b_bar+(1-z(1))*(1-m)*x+(1-m)*(1+r)*z(2))*(1/(2*pi*sigma^2))*exp(-((x-w_2_bar)^2)/(2*sigma^2)),x,0,z(3))
+((1+delta)/(1+r))*int(((1-z(1))*k+(1+r)*z(2))*(1/(2*pi*sigma^2))*exp(-((k-w_2_bar)^2)/(2*sigma^2)),k,z(3),inf);
z(1)*w_1+int(z(1)*h*(1/(2*pi*sigma^2))*exp(-((h-w_2_bar)^2)/(2*sigma^2)),h,0,z(3))
-int((b_bar-m*((1-z(1))*y+(1+r)*z(2)))*(1/(2*pi*sigma^2))*exp(-((y-w_2_bar)^2)/(2*sigma^2)),y,0,z(3));
z(3)-(1/(1-z(1))*y_bar+((1+r)/(1-z(1)))*z(2))]),z_0);
Hope there is one who can correct my code so it works, thank you very much if it is you. :)

  0 Comments

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 9 Apr 2019
Code attached.
It is not fast code. You could improve the performance by changing the symbolic integrations into numeric integrations.
Or, since you are using symbolic integration, you could process your function once with symbolic z variables, and simplify. vpasolve() does a good job on what is left, or you can matlabFunction() and fsolve()

  1 Comment

Walter Roberson
Walter Roberson on 9 Apr 2019
Having the int() inside the fsolve() is slow, and the int() have closed form representations if you use symbolic z. The performance gains from executing once symbolically and using the result with vpasolve() or fsolve() is substantial.

Sign in to comment.

More Answers (1)

Torsten
Torsten on 9 Apr 2019
Edited: Torsten on 9 Apr 2019
function main
z_0 = [0.5,3,4];
options = optimset('TolFun',1e-8,'TolX',1e-8);
sol = fminsearch(@fun,z_0,options)
end
function res = fun(z)
m = 0.31;
r = 0.05;
delta = 0.02;
y_bar = 2;
w_1 = 5;
w_2_bar = 5;
sigma = 1;
b_bar = 3;
re(1) = w_1*(1-z(1))-z(2)-((1+delta)/((1-m)*(1+r)))*integral(@(x)(b_bar+(1-z(1))*(1-m)*x+(1-m)*(1+r)*z(2))*(1/(2*pi*sigma^2))*exp(-((x-w_2_bar).^2)/(2*sigma^2)),0,z(3),'ArrayValued',1)...
+((1+delta)/(1+r))*integral(@(k)((1-z(1))*k+(1+r)*z(2))*(1/(2*pi*sigma^2)).*exp(-((k-w_2_bar).^2)/(2*sigma^2)),z(3),inf,'ArrayValued',1);
re(2) = z(1)*w_1+integral(@(h)z(1)*h*(1/(2*pi*sigma^2)).*exp(-((h-w_2_bar).^2)/(2*sigma^2)),0,z(3),'ArrayValued',1)...
-integral(@(y)(b_bar-m*((1-z(1))*y+(1+r)*z(2)))*(1/(2*pi*sigma^2)).*exp(-((y-w_2_bar).^2)/(2*sigma^2)),0,z(3),'ArrayValued',1);
re(3) = z(3)-(1/(1-z(1))*y_bar+((1+r)/(1-z(1)))*z(2));
res = norm(re)
end

  0 Comments

Sign in to comment.

Sign in to answer this question.

Products


Release

R2018b