How to solve a double complex ode using ode45 in MATLAB

조회 수: 6 (최근 30일)
嘉杰 程
嘉杰 程 . 2021년 10월 16일
댓글: 嘉杰 程 . 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.

채택된 답변

Walter Roberson
Walter Roberson 2021년 10월 16일
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:
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.
semilogy(t, xy(:,1), 'k+-')
plot(t, xy(:,2), 'k*-')
  댓글 수: 5
嘉杰 程
嘉杰 程 2021년 10월 17일
thank you!

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

Cris LaPierre
Cris LaPierre 2021년 10월 16일
You can solve a symbolic ODE with initial/boundary conditions using dsolve


Help CenterFile 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