Coupled Differential Equations - derivatives on both sides of equations

Question

razorfin 2017-8-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/354174-coupled-differential-equations-derivatives-on-both-sides-of-equations

评论： razorfin 2017-8-28

I am trying to solve two coupled differential equations, but I cannot seem to figure out how to get into a form that can be solved with an ode solver. Specifically, I can't get all the derivatives to the left hand side of the equations.

The equations are:

 dx/dt=[-1.005*x+0.0000265*dy/dt+0.01*y+f(t)*g1(x,t)]*1/(0.000053+g1(x,t))
 dy/dt=[-0.872*y+0.0000265*dx/dt+0.01*x-f(t)*g2(y,t)]*1/(0.000053+g2(y,t))

Where:

 f(t)=44100*cos(1.48-377*t)-3888*exp(-33.3*t)
 g1(x,t)=24800/(126+.55*abs(117*sin(377*t-1.48)-117*sin(-1.48)*exp(-t/0.03)-x)*240)^2
 g2(y,t)=24800/(126+.55*abs(117*sin(377*t-1.48)-117*sin(-1.48)*exp(-t/0.03)+y)*240)^2

Any help would be greatly appreciated.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Andrew Newell 2017-8-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/354174-coupled-differential-equations-derivatives-on-both-sides-of-equations#answer_279327

编辑：Andrew Newell 2017-8-27

在 MATLAB Online 中打开

Let's simplify things a little bit so we can better see what is going on. Let

Y = [x; y];
G1(t,Y) = 0.000053+g1(Y(1),t);
G2(t,Y) = 0.000053+g2(Y(2),t);

Then a little rearrangement gives

G1(t,Y)*(d/dt)Y(1) - 0.0000265*(d/dt)Y(2) = F1(t,Y)
G2(t,Y)*(d/dt)Y(2) - 0.0000265*(d/dt)Y(1) = F2(t,Y)

where F1(t,Y) and F2(t,Y) bundle all the remaining terms together. If we now define a "mass matrix"

M(t,Y) = [G1(t,Y) -0.0000265; -0.0000265 G2(2,Y)]

(see the odeset documentation) and a function

F(t,Y) = [F1(Y,t); F2(Y,t)];

then we can define

options = odeset('Mass',M);

and solve something like

[t,y] = ode45(F,tspan,y0)

(see the ode45 documentation). You'll have to define functions F(t,Y) and M(t,Y) to reproduce the above steps. Warning: I have been using a mixture of MATLAB and regular math to keep it simple; you can sort out the details.

3 个评论
显示 1更早的评论隐藏 1更早的评论

razorfin 2017-8-27

在 MATLAB Online 中打开

Thank you very much for helping to set me on the right track. I see how the mass matrix can get the equations into the proper form now. I've spent a day reading examples and the literature and trying to code a solution, but I'm just not quite there yet. Perhaps there is something obvious that I'm doing wrong?

I created one file called "Fn.m" and used the following code (I added some more sig figs than in my original question):

 function Ydot=Fn(Y,t)
   Ydot=[(-Y(1)*.995+0.01*Y(2)+H(t)*G1(t,Y(1)))/(G1(t,Y(1))+2.652582E-5); (-Y(2)*.862+0.01*Y(1)+-H(t)*G2(t,Y(2)))/(G2(t,Y(2))+2.652582E-5)];
 end
 function G1OUT=G1(t,Y1)
   G1OUT=24873.198/(126.313+.54855621*abs(117.85113*sin(377*t-1.48260626)+117.393134086*exp(-t/0.03)-Y1)*240)^2;
 end
function G2OUT=G2(t,Y2)
 G2OUT=24873.198/(126.313+.54855621*abs(117.85113*sin(377*t-1.48260626)+117.393134086*exp(-t/0.03)+Y2)*240)^2;
end
 function H1=H(t)
   H1=44429.9*cos(1.48261-377*t)-3913.1*exp(-33.3333*t);
 end

Then I created another file called "Mass1.m" and added the following code (notice G1 and G2 are duplicated again):

 function M=Mass1(t,Y)
    M=[5.30504E-5+G1(t,Y(1)), -2.6525E-5; -2.6525E-5, 5.30504E-5+G2(t,Y(2))];
  end
 function G1OUT=G1(t,Y1)
    G1OUT=24873.198/(126.313+.54855621*abs(117.85113*sin(377*t-1.48260626)+117.393134086*exp(-t/0.03)-Y1)*240)^2;
 end
function G2OUT=G2(t,Y2)
  G2OUT=24873.198/(126.313+.54855621*abs(117.85113*sin(377*t-1.48260626)+117.393134086*exp(-t/0.03)+Y2)*240)^2;
end

Then I tried to solve in the command windows with these commands:

 t=linspace(0,0.1,10000)';
 options=odeset('Mass',@Mass1)
 [t,y]=ode45("Fn",t,[0;0],options);

The program complains that Y(2) is out of bounds, I think it must be something simple but I haven't found too many other mass matrix examples to be able to see how this is supposed to be done. Thank you.

Andrew Newell 2017-8-28

在 MATLAB Online 中打开

I made two modifications to your code. The first is to change the definition of Fn to

function Ydot=Fn(t,Y)

and the second was to change the final command to

[t,y]=ode45(@Fn,t,[0;0],options);

Then it ran fine.

razorfin 2017-8-28

Thank you, this has been a huge help to my understanding. Have a great day!

请先登录，再进行评论。

Coupled Differential Equations - derivatives on both sides of equations

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

3 个评论
显示 1更早的评论隐藏 1更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Coupled Differential Equations - derivatives on both sides of equations

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

3 个评论 显示 1更早的评论隐藏 1更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论