求一个未知的常微分方程
Finding an unknown ordinary differential equation
给定
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 位有效数字。
除了暴力破解还有什么方法可以解决吗?
据我所知,无法按照说明解决此问题。
这是我所做的。我以一种相当通用的方式解决了你的问题:
%{
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
第一部分给了我这个数字:
一些进一步的调查表明,这只会随着 a
的增大而增长。
这个数字产生了初步估计 a = [0, 1.4, 3.2]
,我通过 fminsearch()
对其进行了改进并绘制了以下解决方案:
所以,这可能会让你交作业 :)
然而,为什么我说不可能回答所陈述的问题,是因为这是第一个情节 negative a
:
振荡行为似乎无限期地持续下去,零点之间的间距似乎以一种不可预测的方式减小。
现在,我的大学时代已经过去很久了,我对 ODE 理论也不是那么精通了。也许 是 的一种模式,只是由于数值问题而没有显示出来。或者振荡可能会在某个值后停止,再也不会 return。或者也许另一个零出现在 a = +1053462664212.25
。
我无法证明这些事情,我只知道如何暴力破解;剩下的由你决定。
给定
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 位有效数字。
除了暴力破解还有什么方法可以解决吗?
据我所知,无法按照说明解决此问题。
这是我所做的。我以一种相当通用的方式解决了你的问题:
%{
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
第一部分给了我这个数字:
一些进一步的调查表明,这只会随着 a
的增大而增长。
这个数字产生了初步估计 a = [0, 1.4, 3.2]
,我通过 fminsearch()
对其进行了改进并绘制了以下解决方案:
所以,这可能会让你交作业 :)
然而,为什么我说不可能回答所陈述的问题,是因为这是第一个情节 negative a
:
振荡行为似乎无限期地持续下去,零点之间的间距似乎以一种不可预测的方式减小。
现在,我的大学时代已经过去很久了,我对 ODE 理论也不是那么精通了。也许 是 的一种模式,只是由于数值问题而没有显示出来。或者振荡可能会在某个值后停止,再也不会 return。或者也许另一个零出现在 a = +1053462664212.25
。
我无法证明这些事情,我只知道如何暴力破解;剩下的由你决定。