the function ode45 isn't accurate
40 views (last 30 days)
Show older comments
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 ?
2 Comments
John D'Errico
on 25 Jun 2014
You use a stiff solver to solve stiff problems. How can you possibly be surprised at this?
Answers (1)
Star Strider
on 25 Jun 2014
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.
4 Comments
Star Strider
on 29 Jun 2014
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);
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!