Search code examples
matlabsystemodedsolve

Error in dsolve when variable is multiplied by a constant (R2011a)

I am trying to use the following code:

ode1='D2y1=-1256.4*y1-5*Dy1+255.1*y2+182.781';   
ode2='D2y2=-151.5*y2-5*Dy2+255.1*y1-14.0459';   
CI='y1(0)=2,y2(0)=-2,Dy1(0)=0,Dy2(0)=0';

sol=dsolve(ode1,ode2,CI,'t');
sol.y1
sol.y2

and matlab returns an error:

??? Error using ==> mupadengine.mupadengine>mupadengine.feval at 144 MuPAD error: Error: Division by zero [_invert];

during evaluation of 'stdlib::normalNoExpand'

Error in ==> dsolve>mupadDsolve at 215 T = feval(symengine,'symobj::dsolve',sys,x,ignoreConstraints);

Error in ==> dsolve at 96 sol = mupadDsolve(ignoreConstraints,varargin{1:narg});

Error in ==> maglevsol at 7 sol=dsolve(ode1,ode2,CI,'t');

However, if i run this

ode1='D2y1=-y1-5*Dy1+255.1*y2+182.781';   
ode2='D2y2=-y2-5*Dy2+255.1*y1-14.0459';   
CI='y1(0)=2,y2(0)=-2,Dy1(0)=0,Dy2(0)=0';

sol=dsolve(ode1,ode2,CI,'t');
sol.y1
sol.y2

I get no errors at all. What's wrong? My version is r2011a

Solution

  • It seems Matlab is not able to solve this, due to issue with IC. (non numerical value for derivative at t=0 given from trying a numerical solver). I could not find a way around it (Matlab 2015a hangs actually on this, for analytical solution).

    But Maple is able to solve it (Maple is very good in ODE's). Here is the solution and the plots (it is also a good idea to use symbolic for all those hardcoded numerical values and then use subs at the end. No need to have magic numbers all over the place). You might want to report this to Mathworks

    restart;
ode1 := diff(y1(t),t$2) = a*y1(t) + b*diff(y1(t),t)+ c*y2(t) + d;
ode2 := diff(y2(t),t$2)= e*y2(t) + b*diff(y2(t),t)+ c*y1(t) + f;
ic := y1(0) = 2, y2(0) = -2, D(y1)(0) = 0, D(y2)(0) = 0;
sol:=dsolve({ode1, ode2, ic}, {y1(t), y2(t)}):
sol:=subs( {a = -1256.4, b = -5, c= 182.781, e = -151.5, f = -14.0459,d=182.781},sol);


sol := {y1(t) = (-.1922320350+0.5088536230e-1*I)*
exp((-2.500000000+9.444367980*I)*t)+(1.096005949-0.7581365270e-1*I)*
exp((-2.500000000+36.14144316*I)*t)+(-.1922320331-0.5088536185e-1*I)*
exp((-2.500000000-9.444367980*I)*t)+(1.096005952+0.7581365750e-1*I)*
exp((-2.500000000-36.14144316*I)*t)+.1924521734, 

y2(t) = (-.8748433300+.2315780499*I)*exp((-2.500000000+9.444367980*I)*t)+
(-.2408287839+0.1665876887e-1*I)*exp((-2.500000000+36.14144316*I)*t)+
(-.8748433212-.2315780478*I)*exp((-2.500000000-9.444367980*I)*t)+
(-.2408287844-0.1665876992e1*I)*
exp((-2.500000000-36.14144316*I)*t)+.2313442207}

    Here is the solution plots of y1(t) and y2(t) for 2 seconds

    Mathematica graphics

    Mathematica graphics