How to solve a double complex ode using ode45 in MATLAB

Question

嘉杰程 2021-10-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1564996-how-to-solve-a-double-complex-ode-using-ode45-in-matlab

评论：嘉杰程 2021-10-17

The odes are:

ode(1) = diff(x,t) == (sqrt(1+x^2)*y)/(x*(t^2));
ode(2) = diff(y,t) == (x^3)*(t^2);
given x(0)=infty; y(0)=0.

I want to draw two picture to reflect x^3(u)-y relation and the u axis should be drawn by log(u) scale.

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

请先登录，再进行评论。

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

Answer 1

Walter Roberson 2021-10-16

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1564996-how-to-solve-a-double-complex-ode-using-ode45-in-matlab#answer_809816

在 MATLAB Online 中打开

The infinite boundary condition is a problem for symbolic solution: dsolve() rejects it.

The 0 initial time is a problem: you divide by time, so you generate a numeric infinity at the initial time, which numeric solvers cannot recover from.

If you set initial time above 0, and set the boundary condition to be finite, then you get a singularity. The time of singularity depends upon the boundary condition you set -- with the failure time being approximately 1 over the boundary condition.

syms x(t) y(t)

dx = diff(x)

dx(t) =

dy = diff(y)

dy(t) =

odes = [dx == (sqrt(1+x^2)*y)/(x*(t^2));

dy == (x^3)*(t^2)]

odes(t) =

ic = [x(0) == 1e8, y(0) == 0]

ic =

sol = dsolve(odes, ic)

Warning: Unable to find symbolic solution.

sol = [ empty sym ]

[eqs,vars] = reduceDifferentialOrder(odes,[x(t), y(t)])

eqs =

vars =

[M,F] = massMatrixForm(eqs,vars)

M =

F =

f = M\F

f =

odefun = odeFunction(f,vars)

odefun = function_handle with value:

@(t,in2)[(1.0./t.^2.*in2(2,:).*sqrt(in2(1,:).^2+1.0))./in2(1,:);t.^2.*in2(1,:).^3]

xy0 = double(rhs(ic))

xy0 = 1×2

100000000 0

tspan = [1e-50 100]

tspan = 1×2

0.0000 100.0000

[t, xy] = ode45(odefun, tspan, xy0);

Warning: Failure at t=2.575086e-08. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (5.293956e-23) at time t.

figure

semilogy(t, xy(:,1), 'k+-')

title('x')

figure

plot(t, xy(:,2), 'k*-')

title('y')

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

Walter Roberson 2021-10-17

在 MATLAB Online 中打开

syms x(t) y(t)

dx = diff(x)

dx(t) =

dy = diff(y)

dy(t) =

odes = [dx == (sqrt(1+x^2)*y)/(x*(t^2));

dy == (x^3)*(t^2)]

odes(t) =

ic = [x(0) == 1e8, y(0) == 0]

ic =

[eqs,vars] = reduceDifferentialOrder(odes,[x(t), y(t)])

eqs =

vars =

[M,F] = massMatrixForm(eqs,vars)

M =

F =

f = M\F

f =

odefun = odeFunction(f,vars)

odefun = function_handle with value:

@(t,in2)[(1.0./t.^2.*in2(2,:).*sqrt(in2(1,:).^2+1.0))./in2(1,:);t.^2.*in2(1,:).^3]

xy0 = double(rhs(ic))

xy0 = 1×2

100000000 0

tspan = [1e-50 100]

tspan = 1×2

0.0000 100.0000

[t, xy] = ode45(odefun, tspan, xy0);

Warning: Failure at t=2.575086e-08. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (5.293956e-23) at time t.

y = xy(:,2);

u = xy(:,1).^3;

ax = subplot(2,1,1)

ax =

Axes with properties: XLim: [0 1] YLim: [0 1] XScale: 'linear' YScale: 'linear' GridLineStyle: '-' Position: [0.1300 0.5838 0.7750 0.3412] Units: 'normalized' Show all properties

plot(ax, u, y);

xlabel(ax, 'u'); ylabel(ax, 'y')

ax.XScale = 'log';

ax.YScale = 'log';

ax = subplot(2,1,2)

ax =

Axes with properties: XLim: [0 1] YLim: [0 1] XScale: 'linear' YScale: 'linear' GridLineStyle: '-' Position: [0.1300 0.1100 0.7750 0.3412] Units: 'normalized' Show all properties

plot(ax, y, u);

xlabel(ax, 'y'); ylabel(ax, 'u')

ax.XScale = 'log';

ax.YScale = 'log';