Azzera filtri
Azzera filtri

BVP Error, help with explaination

27 visualizzazioni (ultimi 30 giorni)
Thanh Hoang
Thanh Hoang il 23 Mag 2024 alle 10:32
Commentato: Ayush Anand il 24 Mag 2024 alle 4:52
Hello all,
I am struggeling a bit with understanding why I get an error message here. The code works fine with bigger numbers for A or with smaller numbers for B but is crashing for example with the given parameters. I already figured that y(1) is becoming NaN at some point which causes the code to crash. Could somebody help me understand why this is happening?
bvp()
Error using bvp4c (line 196)
Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in solution>bvp (line 8)
sol = bvp4c(@ode, @bcfun, sol_init);
function bvp
A = 50e-3;
B = 1;
xmesh = linspace(0, 1, 100);
sol_init = bvpinit(xmesh, @init);
sol = bvp4c(@ode, @bcfun, sol_init);
function dydx = ode(x, y)
dydx(1) = y(2);
dydx(2) = exp(y(1)*log(10)/A);
%disp(y(1))
end
% Define the boundary conditions
function res = bcfun(ya, yb)
res(1) = ya(1) - B;
res(2) = yb(1);
end
% Define the initial guess for the solution
function guess = init(x)
guess(1) = B;
guess(2) = 0;
end
end

Risposte (2)

Ayush Anand
Ayush Anand il 23 Mag 2024 alle 10:53
Modificato: Ayush Anand il 23 Mag 2024 alle 11:08
The magnitude of the exponential function exp(y(1)*log(10)/A) blows up very rapidly as A is . This is most probably causing the overflow, leading to an Inf (infinity) value in MATLAB. Further calculations involving Inf can lead to NaN (Not a Number), as these are operations that MATLAB cannot numerically evaluate to a real or infinite number.
Also, the solver bvp4c iteratively adjusts y(1) and y(2) to satisfy both the differential equation and the boundary conditions. If during its iterations, y(1) becomes large enough, exp(y(1)*log(10)/A) will overflow. Once y(1) or any subsequent value becomes NaN, the solver cannot proceed with numeric computations and thus the error you see.
  2 Commenti
Thanh Hoang
Thanh Hoang il 23 Mag 2024 alle 11:05
Thanks a lot Ayush!
Is there a way to overcome this? Or is this system not solvable numerically?
Ayush Anand
Ayush Anand il 24 Mag 2024 alle 4:52
I think as mentioned by Torsten in the other answer, the solver reaches the limits of computation/precision possible on the machine and hence can't proceed with any further feasible steps. It seems to be associated with the system of equations itself and you may want to change the formulation of the system.

Accedi per commentare.


Torsten
Torsten il 23 Mag 2024 alle 11:37
Modificato: Torsten il 23 Mag 2024 alle 13:15
The solver is no longer able to capture the steep gradient of your solution at x=0.
Look at sol.x in the last step of the continuation method from below and how fine the mesh has to be chosen around x = 0.
bvp()
ans = 1x64
1.0e+00 * 0 0.000000000125289 0.000000000250578 0.000000000375868 0.000000000751735 0.000000001127603 0.000000002255205 0.000000006201814 0.000000010148424 0.000000025371059 0.000000040593694 0.000000060890541 0.000000213116894 0.000000365343247 0.000000548014871 0.000000730686495 0.000001803882283 0.000002877078072 0.000005822658004 0.000008768237935 0.000017536475870 0.000026304713805 0.000048225308642 0.000092546011329 0.000136866714015 0.000199340409301 0.000329630944865 0.000459921480429 0.000717008956755 0.000974096433081
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
function bvp
A_array = [5e-1,3e-1,2e-1,1.5e-1,1e-1,0.9e-1,0.8e-1,0.7e-1,0.6e-1,0.55e-1,0.5e-1];
B = 1;
xmesh = linspace(0, 1, 100);
sol_init = bvpinit(xmesh, [0 B]);
hold on
for i = 1:numel(A_array)
A = A_array(i);
sol = bvp4c(@ode, @bcfun, sol_init);
plot(sol.x,sol.y(1,:))
init = @(x)interp1(sol.x.',sol.y.',x);
sol_init = bvpinit(sol.x,init);
end
format long
sol.x
hold off
grid on
function dydx = ode(x, y)
dydx(1) = y(2);
dydx(2) = exp(y(1)*log(10)/A);
%disp(y(1))
end
% Define the boundary conditions
function res = bcfun(ya, yb)
res(1) = ya(1) - B;
res(2) = yb(1);
end
end

Categorie

Scopri di più su Numerical Integration and Differential Equations in Help Center e File Exchange

Prodotti


Release

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by