ODE求解未知限制

Question

我試圖解決從Y_1到Y_2常微分方程的系統，說

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

與初始條件

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

我的問題是，我做的不知道y_1和y_2，爲了解決這個問題，我當然需要兩個額外的等式，即

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

所以遠遠我已經嘗試使用fsolve（我猜新名稱是fzero）來實現它以找到零，然後從函數F_1和F_2我調用ode45從y_1到y_2求解以便計算函數。但它不起作用，我似乎無法找到任何錯誤。因此，我尋找新的方法/想法，我發現了方法bvp4c，但我不知道我是否可以在我的情況下使用它。有沒有人有與bvp4c的經驗，並知道是否和如何使用它，或者你有其他的想法？ 任何幫助表示讚賞。

代碼以供參考：

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]); 

Answer 1

假設定義域是足夠大，你有很好的初始值，牛頓實現與雅可比的建設分差看起來是這樣的：

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