the function ode45 isn't accurate

when i use ode45 to solve two differential equation simultaneously it works well when i use small constants in the equation but when i use constant in the range 10^6 it doesn't work well and the answer isn't correct , i also try to use other ode functions but i get the same problem

can any one help me , or suggest another method to solve stiff differential equations correctly ?

[John D'Errico]

You use a stiff solver to solve stiff problems. How can you possibly be surprised at this?

[Jonel Ortiz]

Wow, take it easy. Maybe he/just doesn't know about the stiff solvers.

Answer 1

[Star Strider]

Use ode15s. If it doesn’t produce the results you want, search through the other solvers. In my experience, if ode45 nas problems, ode15s usually succeeds. Also, avoid discontinuities if at all possible. The ODE solvers don’t do well with non-differentiable functions.

[mai]

i try also ode15s but it doesn't work correctly , is there any another method or program to solve this types of problem correctly ?

[Star Strider]

If your ODE has discontinuities or singularities, all solvers are going to have problems. Without knowing the system you are integrating, it’s impossible to suggest other solutions.

[mai]

here is the two equations i try to solve (in the image ):

where b1,b2,k are constant and i^2=-1

when i set the value of b1,b2 in the range of 10^6 b1=5.84004*10^6, and b2=5.84*10^6 , k=0.1 initial values also x1(0)=0 ,x2(0)=1 , z=[0 10], the solver takes long time and the answer is not correct

[Star Strider]

I could not get it to integrate at all with ode15s (R2014a) with this code:

b1 = 5.84E+6;
b2 = b1;
k = 0.1;
ODESys = @(z,x) [-1i.*b1*x(1) - 1i*k*x(2);  -1i.*b2*x(2) - 1i*k*x(1)];
zv = linspace(1, 10, 25)
[z2, X1X2] = ode15s(ODESys, zv, [0 1]);

it stopped with:

Warning: Failure at t=1.260800e+00.  Unable to meet integration tolerances without reducing the step size below the
smallest value allowed (3.552714e-15) at time t. 
> In ode15s at 668

so I suspect you may have erred in coding it.

However solving it analytically with the Symbolic Math Toolbox, it behaved quite well:

x1 = @(X10,X20,b1,b2,k,z)1.0./sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(X10.*exp(z.*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(-1.0./2.0)).*exp(b1.*z.*-5.0e-1i).*exp(b2.*z.*-5.0e-1i).*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(1.0./2.0)+X10.*exp(z.*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(1.0./2.0)).*exp(b1.*z.*-5.0e-1i).*exp(b2.*z.*-5.0e-1i).*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(1.0./2.0)+X10.*b1.*exp(z.*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(-1.0./2.0)).*exp(b1.*z.*-5.0e-1i).*exp(b2.*z.*-5.0e-1i).*5.0e-1i-X10.*b1.*exp(z.*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(1.0./2.0)).*exp(b1.*z.*-5.0e-1i).*exp(b2.*z.*-5.0e-1i).*5.0e-1i-X10.*b2.*exp(z.*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(-1.0./2.0)).*exp(b1.*z.*-5.0e-1i).*exp(b2.*z.*-5.0e-1i).*5.0e-1i+X10.*b2.*exp(z.*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(1.0./2.0)).*exp(b1.*z.*-5.0e-1i).*exp(b2.*z.*-5.0e-1i).*5.0e-1i+X20.*k.*exp(z.*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(-1.0./2.0)).*exp(b1.*z.*-5.0e-1i).*exp(b2.*z.*-5.0e-1i).*1i-X20.*k.*exp(z.*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(1.0./2.0)).*exp(b1.*z.*-5.0e-1i).*exp(b2.*z.*-5.0e-1i).*1i)

x2 = @(X10,X20,b1,b2,k,z)X20.*exp(z.*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(-1.0./2.0)).*exp(b1.*z.*-1i).*exp(b2.*z.*-1i).*exp(b1.*z.*5.0e-1i).*exp(b2.*z.*5.0e-1i).*(1.0./2.0)+X20.*exp(z.*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(1.0./2.0)).*exp(b1.*z.*-1i).*exp(b2.*z.*-1i).*exp(b1.*z.*5.0e-1i).*exp(b2.*z.*5.0e-1i).*(1.0./2.0)-X20.*b1.*exp(z.*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(-1.0./2.0)).*exp(b1.*z.*-1i).*exp(b2.*z.*-1i).*exp(b1.*z.*5.0e-1i).*exp(b2.*z.*5.0e-1i).*1.0./sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*5.0e-1i+X20.*b1.*exp(z.*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(1.0./2.0)).*exp(b1.*z.*-1i).*exp(b2.*z.*-1i).*exp(b1.*z.*5.0e-1i).*exp(b2.*z.*5.0e-1i).*1.0./sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*5.0e-1i+X20.*b2.*exp(z.*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(-1.0./2.0)).*exp(b1.*z.*-1i).*exp(b2.*z.*-1i).*exp(b1.*z.*5.0e-1i).*exp(b2.*z.*5.0e-1i).*1.0./sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*5.0e-1i-X20.*b2.*exp(z.*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(1.0./2.0)).*exp(b1.*z.*-1i).*exp(b2.*z.*-1i).*exp(b1.*z.*5.0e-1i).*exp(b2.*z.*5.0e-1i).*1.0./sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*5.0e-1i+X10.*k.*exp(z.*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(-1.0./2.0)).*exp(b1.*z.*-1i).*exp(b2.*z.*-1i).*exp(b1.*z.*5.0e-1i).*exp(b2.*z.*5.0e-1i).*1.0./sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*1i-X10.*k.*exp(z.*sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*(1.0./2.0)).*exp(b1.*z.*-1i).*exp(b2.*z.*-1i).*exp(b1.*z.*5.0e-1i).*exp(b2.*z.*5.0e-1i).*1.0./sqrt(b1.*b2.*2.0-b1.^2-b2.^2-k.^2.*4.0).*1i

z1 = linspace(0,10);
fx1 = x1(0.0, 1.0, b1, b2, k, z1);
fx2 = x2(0.0, 1.0, b1, b2, k, z1);