How to solve a double complex ode using ode45 in MATLAB

2021 10 16
2 个回答
回答已采纳
更新时间：2021 10 17
5 次查看（30 天）

0 个投票

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.

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

 采纳的回答

Walter Roberson
Walter Roberson 2021-10-16
Ran in:

1 个投票

Ran in:
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更早的评论

嘉杰 程
嘉杰 程 2021-10-16
编辑：嘉杰 程 2021-10-16
your method works, thanks, but there's another question is that could you draw an asymptotic line when x(u) approaches +infty, y=?
嘉杰 程
嘉杰 程 2021-10-16
and maybe you didn't notice that I want to plot y-u relation (u=x^3), can you do that?
嘉杰 程
嘉杰 程 2021-10-17
this is y-x relation, and how to print y-u relation (u=x^3)
Walter Roberson
Walter Roberson 2021-10-17
Ran in:
Ran in:
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';

请先登录，再进行评论。

更多回答（1 个）

Cris LaPierre
Cris LaPierre 2021-10-16

0 个投票

You can solve a symbolic ODE with initial/boundary conditions using dsolve

类别

在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

产品

标签

Community Treasure Hunt

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

Start Hunting!

Translated by