Please help solve a system of differential equation

1 次查看（过去 30 天）

Jayesh Jawandhia 2023-3-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1929020-please-help-solve-a-system-of-differential-equation

评论： Torsten 2023-4-4

I have a system of differential equation in x,y with boundary conditions. I want to solve it analytically and numerically. The analytical solutions coming in the form of error functions. Can someone please help me find analytical and numerical solution for this, I have been trying for weeks now, really need to solve this for my thesis. Please help.

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

请先登录，再进行评论。

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

回答（1 个）

Alan Stevens 2023-3-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1929020-please-help-solve-a-system-of-differential-equation#answer_1193540

在 MATLAB Online 中打开

Here's a numerical approach:

First manipulate the equations

tspan = [0, 120];

ic = [80000, 0.001];

[t, xy] = ode45(@rate, tspan, ic);

x = xy(:,1); y = xy(:,2);

subplot(2,1,1)

plot(t,x),grid

xlabel('t'), ylabel('x')

subplot(2,1,2)

plot(t,y),grid

xlabel('t'), ylabel('y')

function dxydt = rate(~,xy)

a = 0.1; b = 0.2; c = -0.1; % Replace with desired values

x = xy(1); y = xy(2);

dxydt = [a-b+c*x*y; % eqn (3)

c*y*(1-y) - a*y/x]; % eqn (5)

end

19 个评论
显示 17更早的评论隐藏 17更早的评论

Alan Stevens 2023-3-16

在 MATLAB Online 中打开

This data looks more like it's related to the logistic function than the error function. The following shows a purely empirical fit:

t = 5:5:120;

x = [15.6700 18.0000 19.0000 19.3300 19.5000 19.6700 19.3300 ...

18.7700 18.6000 18.2000 17.9300 17.9300 17.6700 17.3300 ...

17.1700 16.8300 16.8300 16.1700 15.5000 15.5000 15.8300 ...

15.6700 15.6700 15.6700];

y = [53.4467 61.4233 63.9467 64.6933 65.6633 65.7567 65.6100 ...

65.3000 65.6133 65.9367 66.2833 66.0300 67.0500 65.9600 ...

66.3733 66.7867 66.5233 66.3533 66.6567 66.6400 66.4267 ...

67.3067 67.0300 67.3133];

a = 67; b = 0.2; c = 0.01; abc = [a, b, c];

abc = fminsearch(@(abc) fny(abc, t, y), abc);

a = abc(1); b = abc(2); c = abc(3);

Y = a.*exp(-c*t)./(1+exp(-b*t));

A = 20; B = 0.2; C = 0.01; ABC = [A, B, C];

ABC = fminsearch(@(ABC) fnx(ABC, t, x), ABC);

A = ABC(1);B = ABC(2); C = ABC(3);

X = A.*exp(C*t)./(1+exp(-B*t));

plot(t,x,'ko',t,y,'rs',t,Y,'k',t,X), grid

xlabel('t'),ylabel('x and y')

legend('x','y','yfit','xfit')

axis([0 120 0 90])

hold on

function Fy = fny(abc, t, y)

a = abc(1); b = abc(2); c = abc(3);

Y = a.*exp(-c*t)./(1+exp(-b*t));

d = Y-y;

Fy = sum(norm(d));

end

function Fx = fnx(ABC, t, x)

A = ABC(1); B = ABC(2); C = ABC(3);

X = A.*exp(C*t)./(1+exp(-B*t));

d = X - x;

Fx = sum(norm(d));

end

Alan Stevens 2023-3-19

在 MATLAB Online 中打开

@Jayesh Jawandhia Ok. Here is the code I put together to have a quick look. Change the values of a, b and c to your heart's content! I'm not convinced that your odes describe your data, but more than happy to be proved wrong. If you are successful, post your results here. Good luck..

t = 0:300:7200;
x = [15.0000   17.2300   18.1900   18.5000   18.6600   18.8200   18.5000 ...
     17.9600   17.8000   17.4200   17.1600   17.1600   16.9100   16.5900 ...
     16.4300   16.1100   16.1100   15.4700   14.8400   14.8400   15.1500 ...
     15.0000   15.0000   15.0000];
 
 y = [0.5345    0.6142    0.6395    0.6469    0.6566    0.6576    0.6561 ...
      0.6530    0.6561    0.6594    0.6628    0.6603    0.6705    0.6596 ...
      0.6637    0.6679    0.6652    0.6635    0.6666    0.6664    0.6643 ...
      0.6731    0.6703    0.6731];
 
 % Add initial values to data
 x = [12.5 x];
 y = [0.001 y];
 
 % Initial guesses for a, b and c
 %      [a, b, c]
  abc = [1, 1, 1];
 
  % Try fminsearch to estimate better values of a, b and c
  ABC = fminsearch(@(abc) compare(abc, t, x, y), abc);
  
  disp(ABC)  % Display resulting values
  
  % Numerically integrate the odes with resulting values
  [t, XY] = runode(ABC,t);
  
  % Extract fitted values and plot fits vs data
  X = XY(:,1);  Y = XY(:,2);
  subplot(2,1,1)
  plot(t, X, 'r', t, x, 'ro'),grid
  ylabel('x')
  axis([0 7200 0 20])
  subplot(2,1,2)
  plot(t, Y, 'r', t, y, 'ro'),grid
  ylabel('y')
  axis([0 7200 0 1])
 
  % Calculate sum of squared residuals
  function F = compare(abc, t, x, y)
  
           [~, XY] = runode(abc, t);
           X = XY(:,1); Y = XY(:,2);
           dx = x - X;
           dy = y - Y;
           F = norm(dx)+norm(dy);
           
  end
  
  % Call the odes with current values of a, b and c
  function [t, XY] = runode(abc, t)
  
  ic = [12.5, 0.001];
  [t, XY] = ode45(@(t,xy) odefn(t, xy, abc), t, ic);
  end
  
% ODE equations
 function dxydt = odefn(~,xy, abc)
           a = abc(1); b = abc(2); c = abc(3);
           x = xy(1); y = xy(2);
           dxydt = [a-b+c*x*y;
                    c*y*(1-y)-a*y/x];
 
 end

Jayesh Jawandhia 2023-4-4

Hey @Alan Stevens

Hope you are doing well. I have been reading up on these and your solution has been really helpful. I just had a follow-up. We are trying to find a better fit by making a and b as a function of t. I haven't been able to figure out how do I add these to my ode45 solver in Matlab just like what you did here.

Can you please help with this? For simplicity, I am assuming both a and b second order in t, i.e.:

a = mt^2+nt+r

b = pt^2+qt+s

Really thankful for all your help!

Torsten 2023-4-4

在 MATLAB Online 中打开

What's the problem ?

function dxydt = rate(t,xy);
  m = ...;
  n = ...;
  r = ...;
  p = ...;
  q = ...;
  s = ...;
  a = m*t^2+n*t+r;
  b = p*t^2+q*t+s;
  x = xy(1);
  y = xy(2);
  dxydt = [a-b+c*x*y;c*y*(1-y)-a*y/x];
end

请先登录，再进行评论。

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

类别

MATLAB Mathematics Numerical Integration and Differential Equations

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

产品

MATLAB

Community Treasure Hunt

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

Start Hunting!

Translated by

Please help solve a system of differential equation

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

回答（1 个）

19 个评论
显示 17更早的评论隐藏 17更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

Please help solve a system of differential equation

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

回答（1 个）

19 个评论 显示 17更早的评论隐藏 17更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

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

19 个评论
显示 17更早的评论隐藏 17更早的评论