I have to solve this problem

$$\left\{\begin{matrix} y''=4y'+\frac{1}{h}-7y\\ y(1)=1\\ y'(1)=1 \end{matrix}\right. \text{where} \; h=5y'-2y$$

as it were a DAE (I know I could just substitute $h$ into the equation, but this is just an example, because in reality the problem I have to solve is a DAE and more complex than this). When I use ode45 and treat the problem as a second order differential equation, the graph $t$ Vs. $y$ is

enter image description here

but when I treat it as a DAE, the graph is completely different and I do not understand why. Here is my code:

ode45 second order differential equation

function yp = dae_normale(t,y)
yp = zeros(2,1);
yp(1) = y(2);
yp(2) = 4*y(2) + 1/(5*y(2) - 2*y(1) ) - 7*y(1);

ode45 second order differential equation run

[t,y] = ode45('dae_normale',[1,5],[1,1]);

DAE ode15s

function out = dae(t,y)
out = [y(2)
   4*y(2) + 1/y(3) - 7*y(1)
   y(3) - 5*y(2) + 2*y(1) ];

DAE ode15s run

y0 = [1; 1; 3];
M = [1 0 0; 0 1 0; 0 0 0];
options = odeset('Mass',M);
[t,y] = ode15s(@dae,[1 5],y0,options);

ode45 is an explicit (Runge-Kutta) ODE solver and, therefore, is only conditionally stable. Normally, ode45 will automatically reduce the step size to maintain stability. But, as discussed in this article, Error Control Matters, sometimes it does not because the error tolerances are too large. Try modifying your code as follows to reduce the error tolerances:

opts=odeset('RelTol', 1e-4, 'AbsTol', 1e-7);
[t,y] = ode45('dae_normale',[1,5],[1,1], opts);

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.