Блог

Я написал код MATLAB для решения следующих систем дифференциальных уравнений.

с где и z₂ = x₂ (1 a)x₁

 a = 2;
k = 1 a;  
b = 3;
ca = 5;
cb = 2; 

theta1t = 0:.1:10;
theta1 = ca*normpdf(theta1t-5);

theta2t = 0:.1:10;
theta2 = cb*ones(1,101);

h = 0.05;
t = 1:h:10;

y = zeros(2,length(t));  

y(1,1) = 1;  % <-- The initial value of y  at time 1
y(2,1) = 0;  % <-- The initial value of y' at time 1

f = @(t,y) [y(2) interp1(theta1t,theta1,t,'spline')*y(1)*sin(y(2));
            interp1(theta2t,theta2,t,'spline')*(y(2)^2) y(1)-y(1)-y(1)-(1 a)*y(2)-k*(y(2) (1 a)*y(1))];

for i=1:(length(t)-1)  % At each step in the loop below, changed y(i) to y(:,i) to accommodate multi results
   k1 = f( t(i)      , y(:,i)         );
   k2 = f( t(i) 0.5*h, y(:,i) 0.5*h*k1);
   k3 = f( t(i) 0.5*h, y(:,i) 0.5*h*k2);
   k4 = f( t(i)     h, y(:,i)     h*k3);
   y(:,i 1) = y(:,i)   (1/6)*(k1   2*k2   2*k3   k4)*h; 
end 

plot(t,y(:,:),'r','LineWidth',2);
legend('RK4');
xlabel('Time')
ylabel('y') 
  

Теперь то, что я хочу сделать, это определить интерполяции / экстраполяции вне определения функции, например

 theta1_interp = interp1(theta1t,theta1,t,'spline');
theta2_interp = interp1(theta2t,theta2,t,'spline');
f = @(t,y) [y(2) theta1_interp*y(1)*sin(y(2));
           theta2_interp*(y(2)^2) y(1)-y(1)-y(1)-(1 a)*y(2)-k*(y(2) (1 a)*y(1))];
  

Но это выдает ошибку

Пожалуйста, предложите решение этой проблемы.

Обратите внимание, что в вашем исходном коде:

 f = @(t,y) [y(2) interp1(theta1t,theta1,t,'spline')*y(1)*sin(y(2));
            interp1(theta2t,theta2,t,'spline')*(y(2)^2) y(1)-y(1)-y(1)-(1 a)*y(2)-k*(y(2) (1 a)*y(1))];
  

вызов interp1 использует входную переменную t . t внутренняя часть этой анонимной функции отличается от t внешней, где она определяется как вектор.

Это означает, что, когда вы делаете

 theta1_interp = interp1(theta1t,theta1,t,'spline');
  

тогда theta1_interp это вектор, содержащий интерполированные значения для всех ваших t , а не только для одного. Один из способов обойти это — создать больше анонимных функций:

 theta1_interp = @(t) interp1(theta1t,theta1,t,'spline');
theta2_interp = @(t) interp1(theta2t,theta2,t,'spline');
f = @(t,y) [y(2) theta1_interp(t)*y(1)*sin(y(2));
            theta2_interp(t)*(y(2)^2) y(1)-y(1)-y(1)-(1 a)*y(2)-k*(y(2) (1 a)*y(1))];
  

Хотя на самом деле это никоим образом не улучшает ваш код по сравнению с оригиналом.