Solving ODE with unknown limits

Question

I'm trying to solve a system of ordinary differential equations from y_1 to y_2, say

G' = D
D' = f(y,G,D)

with the initial conditions

G(y_1) = 0
D(y_1) = 0

My problem is that I do not know y_1 and y_2, to counter for this I of course need two additional equations which is

F_1(y_1,y_2,G(y_2)) = 0
F_2(y_1,y_2,G(y_2)) = 0

So far I have tried to implement it using fsolve (i guess the new name is fzero) to find the zeros, then from the functions F_1 and F_2 i call the ode45 to solve from y_1 to y_2 in order to calculate the functions. However it does not work and I can not seem to find any mistakes. Therefore i looked for new methods/ideas and I found the method bvp4c, but I'm not sure whether i can use it in my case. Does anyone have experience with bvp4c and know whether and how I should use it, or do you have other ideas? Any help is appreciated.

Code for reference:

function [Fval,sol,t,G] = EntrepreneurialFinanceNondiversifiableRisk
global sigma rf tau b epsilon gamma theta1 theta2 alpha ya K I taug mu eta...    
omega
sigma=0.2236; rf=0.03; tau=0.1129; b=0.85; epsilon=0.2; gamma=2;
theta1 = -0.704; alpha = 0.6; theta2=1.704; ya=0.1438; K=27; I = 10;
taug = 0.1; mu = 0.04; eta = 0.4; omega = 0.1;
tau = 0;

option = optimset('Display','iter');
sol = fsolve(@f,[0.1483,2.8],option);
Fval = f(sol);
[t,G] = solvediff(sol(1),sol(2));

function F = f(x)
yd = x(1);
yu = x(2);
global rf tau b theta1 alpha theta2 K I taug
Vstar =@ (y) (1-tau+tau*(1-theta1-(1-alpha)*(1-tau)* theta1/tau)^...
(1/theta1))*y/rf;
VstarPrime = (1-tau+tau*(1-theta1-(1-alpha)*(1-tau)...
*theta1/tau)^(1/theta1))*1/rf;
qbar =@(yd,yu) (yd^theta1-yd^theta2)/(yu^theta2*yd^theta1-yu^theta1*...
yd^theta2);
qunderbar =@ (yd,yu) (yu^theta2-yu^theta1)/(yu^theta2*yd^theta1-...
yu^theta1*yd^theta2);
A =@ (y) (1-tau)*(y/rf);
V = @(y,yd) A(y) + (tau * b)/rf * (1 - (y/yd)^(theta2)) - (1-alpha)*A(yd)*...
(y/yd)^(theta2); 
F0 =@ (yd,yu) b/rf-(b/rf-alpha*A(yd))*qunderbar(yd,yu)/(1-qbar(yd,yu));
[t,G] = solvediff(yd,yu);
F = zeros(2,1);

F(1) = Vstar(yu)-F0(yd,yu)-K-taug*(Vstar(yu)-K-I)-G(end,2);
F(2) = (1-taug)*VstarPrime-G(end,2);

function [t,G] = solvediff(yd,yu)
[t,G] = ode45(@diff,[yd,yu],[0,0]);

Could we take a look at the code you are using? What is the function you are trying to solve? — brodoll
I have changed the description to now include the code. But I am curious how to solve this in general. — Nicky Mattsson
I would suggest an iterative procedure trying to find a convergence between ODE and algebraic equations involving ODE solutions and initial conditions. IN other word, set one solve another and vice versa. — freude
But I can't split the two algebraic equations, so if I set one and solve for the other one I'm not sure that there exist any solution for that variable. Or have I misunderstood anything? — Nicky Mattsson

Lutz Lehmann Lutz Lehmann · Accepted Answer · 2015-05-04T15:38:40

Assuming the domain of definition is big enough and you have good initial values, a Newton implementation with divided differences for the construction of the Jacobian looks like this:

h=1e-5

for k=1:10 do
    Fval = F(yd, yu)
    dFdyd = (F(yd+h, yu)-F(yd-h, yu))/(2*h)
    dFdyu = (F(yd, yu+h)-F(yd, yu-h))/(2*h)

    dF = [ dFdyd dFdyu ]

    ynext = [ yd; yu ] - dF^(-1)*Fval;

    yd = ynext(1); yu = ynext(2);

end

Solving ODE with unknown limits

1 Answers