Использование интерполяции вне определения функции в решении с использованием Рунге-Кутты 4-го порядка - PullRequest
Новая
       28

Использование интерполяции вне определения функции в решении с использованием Рунге-Кутты 4-го порядка

0 голосов
/ 08 апреля 2019

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

enter image description here
с enter image description here где enter image description here и z 2 = x 2 + (1 + a ) x 1 

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))];

Но это дает ошибку enter image description here

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

1 Ответ

2 голосов
/ 08 апреля 2019

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

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))];

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

...