решение системы оды с помощью Matlab

У меня есть 9 уравнений с зависящим от времени коэффициентом g

% MY M file
function dy =tarak(t,y)
G= 3.16;
g =  0.1*exp(-((t-200)/90).^2);
dy=zeros(9,1);
dy(1)=-2*2*y(1)+2*G*y(5)+2*g*y(7);
dy(2)=2*y(1)-2*G*y(5);
dy(3)=2*y(1)-2*g*y(7);
dy(4)=-2*y(4)+g*y(9);
dy(5)=-2*y(5)+G*(y(2)-y(1))+g*y(8);
dy(6)=-2*y(6)-G*y(9);
dy(7)=-2*y(7)+g*(y(3)-y(1))+G*y(8);
dy(8)=-G*y(7)-g*y(5);
dy(9)=G*y(6)-g*y(4);

затем в командном окне:

[T,Y] = ode45(@tarak,[0 ,500],[0 0 1 0 0 0 0 0 0])

где коэффициент G = 3.16 и g = 0.1*exp(-((t-200)/90).^2) является зависящим от времени коэффициентом и временем t = 0:500; Исходное состояние [0 0 1 0 0 0 0 0 0].

Я получаю НЕПРАВИЛЬНЫЕ отрицательные значения для выходных данных y(1), y(2). Может кто-нибудь, пожалуйста, попробуйте решить приведенные выше уравнения с ode45, чтобы я мог сравнить результаты.


matlab ode differential-equations
person harman singh    schedule 01.04.2015    source источник
comment
Откуда вы взяли, что эти компоненты должны оставаться положительными? Разве это не зависит от y(5) и y(7)?   -  person Lutz Lehmann    schedule 01.04.2015
comment
привет, Лутцл, Y (1), Y (2), Y (3) - это вероятности, поэтому они не могут быть отрицательными, мой друг решил приведенные выше уравнения в C и понял это правильно. вы тоже получаете отрицательное значение?   -  person harman singh    schedule 01.04.2015
comment
Нет, см. ответ. Код был на питоне, но это не должно иметь значения. Вы должны опубликовать более полный код, чтобы увидеть, есть ли какие-либо легко различимые проблемы.   -  person Lutz Lehmann    schedule 01.04.2015
comment
дорогой Лутцль, Спасибо! ваши графики согласуются с результатами, которые я ищу. Вы использовали Odeint в python?   -  person harman singh    schedule 01.04.2015
comment
Lutzl, я отредактировал код в своем первом сообщении, чтобы дать полный код Matlab, пожалуйста, посмотрите   -  person harman singh    schedule 01.04.2015
comment
Нет, просто простая кулинарная реализация RK4. Я добавил код в ответ. Обратите внимание, что массивы Python отсчитываются от нуля, поэтому сдвиг индекса.   -  person Lutz Lehmann    schedule 01.04.2015
comment
Спасибо, сэр, за ценную информацию, я пройдусь по ней, что-то не так в моем коде Matlab? еще раз спасибо   -  person harman singh    schedule 01.04.2015
comment
Код Matlab в порядке, но вам нужно поиграть с допусками, потому что ваше первое состояние имеет порядок 1e-8, а допуск по умолчанию для ode45 равен 1e-6.   -  person remus    schedule 01.04.2015
comment
Нет, не видно ничего, что могло бы быть неправильным. Единственное, что немного странно, это точка в g=...).^2), так как t в этой формуле является скаляром, поэтому нет необходимости в векторизованных операциях. Но это не должно изменить числовые значения.   -  person Lutz Lehmann    schedule 01.04.2015
comment
Спасибо обоим @remus , @Lutzl .   -  person harman singh    schedule 01.04.2015

Ответы (2)


arrow_upward
1
arrow_downward

С помощью простого приложения RK4 я получаю эту картину

введите здесь описание изображения

приятный позитив, с одним странным начальным скачком в компоненте y(1). Но обратите внимание на масштаб, в целом y(1) довольно маленький. Кажется, что система на данный момент жесткая, поэтому у rk45 могут быть проблемы, лучше использовать неявный метод Рунге-Кутты.

И зум начальных колебаний

введите здесь описание изображения

Код Python

import numpy as np
import matplotlib.pyplot as plt
from math import exp

def dydt(t,y):
    dy = np.array(y);

    G = 3.16;
    g = 0.1*exp(-((t-200)/90)**2);

    dy[0]=-2*2*y[0]+2*G*y[4]+2*g*y[6];
    dy[1]=   2*y[0]-2*G*y[4];
    dy[2]=   2*y[0]-2*g*y[6];
    dy[3]=  -2*y[3]+  g*y[8];
    dy[4]=  -2*y[4]+  G*(y[1]-y[0])+g*y[7];
    dy[5]=  -2*y[5]-  G*y[8];
    dy[6]=  -2*y[6]+  g*(y[2]-y[0])+G*y[7];
    dy[7]=  -G*y[6]-  g*y[4];
    dy[8]=   G*y[5]-  g*y[3];
    return dy;

def RK4Step(f,x,y,h):
    k1=f(x      , y         )
    k2=f(x+0.5*h, y+0.5*h*k1)
    k3=f(x+0.5*h, y+0.5*h*k2)
    k4=f(x+    h, y+    h*k3)
    return (k1+2*(k2+k3)+k4)/6.0


t= np.linspace(0,500,200+1);
dt = t[1]-t[0];
y0=np.array([0, 0, 1, 0, 0, 0, 0, 0, 0]);

y = [y0]

for t0 in t[0:-1]:
    N=200;
    h = dt/N;
    for i in range(N):
        y0 = y0 + h*RK4Step(dydt,t0+i*h,y0,h);
    y.append(y0);

y = np.array(y);

plt.subplot(121);
plt.title("y(1)")
plt.plot(t,y[:,0],"b.--")
plt.subplot(122);
plt.title("y(2)")
plt.plot(t,y[:,1],"b-..")
plt.show()
person Lutz Lehmann    schedule 01.04.2015

arrow_upward
1
arrow_downward

И в Matlab:

options = odeset('AbsTol', 1e-12);
[T,Y] = ode45(@tarak, [0, 500], [0 0 1 0 0 0 0 0 0], options);

введите здесь описание изображения

person remus    schedule 01.04.2015
comment
@harmansingh Пожалуйста. Вы должны принять это как ответ, если он решит вашу проблему. - person remus; 01.04.2015