the function ode45 isn't accurate

20 次查看（过去 30 天）
mai
mai 2014-6-25
评论： Jonel Ortiz 2017-11-14
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 个评论
John D'Errico
John D'Errico 2014-6-25
You use a stiff solver to solve stiff problems. How can you possibly be surprised at this?
Jonel Ortiz
Jonel Ortiz 2017-11-14
Wow, take it easy. Maybe he/just doesn't know about the stiff solvers.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Star Strider
Star Strider 2014-6-25
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 个评论
显示 2更早的评论
mai
mai 2014-6-25
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
Star Strider 2014-6-29

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);

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by