How to solve a set of 4 first-order non-linear coupled ODEs?

1 回表示 (過去 30 日間)
Shashank
Shashank 2013 年 7 月 3 日
I am trying to solve a 4th order differential equation using shooting method by disintegrating the ODE into four coupled first order ODEs. I only have initial conditions. So I am trying to perturb y(1) and y(4) and use them as a parameter to achieve a target value for y(3). But all I get is NaNs on my output matrices for y(1), y(2), y(3) and y(4). Any inputs will be appreciated.
Here is my code:
%%Constants
A = 1e-20; % Hamaker Constant
R_FC = 24.597; % Gas Constant of FC-72 (J/kg-K)
c = 1; % Accomodation Coefficient
sigma = 0.01; % Surface Tension Coefficient (N/m)
rho_l = 1593.84; % Liquid Density (kg/m3)
k_l = 0.0544; % Liquid Conductivity (W/m-K)
mu_l = 0.0004683; % Liquid Viscosity (kg/m2-s)
rho_v = 11.61; % Vapor Density (kg/m3)
h_lv = 88000; % Enthalpy of Vaporization (J/kg)
Twall = 334; % Wall Temperature (K), 10 degrees superheat
Tsat = 329; % Saturation Temperature of FC-72 at 1 atm (K)
R_int = ((2 - c)/(2*c))*(Tsat^(1.5))*((2*pi*R_FC)^(0.5))/(rho_v*(h_lv^2));
delta_ad = (A/(rho_l*h_lv*((Twall/Tsat) - 1)))^(1/3);
e1 = delta_ad/1000;
deltaP_ad = A/((delta_ad)^3);
deltaP_ad_corr = A/((delta_ad+e1)^3);
eQ = 1;
%%Function
tpcl = @(xi,y) [y(2); (1/sigma)*((1+(y(2)*y(2)))^(1.5))*(y(3)-(A/(y(1)^3))); (3*mu_l/(rho_l*h_lv))*(-y(4)/(y(1)^3)); (Twall-Tsat-(Tsat*y(3)/(rho_l*h_lv)))/((y(1)/k_l)+R_int)];
xi_end = 5e-7;
N = 1000;
[xi,y1] = rk4(tpcl,[0 xi_end],[delta_ad 0 deltaP_ad 0],N);
[xi,y2] = rk4(tpcl,[0 xi_end],[delta_ad+e1 0 deltaP_ad_corr eQ],N);
delta_a = delta_ad;
delta_b = delta_ad + e1;
Q_a = 0;
Q_b = eQ;
slope_a = y1(2,end);
slope_b = y2(2,end);
%%Shooting Method
tol = 1e-5;
iter = 10;
slope_target = tan(pi/18);
for i = 1:iter
Q_new = Q_a + (Q_b-Q_a)/(slope_b-slope_a)*(slope_target-slope_a);
delta_new = delta_a + (delta_b-delta_a)/(slope_b-slope_a)*(slope_target-slope_a);
deltaP_new = A/((delta_new)^3);
[xi,y] = rk4(tpcl,[0 xi_end],[delta_new 0 deltaP_new 0],N);
fprintf('iter:%2d, delta(0)=%17.15f, slope(xi_end)=%17.15f\n',i,delta_new,y(2,end));
if (abs(y(2,end)-slope_target) <= tol)
break;
end
Q_a = Q_b;
Q_b = Q_new;
delta_a = delta_b;
delta_b = delta_new;
slope_a = slope_b;
slope_b = y(2,end);
end

採用された回答

Jan
Jan 2013 年 7 月 4 日
編集済み: Jan 2013 年 7 月 4 日
Anonymous functions are hard to debug. Better define the ODE as a function and use http://www.mathworks.de/matlabcentral/answers/1971 to define the constants on demand.
Then you can use the debugger to step through the program line by line until you find the source of the NaN's by your own. This is more efficient to let all forum users debug the code for you, especially when the function rk4() is not known.
  5 件のコメント
Shashank
Shashank 2013 年 7 月 6 日
The NaNs are coming after a few iterations in y1 and y2.
Jan
Jan 2013 年 7 月 7 日
This could have different reasons: A bug in the code, a division by zero or another mathematical effect. It is your turn to find the exact cause for the creation of the NaNs. E.g. this might be useful:
dbstop if naninf

サインインしてコメントする。

その他の回答 (1 件)

Gentian Zavalani
Gentian Zavalani 2013 年 7 月 7 日
When you construct an anonymous function, the part directly after the @ must be pure variable names and not expressions or indexed variables. Your code is tpcl = @(xi,y) [y(2); (1/sigma)*((1+(y(2)*y(2)))^(1.5))*(y(3)-(A/(y(1)^3))); (3*mu_l/(rho_l*h_lv))*(-y(4)/(y(1)^3)); (Twall-Tsat-(Tsat*y(3)/(rho_l*h_lv)))/((y(1)/k_l)+R_int)];

カテゴリ

Help Center および File ExchangeOrdinary Differential Equations についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by