ode45 failing to produce a solution

I'm looking for some help with ode45 solver.

I have the following set of ODEs:

function Y=f(~,X)

n=1.3;

dF1deta=X(2);

dF2deta=n^(-1)*((X(2)^(2)+X(4)^(2))^((n-1)/2))^(-1)*((X(1)^(2)-X(3)^(2)+(-X(5))*X(2))*(1+(n-1)*(X(2)^(2)+X(4)^2)^(-1)*X(4)^(2))-(n-1)*X(2)*X(4)*(X(2)^(2)+X(4)^(2))^(-1)*(2*X(1)*X(3)+(-X(5))*X(4)));

dG1deta=X(4);

dG2deta=n^(-1)*((X(2)^(2)+X(4)^(2))^((n-1)/2))^(-1)*((2*X(1)*X(3)+(-X(5))*X(4))*(1+(n-1)*(X(2)^(2)+X(4)^2)^(-1)*X(2)^(2))-(n-1)*X(2)*X(4)*(X(2)^(2)+X(4)^(2))^(-1)*(X(1)^(2)-X(3)^(2)+(-X(5))*X(2)));

dH1deta=(3*n+1)/(n+1)*X(1);

Y = [dF1deta; dF2deta; dG1deta; dG2deta; dH1deta];

I need to solve this as such:

etaspan = [0 20];

Xinit = [0; 0; 0; 0; -0.73591509];

options = odeset('AbsTol', 1e-10, 'RelTol', 1e-8);

[eta,X]=ode45(@f,etaspan,Xinit,options);

figure;
hold all;
plot(eta,X(:,1)); 
plot(eta,X(:,3)); 
plot(eta,-X(:,5)); 
hold off;
xlabel('\eta')
hleg = legend('F','G','-H','Location','NorthWest');

In this case ode45 fails to produce a solution. I'm guessing this is because of the four initial conditions being set to zero? Is there another solver that will be able to handle this?

Thanks.

[Jan]

What does "ODE45 fails" explicitly mean? Do you get an error message or do the results differ from your expectations? Please explain such important details.

An absolute tolerance of 1e-10 is really very tiny. I cannot imagine any scientific field, where an ODE45 integration with such a tiny tolerance has any kind of power and meaning.

[pxg882]

Granted the tolerance is high, maybe that should be altered. Running the script produces a blank plot, I am assuming this is due to the initial values being set to zero. It would appear that since these values are zero at the initial state they remain at zero for all eta. Looking at an example solution to this problem I know this isn't the case though.

Answer 1

[Jan]

This does not concern your problem, but I like clean code, because it makes debugging much easier.

• 1/X is much cheaper than x^-1

n = 1.3;
XX = X .^ 2;
dF2deta = 1 / (n * ((XX(2) + XX(4))^((n-1)/2)) * ...
              ((XX(1)-XX(3) - X(5)*X(2)) * ...
               (1+(n-1) *(XX(2)+XX(4)) * XX(4)) - ...
               (n-1) * X(2)*X(4) / (XX(2)+XX(4)) * ...
              (2*X(1) * X(3) - X(5) * X(4))));

Without an editor for the matching of parenthesis the cleanup was not secure. So please test this, or better apply the shown simplifications by your own. Do you see the optical improvements? Usually theye are accompanied by a faster processing also. But this can still be improved.

Answer 2

[Jan]

You can use the debugger to find out, if your function replies Y value differing from 0. Simply set a breakpoint in a lines of your function "f" and look, what's going on.

Again I repeat, that in my opinion lines like these cannot be debugged by a human - at least not with a reasonable effort:

dG2deta=n^(-1)*((X(2)^(2)+X(4)^(2))^((n-1)/2))^(-1)*((2*X(1)*X(3)+(-X(5))*X(4))*(1+(n-1)*(X(2)^(2)+X(4)^2)^(-1)*X(2)^(2))-(n-1)*X(2)*X(4)*(X(2)^(2)+X(4)^(2))^(-1)*(X(1)^(2)-X(3)^(2)+(-X(5))*X(2)));

Simplify your code to simplify your life!

[pxg882]

O.k I've tried my best to simplify things a bit, this is what I have arrived at now:

function Y=critical(x,X)

n=1.3;

A=1/(X(2)^2+X(4)^2)^((n-1)/2);

B=(n-1)*1/(X(2)^2+X(4)^2)^2;

C=-(n-1)*X(2)*X(4)*1/(X(2)^2+X(4)^2)^2;

dF1deta=X(2);

dF2deta= (A / n)*((X(1)^2-X(3)^2-X(5)*X(2))*(1+B*X(4)^2)+C*(2*X(1)*X(3)-X(5)*X(4)));

dG1deta=X(4);

dG2deta= (A / n)*((2*X(1)*X(3)-X(5)*X(4))*(1+B*X(2)^2)+C*(X(1)^2-X(3)^2-X(5)*X(2)));

dH1deta=-2*X(1)-(1-n)/(n+1)*x*X(2);

Y = [dF1deta; dF2deta; dG1deta; dG2deta; dH1deta];

Does this make things easier to look at?

[EDITED, Jan, X^(2) -> X^2]

[Jan]

What is the connection between this function "Y=critical(x,X)" and the function in the original post "Y=f(~,X)"?

[pxg882]

That's my fault, should be the same thing!

[Jan]

Ok. Now what does "fails to produce a solution" mean? Do you get an error message of a trajectory, which stays on the zero value?

[pxg882]

Right o.k. I've been confusing myself here a bit as well (trying to do too many things at once!).

So the re-written function to be solved is this:

function Y=f(~,X)

n=1.3;

A=1/(X(2)^2+X(4)^2)^((n-1)/2);
B=(n-1)*1/(X(2)^2+X(4)^2)^2;
C=-(n-1)*X(2)*X(4)*1/(X(2)^2+X(4)^2)^2;

dF1deta=X(2);

dF2deta=(A/n)*((X(1)^2-X(3)^2-X(5)*X(2))*(1+B*X(4)^2)+C*(2*X(1)*X(3)-X(5)*X(4)));

dG1deta=X(4);

dG2deta=(A/n)*((2*X(1)*X(3)-X(5)*X(4))*(1+B*X(2)^2)+C*(X(1)^2-X(3)^2-X(5)*X(2)));

dH1deta=(3*n+1)/(n+1)*X(1);

Y = [dF1deta; dF2deta; dG1deta; dG2deta; dH1deta];

Using ode45 as such:

etaspan = [0 6];

Xinit = [0; 0; 0; 0; 0.7359];

options = odeset('AbsTol', 1e-6, 'RelTol', 1e-4);

[eta,X]=ode45(@f,etaspan,Xinit,options);

figure;
hold all;
plot(eta,X(:,1)); 
plot(eta,X(:,3)); 
plot(eta,-X(:,5)); 
hold off;
xlabel('\eta')
hleg = legend('F','G','-H','Location','NorthWest');

Just making sure things are clear.

Running the above script produces a blank plot. I'm assuming this is indicating that all the solutions are identically zero for all values of eta. From the known solution I know this isn't the case though.