Mathematica - проблемы с NDSolve

[Zepp]

Доброго времени суток! В неравном бою с Mathematica 8.0 не могу численно решить следующую систему из 2х ОДУ с параметром V (V~t как частный случай):

k1 = 6000;
k2 = 60000;
Ix = 0.2;
Iz = 0.05;
Ixz = 0.007;
Xc = 0.025;
Zc = 0.1;
ro = 0.7;
l = 0.3;

With[{a11 = k1, a22 = k2, c11 = Ix, c22 = Iz, c12 = -Ixz, c21 = -Ixz, 
  V = 350 + t, M = V/330, 
  k = (0.25 + 0.07*M*HeavisideTheta[3 - M] + 
     0.21*HeavisideTheta[M - 3]), K1 = 2*ro/(sqrt (M^2 - 1)), 
  d011 = (0.25 - 0.5*z)*z^2, 
  d012 = ((0.15 - 0.5*z)/(0.25 - 0.5*z) + ((k*z)/2 - z/2 - k/4 + 
         3/20)/(0.25 - 0.5*z))*(0.25 - 0.5*z)^2*z, 
  d022 = (((0.15 - 0.5*z)/(0.25 - 0.5*z) + ((k*z)/2 - z/2 - k/4 + 
            3/20)/(0.25 - 0.5*z))^2 + 1/12)*(0.25 - 0.5*z)^3, 
  d11 = K1*NIntegrate[d011, {z, 0, l}], 
  d12 = K1*NIntegrate[d012, {z, 0, l}], d21 = d12, 
  d22 = K1*NIntegrate[d022, {z, 0, l}], b012 = (0.25 - 0.5*z)*z, 
  b12 = -K1*NIntegrate[b012, {z, 0, l}], 
  b022 = ((0.15 - 0.5*z)/(0.25 - 0.5*z) + ((k*z)/2 - z/2 - k/4 + 
         3/20)/(0.25 - 0.5*z))*(0.25 - 0.5*z)^2, 
  b22 = -K1*NIntegrate[b022, {z, 0, l}]}
 , sol = NDSolve[{c11*q''[t] + d11*V*q'[t] + a11*q[t] + c12*r''[t] + 
      d12*V*r'[t] + b11*V^2*q[t] + b12*V^2*r[t] == 10 sin (t), 
    c21*q''[t] + d21*V*q'[t] + a22*r[t] + c22*r''[t] + d22*V*r'[t] + 
      b21*V^2*q[t] + b22*V^2*r[t] == 10 sin (t), q[0] == 4, 
    q'[0] == 6, r[0] == -17, r'[0] == -10}, {q, r}, {t, 0, 100}]]

Plot[Evaluate[{q[t], r[t]}], {t, 0, 100}]

Пожалуйста, укажите на ошибки или киньтесь годным работающим примером аналогичной задачи

[Vince Diesel]

Математика различает большие и маленькие буквы. И не будет численно решать, если в уравнениях есть переменные без численного значения.

[FFFF]

Насколько я понял, в параметр V входит переменная t, которая является независимой переменной в NDSolve. Нужно записать все коэффициенты в виде функций от t, или присваивать им значения лениво (через :=), сейчас проблема в том, что Математика пытается численно посчитать интегралы, в которые t входит просто как буковка, без привязанного к ней значения.

[Zepp]

Хорошо, заменим NIntegrate -> Integrate. Теперь ошибка в вычислении производной в начальной точке.

$NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`. >>$

upd1: ни перенос н.у. на другой момент времени, ни добавление параметра SolveDelayed -> True по следам хелпа результата не дало.

[FFFF]

Посмотрите в хелпе, как правильно записывать функции в Mathematica (например, вместо sqrt(x) написать Sqrt[x]), также нужно избавиться от всех коэффициентов в выражении, которое считает NDSolve, я, например, не могу понять, где определяется коэффициент b11?

[Zepp]

b11 остался случайно, забыл потереть в ходе упрощения задачи)

k1 = 6000;
k2 = 60000;
Ix = 0.2;
Iz = 0.05;
Ixz = 0.007;
Xc = 0.025;
Zc = 0.1;
ro = 0.7;
l = 0.3;

With[{a11 = k1, a22 = k2, c11 = Ix, c22 = Iz, c12 = -Ixz, c21 = -Ixz,
  V = 350 + t,
  M = V/330,
  k = (0.25 + 0.07*M*HeavisideTheta[3 - M] + 
     0.21*HeavisideTheta[M - 3]),
  K1 = 2*ro/(Sqrt [M^2 - 1]),
  d011 = (0.25 - 0.5*z)*z^2,
  d012 = ((0.15 - 0.5*z)/(0.25 - 0.5*z) + ((k*z)/2 - z/2 - k/4 + 
         3/20)/(0.25 - 0.5*z))*(0.25 - 0.5*z)^2*z, 
  d022 = (((0.15 - 0.5*z)/(0.25 - 0.5*z) + ((k*z)/2 - z/2 - k/4 + 
            3/20)/(0.25 - 0.5*z))^2 + 1/12)*(0.25 - 0.5*z)^3,
  b022 = ((0.15 - 0.5*z)/(0.25 - 0.5*z) + ((k*z)/2 - z/2 - k/4 + 
         3/20)/(0.25 - 0.5*z))*(0.25 - 0.5*z)^2,
  b012 = (0.25 - 0.5*z)*z,
  d11 = K1*Integrate[d011, {z, 0, l}],
  d12 = K1*Integrate[d012, {z, 0, l}],
  d21 = d12,
  d22 = K1*Integrate[d022, {z, 0, l}],
  b12 = -K1*Integrate[b012, {z, 0, l}],
  b22 = -K1*Integrate[b022, {z, 0, l}]
  }, sol = 
  NDSolve[{c11*q''[t] + d11*V*q'[t] + a11*q[t] + c12*r''[t] + 
      d12*V*r'[t] + b12*V^2*r[t] == 10 sin (t), 
    c21*q''[t] + d21*V*q'[t] + a22*r[t] + c22*r''[t] + d22*V*r'[t] + 
      b22*V^2*r[t] == 10 sin (t), q[1] == 1, q'[1] == 2, r[1] == 3, 
    r'[1] == 4}, {q, r}, {t, 1, 100}, SolveDelayed -> True]]

Plot[Evaluate[{q[t], r[t]}], {t, 1, 100}]

NDSolve::nlnum: The function value {0.,6000. -110881. b012 K1+210.6 d011 K1+421.2 d012 K1-10. sin,0.,180000. +702. d12-110881. b022 K1+421.2 d022 K1-10. sin}
 is not a list of numbers with dimensions {4} at {t,q[t],(q^\[Prime])[t],r[t],(r^\[Prime])[t],(q^\[Prime])[t],(q^\[Prime]\[Prime])[t],(r^\[Prime])[t],(r^\[Prime]\[Prime])[t]} = {1.,1.,2.,3.,4.,2.,0.,4.,0.}. >>

Можно ли задать коэффициенты, вычисляющиеся интегралами и зависящие от V(t), неким другим способом?

[FFFF]

Во-первых, синус так и не исправлен, во-вторых, проблема может быть в том, что блоке определения переменных одни переменные определяются через другие (например, сначала определяется M, а потом в дальнейшие коэффициенты входит M, это может не работать), в-третьих, можно попробовать заменить равенства (=) на ленивое присваивание (:=), или явно определить коэффициенты как функции от t (например, d11[x_]:=K1*Integrate[d011, {z, 0, l}]/. t->x)

[Zepp]

Упростил донельзя - теперь подинтегральные выражения вообще не зависят от t, но все равно не выходит :/

k1 = 6000,
k2 = 60000,
Ix = 0.2,
Iz = 0.05,
Ixz = 0.007,
Xc = 0.025,
Zc = 0.1,
ro = 0.7,
l = 0.3,
k = 0.33,
With[{a11 = k1, a22 = k2, c11 = Ix, c22 = Iz, c12 = -Ixz, c21 = -Ixz, 
  ro = 0.7,
  V := 16.5*t + 396,
  M := 1.2 + 0.05*t,
  
  K1 := 2*ro/(Sqrt [(1.2 + 0.05*t)^2 - 1]),
  d011 := (0.25 - 0.5*z)*z^2,
  d012 := ((0.15 - 0.5*z)/(0.25 - 0.5*z) + ((0.33*z)/2 - z/2 - 
         0.33/4 + 3/20)/(0.25 - 0.5*z))*(0.25 - 0.5*z)^2*z,
  d022 := (((0.15 - 0.5*z)/(0.25 - 0.5*z) + ((0.33*z)/2 - z/2 - 
            0.33/4 + 3/20)/(0.25 - 0.5*z))^2 + 1/12)*(0.25 - 0.5*z)^3,
  b022 := ((0.15 - 0.5*z)/(0.25 - 0.5*z) + ((0.33*z)/2 - z/2 - 
         0.33/4 + 3/20)/(0.25 - 0.5*z))*(0.25 - 0.5*z)^2,
  b012 := (0.25 - 0.5*z)*z,
  d11 = NIntegrate[d011, {z, 0, l}],
  d12 = NIntegrate[d012, {z, 0, l}],
  d21 = NIntegrate[d012, {z, 0, l}],
  d22 = NIntegrate[d022, {z, 0, l}],
  b12 = NIntegrate[b012, {z, 0, l}],
  b22 = NIntegrate[b022, {z, 0, l}]
  }, sol = 
  NDSolve[{c11*q''[t] + K1*d11*V*q'[t] + a11*q[t] + c12*r''[t] + 
      K1*d12*V*r'[t] - K1*b12*V^2*r[t] == 10 Sin [t], 
    c21*q''[t] + K1*d21*V*q'[t] + a22*r[t] + c22*r''[t] + 
      K1*d22*V*r'[t] - K1*b22*V^2*r[t] == 10 Sin [t],
    q[1] == 1, q'[1] == 2, r[1] == 3, r'[1] == 4}, {q, r}, {t, 1, 
    100}, SolveDelayed -> True]]

Plot[Evaluate[{q[t], r[t]}], {t, 1, 100}]

NDSolve::nlnum: The function value {0.,24000. +492804. b11+628325. b012 K1+631.8 d011 K1-1053. d012 K1-10. sin,0.,-1.02*10^6+492804. b21+2106. d12+628325. b022 K1-1053. d022 K1-10. sin}
 is not a list of numbers with dimensions {4} at {t,q[t],(q^\[Prime])[t],r[t],(r^\[Prime])[t],(q^\[Prime])[t],(q^\[Prime]\[Prime])[t],(r^\[Prime])[t],(r^\[Prime]\[Prime])[t]} = {1.,4.,6.,-17.,-10.,6.,0.,-10.,0.}. >>

ЧЯДНТ?

[FFFF]

В определении К1 участвует ro. Попробуйте вообще убрать With.

[Zepp]

k1 = 6000,
k2 = 60000,
Ix = 0.2,
Iz = 0.05,
Ixz = 0.007,

k = 0.33,
a11 = k1,
a22 = k2,
c11 = Ix,
c22 = Iz,
c12 = -Ixz,
c21 = -Ixz,
V = 16.5*t + 396,

K1 := 2*0.7/(Sqrt [(1.2 + 0.05*t)^2 - 1]),

d11 = 0.001,
d12 = 0.07,
d21 = 0.43,
d22 = 12,
b12 = 4,
b22 = 3,
sol = NDSolve[{
   c11*q''[t] + K1*d11*V*q'[t] + a11*q[t] + c12*r''[t] + 
     K1*d12*V*r'[t] - K1*b12*V^2*r[t] == 10 Sin [t],
   c21*q''[t] + K1*d21*V*q'[t] + a22*r[t] + c22*r''[t] + 
     K1*d22*V*r'[t] - K1*b22*V^2*r[t] == 10 Sin [t],
   q[1] == 1, q'[1] == 2, r[1] == 3, r'[1] == 4}, {q, r}, {t, 1, 100},
   SolveDelayed -> True]

Plot[Evaluate[{q[t], r[t]}], {t, 1, 100}]

Syntax::tsntxi: "k1=6000,k2=60000,Ix=0.2,Iz=0.05,Ixz=0.007,k=0.33,a11=k1,a22=k2,c11=Ix,c22=Iz,c12=-Ixz,c21=-Ixz,V=16.5*t+396,K1:=2*0.7/(Sqrt[(1.2+0.05*t)^2-1]),d11:=0.001,d12:=0.07,d21:=0.43,d22:=12,b12:=4,b22:=3,sol=NDSolve[<<1>>]" is incomplete; more input is needed.

Чего не хватает?

[Vince Diesel]

Хелп слабо почитать? Там есть вводный курс, где сказано, к примеру, что в конце ставится точка с запятой, а не запятая. Ну и в сообщении NDSolve говорит, что ему передан не численный вектор. Вот всем переменным (буквам), которые там есть, надо присвоить численные значения.