Using numerical methods to plot solution to first-order nonlinear differential equation

Question

Aleem Andrew 2021-11-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1580514-using-numerical-methods-to-plot-solution-to-first-order-nonlinear-differential-equation

编辑： David Goodmanson 2021-11-17

I have a question about plotting x(t), the solution to the following differential equation knowing that dx/dt equals the expression below. The value of x is 0 at t = 0.

syms x
dxdt = -(1.0*(6.84e+45*x^2 + 5.24e+32*x - 2.49e+42))/(2.47e+39*x + 7.12e+37)

I want to plot the solution of this first-order nonlinear differential equation. The analytical solution involves complex numbers so that's not relevant but Matlab can solve the equation using numerical methods and plot it. Can anyone please suggest how to do this? Thank you.

2 个评论
显示无隐藏无

Chris 2021-11-7

I'm not sure how to verify, but does the following seem correct?

fun = @(t,x) -(1.0.*(6.84e+45.*x.^2 + 5.24e+32.*x - 2.49e+42))./(2.47e+39.*x + 7.12e+37);

tspan = [0,5e-6];

x0 = 0;

[t,y] = ode45(fun, tspan,x0);

figure;plot(t,y)

Sahil Jain 2021-11-16

The solution given by @Chris seems correct. For more information, you can refer to the documentation page for "ode45".

请先登录，再进行评论。

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

Answer 1

David Goodmanson 2021-11-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1580514-using-numerical-methods-to-plot-solution-to-first-order-nonlinear-differential-equation#answer_833734

编辑：David Goodmanson 2021-11-17

Hello Aleem,

You are saying that the analytical solution is not relevant because it involves complex numbers, but that's not the case. We have

dx/dt = f(x) --> t = Integral 1/f(x) dx + const

In code,

syms x
t = int(1/(-(1.0*(6.84e+45*x^2 + 5.24e+32*x - 2.49e+42))/(2.47e+39*x + 7.12e+37)))
t = vpa(t,6)
    t = 9.22306e-8*log(x + 0.0190797) - 4.53342e-7*log(x - 0.0190797)

For x near 0, t is complex because the argument of the second log term is negative. But

log(-a) = log(+a) + i*pi,

so a solution that is just as good can be obtained by subtracting off the constant i*pi, resulting in a reversed sign for of the argument of the second term. Then the function is real for 0 < x < 0.0190797. So

fun2 = @(x) 9.22306e-8*log(x + 0.0190797) - 4.53342e-7*log(0.0190797 - x);
y = linspace(0,.0190796,1000);  % pick upper limit for y that is just less than 0.0190797 
t = fun2(y)-fun2(0);          % subtract a constant so that t meets the boundary condition
% analytic solution
ya = y;         
ta = t;
% compare with Chris's solution 
fun = @(t,x) -(1.0.*(6.84e+45.*x.^2 + 5.24e+32.*x - 2.49e+42))./(2.47e+39.*x + 7.12e+37);
tspan = [0,5e-6];
x0 = 0;
[t,y] = ode45(fun, tspan,x0);
figure(1)
plot(t,y,ta,ya+.0002)   % add a small offset, otherwise they overlay
grid on