如何求解常微分方程组...初值问题....参数依赖于时间或自变量?说我有的等式
How solve a system of ordinary differential equation ..an initial value problem ....with parameters dependent on time or independent variable? say the equation I have
Dy(1)/dt=a(t)*y(1)+b(t)*y(2); Dy(2)/dt=-a(t)*y(3)+b(t)*y(1); Dy(3)/dt=a(t)*y(2);其中 a(t) 是一个向量,而 b(t) =c*a(t);其中 a 和 b 的值不以单调方式和每个时间步随时间变化.我尝试使用 这篇文章....但是当我应用相同的原理时...我收到错误消息
where a(t) is a vector and b(t) =c*a(t); where the value of a and b are changing with time not in monotone way and each time step. I tried to solve this using this post....but when I applied the same principle ...I got the error message
"使用 griddedInterpolant 时出错 点坐标不是按严格的单调顺序排列."
"Error using griddedInterpolant The point coordinates are not sequenced in strict monotonic order."
有人可以帮我吗?
推荐答案请阅读到最后,看看答案的第一部分或第二部分是否与您相关:
Please read until the end to see whether the first part or second part of the answer is relevant to you:
第 1 部分:首先创建一个 .m 文件,其中包含一个描述您的计算的函数以及将给出 a 和 b 的函数.例如:创建一个名为 fun_name.m 的文件,其中将包含以下代码:
Part 1: First create an .m file with a function that describe your calculation and functions that will give a and b. For example: create a file called fun_name.m that will contain the following code:
function Dy = fun_name(t,y) Dy=[ a(t)*y(1)+b(t)*y(2); ... -a(t)*y(3)+b(t)*y(1); ... a(t)*y(2)] ; end function fa=a(t); fa=cos(t); % or place whatever you want to place for a(t).. end function fb=b(t); fb=sin(t); % or place whatever you want to place for b(t).. end然后使用带有以下代码的第二个文件:
Then use a second file with the following code:
t_values=linspace(0,10,101); % the time vector you want to use, or use tspan type vector, [0 10] initial_cond=[1 ; 0 ; 0]; [tv,Yv]=ode45('fun_name',t_values,initial_cond); plot(tv,Yv(:,1),'+',tv,Yv(:,2),'x',tv,Yv(:,3),'o'); legend('y1','y2','y3');当然对于我写的 fun_name.m 情况,您不需要为 a(t) 和 b(t) 使用子函数,如果可能的话,您可以在 Dy 中使用显式函数形式(例如 cos(t) 等).
Of course for the fun_name.m case I wrote you need not use sub functions for a(t) and b(t), you can just use the explicit functional form in Dy if that is possible (like cos(t) etc).
第 2 部分:如果 a(t) 、 b(t) 只是您碰巧拥有的无法表达的数字的向量作为 t 的函数(如第 1 部分),那么您还需要有一个时间向量,每个时间向量都会发生,这当然可以与您将使用的时间相同ODE,但不必如此,只要插值有效即可.当它们具有不同的时间跨度或分辨率时,我将处理一般情况.然后您可以执行以下操作,创建 fun_name.m 文件:
Part 2: If a(t) , b(t) are just vectors of numbers you happen to have that cannot be expressed as a function of t (as in part 1), then you'll need to have also a time vector for which each of them happens, this can be of course the same time you'll use for the ODE, but it need not be, as long as an interpolation will work. I'll treat the general case, when they have different time spans or resolutions. Then you can do something of the following, create the fun_name.m file:
function Dy = fun_name(t, y, at, a, bt, b) a = interp1(at, a, t); % Interpolate the data set (at, a) at times t b = interp1(at, b, t); % Interpolate the data set (bt, b) at times t Dy=[ a*y(1)+b*y(2); ... -a*y(3)+b*y(1); ... a*y(2)] ;要使用它,请参阅以下脚本:
In order to use it, see the following script:
%generate bogus `a` ad `b` function vectors with different time vectors `at` and `bt` at= linspace(-1, 11, 74); % Generate t for a in a generic case where their time span and sampling can be different bt= linspace(-3, 33, 122); % Generate t for b a=rand(numel(at,1)); b=rand(numel(bt,1)); % or use those you have, but you also need to pass their time info... t_values=linspace(0,10,101); % the time vector you want to use initial_cond=[1 ; 0 ; 0]; [tv,Yv]= ode45(@(t,y) fun_name(t, y, at, a, bt, b), t_values, initial_cond); % plot(tv,Yv(:,1),'+',tv,Yv(:,2),'x',tv,Yv(:,3),'o'); legend('y1','y2','y3');