ODE15s with non-constant Jacobian

2 ビュー (過去 30 日間)
Berry
Berry 2022 年 8 月 2 日
コメント済み: Berry 2022 年 8 月 3 日
Hi everyone,
i want to solve an ODE by using the Jacobian, that is not constant:
options = odeset( 'Jacobian' , @(t,x)J(u_int, deltaZ, Nz, cFeed), 'RelTol',1e-8, 'AbsTol',1e-8);
[t,V] = ode15s(@(t,v)odeTest(t,v), tspan, v0, options);
function dfdx = J(u_int, deltaZ, Nz, cFeed)
x = sym('x' , [1 Nz+4]);
f1 = (-u_int .* (x(1)^2 - cFeed) ./ deltaZ );
f2 = (-u_int .* (x(2:Nz+4) - x(1:Nz+3)) ./ deltaZ );
f = [f1, f2];
J = jacobian(f,x);
dfdx = (J);
end
%% ODE
function dfdt = odeTest(~, vi)
dfdt = zeros(length(vi),1);
%% Variables
u_int = 0.10;
H = 0.20;
Nz = 100;
deltaZ = H/(Nz+4);
cFeed = 10;
x = vi(1:Nz+4);
%% DGL
dfdt(1) = (-u_int .* (x(1)^2 - cFeed) ./ deltaZ );
dfdt(2: Nz+4) = (-u_int .* (x(2:Nz+4) - x(1:Nz+3)) ./ deltaZ );
end
The problem I am facing is following error:
Conversion to logical from sym is not possible.
Error in ode15s (line 405)
if absh * rh > 1
Error in test_jacobi_upwind (line 39)
[t,V] = ode15s(@(t,v)odeTest(t,v), tspan, v0, options);
I can make it work, if the Jacobian is a constant and does not depend on x(i) - but as soon as J is depending on x(i), I am getting errors.
Any help is appreciated.
Best,
M

採用された回答

Torsten
Torsten 2022 年 8 月 2 日
JacFun = J(u_int, deltaZ, Nz, cFeed);
options = odeset( 'Jacobian' , @(t,x)JacFun(t, x, u_int, deltaZ, Nz, cFeed), 'RelTol',1e-8, 'AbsTol',1e-8);
[t,V] = ode15s(@(t,x)odeTest(t,x,u_int, deltaZ, Nz, cFeed), tspan, v0, options);
function dfdx = J(u_int, deltaZ, Nz, cFeed)
x = sym('x' , [1 Nz+4]);
f1 = (-u_int .* (x(1)^2 - cFeed) ./ deltaZ );
f2 = (-u_int .* (x(2:Nz+4) - x(1:Nz+3)) ./ deltaZ );
f = [f1, f2];
J = jacobian(f,x);
dfdx = matlabFunction(J,'Vars',[t, x, u_int, deltaZ, Nz, cFeed]);
end
function dfdt = odeTest(t, vi, u_int, deltaZ, Nz, cFeed)
dfdt = zeros(length(vi),1);
%% Variables
u_int = 0.10;
H = 0.20;
Nz = 100;
deltaZ = H/(Nz+4);
cFeed = 10;
x = vi(1:Nz+4);
%% DGL
dfdt(1) = (-u_int .* (x(1)^2 - cFeed) ./ deltaZ );
dfdt(2: Nz+4) = (-u_int .* (x(2:Nz+4) - x(1:Nz+3)) ./ deltaZ );
end
  5 件のコメント
Torsten
Torsten 2022 年 8 月 3 日
I didn't use this options command in my answer:
options = odeset( 'Jacobian' , @(t,x)JacFunc(u_int, deltaZ, Nz, epsilon,cFeed), 'RelTol',1e-8, 'AbsTol',1e-8);
Berry
Berry 2022 年 8 月 3 日
My mistake! Thank you very much! Everything works fine now.

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

その他の回答 (0 件)

カテゴリ

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

製品


リリース

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by