2016-10-26 8 views
-1

УчитываяНахождение неизвестного обыкновенного дифференциального уравнения

d²x/dt² + a·dx/dt + 7.9·x³ = 3.2·sin(xt) 

с начальными условиями

x(0)  = +1.2 
dx/dt(0) = −3.3 
x(2.3) = −0.6 

Найти численно все возможные значения a, каждый из которых с точностью не менее 3-х цифр.

Есть ли какой-либо метод, кроме грубой силы для решения этой проблемы?

+0

У Matlab есть [ode solver] (https://www.mathworks.com/help/matlab/ref/ode45.html?requestedDomain=www.mathworks.com), если вы хотите решить это практически, или вы спрашивая академически? – code11

+0

Итак, я имел в виду численное использование Matlab, я знаю об ode45, но должен ли я взять все возможные значения a и найти x (0) для каждого, а затем взять пересечение с линией y = -0.6. Или, есть ли другое решение? –

+0

У Matlab также есть [символическое дифференциальное уравнение] (https://www.mathworks.com/help/symbolic/dsolve.html) решатель, который должен делать то, что вы хотите. Примерно на треть пути вниз показан пример, где они демонстрируют решение (для константы) уравнения второго порядка, подобное тому, которое вы предоставили. – code11

ответ

4

Насколько я вижу, решить эту проблему, как указано, невозможно.

Вот что я сделал. Я осуществил вашу проблему в достаточно общем виде:

%{ 

Find all 'a' for which 

d²x/dt² + a·dx/dt + 7.9·x³ - 3.2·sin(xt) = 0 

with initial conditions 

x(0)  = +1.2 
dx/dt(0) = −3.3 
x(2.3) = −0.6 

%} 

function odetest 

    % See how the function search_a(a) behaves around a = 0: 
    test_as = 0 : 0.1 : 10; 
    da = zeros(size(test_as)); 
    for ii = 1:numel(test_as)   
     da(ii) = search_a(test_as(ii)); end 

    figure(100), clf, hold on 
    plot(test_as, da) 
    axis tight 
    xlabel('a') 
    ylabel('|x(2.3) - 0.6|') 


    % Roughly cherry-pick some positive values, improve the estimate, and 
    % plot the solutions 

    opt = optimset('tolfun',1e-14, 'tolx',1e-12); 

    plot_x(fminsearch(@search_a, 0.0, opt), 1) 
    plot_x(fminsearch(@search_a, 1.4, opt), 2) 
    plot_x(fminsearch(@search_a, 3.2, opt), 3) 

    % Plot single solution 
    function plot_x(a,N) 

     [xt, t] = solve_ode(a); 

     figure(N), clf, hold on 
     plot(t,xt) 
     plot(2.3, -0.6, 'rx', 'markersize', 20) 
     title (['x(t) for a = ' num2str(a)]) 
     xlabel('t') 
     ylabel('x(t)') 
    end 
end 

% Solve the problem for a value a, and return the difference between the 
% actual value and desired value (-0.6) 
function da = search_a(a) 

    a_desired = -0.6; 

    xt = solve_ode(a);  
    da = abs(xt(end) - a_desired); 
end 


% Solve the problem for any given value of a 
function [xt, t] = solve_ode(a) 

    y0  = [1.2 -3.3]; 
    tfinal = 2.3; 

    opt = odeset('AbsTol',1e-12, 'RelTol',1e-6);  
    [t,yt] = ode45(@(y,t) odefun(y,t,a), [0 tfinal], y0, opt);  
    xt  = yt(:,1); % transform back to x(t) 
end 

% Most ODE solvers solve first-order systems. This is not a problem for a 
% second-order system, because if we make the transformation 
% 
% y(t) = [ x (t) 
%   x'(t) ] 
% 
% Then we can solve for 
% 
% y'(t) = [ x' (t) 
%    x''(t) ] <- the second-order homogeneous DE 
% 
function dydt = odefun(t,y,a)  
    dydt = [y(2) 
      -a*y(2) - 7.9*y(1)^3 + 3.2*sin(y(1)*t)];  
end 

Первая часть дала мне эту цифру:

solution for all positive a

Некоторые дальнейшие исследования показывают, что это растет только для больших a.

Эта цифра породило первоначальные оценки a = [0, 1.4, 3.2], который я улучшил через fminsearch() и нанесены решения:

x(t) for a ≈ 0.0 x(t) for a ≈ 1.4 x(t) for a ≈ 3.2

Так, что, вероятно, позволит вам сдать в выполнении домашних заданий:)

Однако, почему я говорю, что невозможно ответить на вопрос, как указано, потому что это то, что первый сюжет выглядит как отрицательноa:

solution for negative a

Колебательное поведение кажется продолжать до бесконечности, а расстояние между нулями, кажется, снижением непредсказуемым образом.

Теперь мои университетские дни давно позади, и я не так хорошо разбираюсь в теории ОДУ. Возможно, там есть шаблон, который просто не показывает из-за числовых проблем. Или, возможно, колебание останавливается после некоторого значения, никогда больше не возвращаться. Или, может быть, еще один ноль появляется на a = +1053462664212.25.

Я не могу доказать ничего из этого, я просто знаю, как это сделать; Остальное зависит от тебя.

+1

На графике должно быть '' | x (2.3) + 0.6 | ''или'' | x (2.3) - (-0.6) | ''. – LutzL

+0

Да! Я закончил тем же, и я проверил его до -6, колебания продолжают происходить. Это не домашнее задание, хотя :) –