import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

x1=0
x2=0
x3=0

def dxdt(X, t):
    x1, x2, x3 = X
    dx1dt = x2
    dx2dt = x3
    dx3dt = 12*x1 - 6*x2 + 9*x3 + 3
    return dx1dt, dx2dt, dx3dt

t = np.linspace(0, 1, 1000000)
X0 = x1, x2, x3
sol = odeint(dxdt, X0, t)
print (sol)


plt.plot(t, sol[:, 0], 'r', label='x1(t)')
plt.plot(t, sol[:, 1], 'g', label='x2(t)')
plt.plot(t, sol[:, 2], 'b', label='x3(t)')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()

Maar als ik het gewooon oplos in Python met odeint, dan kom ik dezelfde oplossing uit als in mijn statemodel.
Dat doet mij vermoeden dat er met het state spacemodel niets mis is.
Kan het zijn dat de fout eerder zit in jullie alternatieve oplossing?


[code]import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

x1=0
x2=0
x3=0

def dxdt(X, t):
x1, x2, x3 = X
dx1dt = x2
dx2dt = x3
dx3dt = 12*x1 - 6*x2 + 9*x3 + 3
return dx1dt, dx2dt, dx3dt

t = np.linspace(0, 1, 1000000)
X0 = x1, x2, x3
sol = odeint(dxdt, X0, t)
print (sol)


plt.plot(t, sol[:, 0], 'r', label='x1(t)')
plt.plot(t, sol[:, 1], 'g', label='x2(t)')
plt.plot(t, sol[:, 2], 'b', label='x3(t)')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()
[/code]



[attachment=0]Figure_1.png[/attachment]

syms x1(t) x2(t) x3(t)

ode1 = diff(x1) == x2;
ode2 = diff(x2) == x3;
ode3 = diff(x3) == 12*x1 - 6*x2 + 9*x3 + 3;
odes = [ode1; ode2;, ode3]

S=dsolve(odes)
simplify(S.x1)
simplify(S.x2)
simplify(S.x3)

Probleem is dat op het niet-homogene systeem matlab blijft hangen. Ik heb even gewerkt met 3 ipv met een stapfunctie om het iets eenvoudiger te maken.

[code]syms x1(t) x2(t) x3(t)

ode1 = diff(x1) == x2;
ode2 = diff(x2) == x3;
ode3 = diff(x3) == 12*x1 - 6*x2 + 9*x3 + 3;
odes = [ode1; ode2;, ode3]

S=dsolve(odes)
simplify(S.x1)
simplify(S.x2)
simplify(S.x3)[/code]

Als je dit met de hand wil doen, dan verloopt dat via de eigenwaarden en eigenvectoren van A, die je kan herkennen in de oplossing.

[attachment=0]dif1.png[/attachment]

De positieve reële delen van de eigenwaarden wijzen erop dat het systeem instabiel is.
Om de niet homogene DV op te lossen kan je werken met bvb variatie van parameters. Dat principe is uitbreidbaar naar stelsels van DV's. Die integralen gaan voor een stapfunctie u in het rechterlid volgens mij analytisch oplosbaar zijn.

syms x1(t) x2(t) x3(t)

ode1 = diff(x1) == x2;
ode2 = diff(x2) == x3;
ode3 = diff(x3) == 12*x1 - 6*x2 + 9*x3;
odes = [ode1; ode2;, ode3]

S=dsolve(odes)
simplify(S.x1)
simplify(S.x2)
simplify(S.x3)

Eventjes symbolisch met matlab het homogene systeem opgelost.

[code]
syms x1(t) x2(t) x3(t)

ode1 = diff(x1) == x2;
ode2 = diff(x2) == x3;
ode3 = diff(x3) == 12*x1 - 6*x2 + 9*x3;
odes = [ode1; ode2;, ode3]

S=dsolve(odes)
simplify(S.x1)
simplify(S.x2)
simplify(S.x3)
[/code]

Het zijn wel hele lange antwoorden.


x_1 =
C2*(sin((3^(1/2)*t*((233^(1/2) + 24)^(1/3) - (24 - 233^(1/2))^(1/3)))/2)*exp(-t*((233^(1/2) + 24)^(1/3)/2 + (24 - 233^(1/2))^(1/3)/2 - 3))*((699^(1/2)*(24 - 233^(1/2))^(2/3))/2352 - (25*3^(1/2)*(24 - 233^(1/2))^(2/3))/2352 + ((17*3^(1/2) + 699^(1/2))*(24 - 233^(1/2))^(1/3))/336) - exp(-t*((233^(1/2) + 24)^(1/3)/2 + (24 - 233^(1/2))^(1/3)/2 - 3))*cos((3^(1/2)*t*((233^(1/2) + 24)^(1/3) - (24 - 233^(1/2))^(1/3)))/2)*((233^(1/2)*(24 - 233^(1/2))^(1/3))/336 - (233^(1/2)*(24 - 233^(1/2))^(2/3))/2352 + (17*(24 - 233^(1/2))^(1/3))/336 + (25*(24 - 233^(1/2))^(2/3))/2352 + 5/12)) - C3*(exp(-t*((233^(1/2) + 24)^(1/3)/2 + (24 - 233^(1/2))^(1/3)/2 - 3))*cos((3^(1/2)*t*((233^(1/2) + 24)^(1/3) - (24 - 233^(1/2))^(1/3)))/2)*((699^(1/2)*(24 - 233^(1/2))^(2/3))/2352 - (25*3^(1/2)*(24 - 233^(1/2))^(2/3))/2352 + ((17*3^(1/2) + 699^(1/2))*(24 - 233^(1/2))^(1/3))/336) + sin((3^(1/2)*t*((233^(1/2) + 24)^(1/3) - (24 - 233^(1/2))^(1/3)))/2)*exp(-t*((233^(1/2) + 24)^(1/3)/2 + (24 - 233^(1/2))^(1/3)/2 - 3))*((233^(1/2)*(24 - 233^(1/2))^(1/3))/336 - (233^(1/2)*(24 - 233^(1/2))^(2/3))/2352 + (17*(24 - 233^(1/2))^(1/3))/336 + (25*(24 - 233^(1/2))^(2/3))/2352 + 5/12)) + C4*exp(t*((233^(1/2) + 24)^(1/3) + (24 - 233^(1/2))^(1/3) + 3))*((233^(1/2)*(24 - 233^(1/2))^(1/3))/168 - (233^(1/2)*(24 - 233^(1/2))^(2/3))/1176 + (17*(24 - 233^(1/2))^(1/3))/168 + (25*(24 - 233^(1/2))^(2/3))/1176 - 5/12)


x_2=
C2*(sin((3^(1/2)*t*((233^(1/2) + 24)^(1/3) - (24 - 233^(1/2))^(1/3)))/2)*exp(-t*((233^(1/2) + 24)^(1/3)/2 + (24 - 233^(1/2))^(1/3)/2 - 3))*((23*3^(1/2)*(24 - 233^(1/2))^(2/3))/1176 + (699^(1/2)*(24 - 233^(1/2))^(2/3))/392 + ((3*3^(1/2) + 699^(1/2))*(24 - 233^(1/2))^(1/3))/168) + exp(-t*((233^(1/2) + 24)^(1/3)/2 + (24 - 233^(1/2))^(1/3)/2 - 3))*cos((3^(1/2)*t*((233^(1/2) + 24)^(1/3) - (24 - 233^(1/2))^(1/3)))/2)*((233^(1/2)*(24 - 233^(1/2))^(2/3))/392 - (233^(1/2)*(24 - 233^(1/2))^(1/3))/168 - (24 - 233^(1/2))^(1/3)/56 + (23*(24 - 233^(1/2))^(2/3))/1176 + 1/6)) - C3*(exp(-t*((233^(1/2) + 24)^(1/3)/2 + (24 - 233^(1/2))^(1/3)/2 - 3))*cos((3^(1/2)*t*((233^(1/2) + 24)^(1/3) - (24 - 233^(1/2))^(1/3)))/2)*((23*3^(1/2)*(24 - 233^(1/2))^(2/3))/1176 + (699^(1/2)*(24 - 233^(1/2))^(2/3))/392 + ((3*3^(1/2) + 699^(1/2))*(24 - 233^(1/2))^(1/3))/168) - sin((3^(1/2)*t*((233^(1/2) + 24)^(1/3) - (24 - 233^(1/2))^(1/3)))/2)*exp(-t*((233^(1/2) + 24)^(1/3)/2 + (24 - 233^(1/2))^(1/3)/2 - 3))*((233^(1/2)*(24 - 233^(1/2))^(2/3))/392 - (233^(1/2)*(24 - 233^(1/2))^(1/3))/168 - (24 - 233^(1/2))^(1/3)/56 + (23*(24 - 233^(1/2))^(2/3))/1176 + 1/6)) + C4*exp(t*((233^(1/2) + 24)^(1/3) + (24 - 233^(1/2))^(1/3) + 3))*((233^(1/2)*(24 - 233^(1/2))^(1/3))/84 - (233^(1/2)*(24 - 233^(1/2))^(2/3))/196 + (24 - 233^(1/2))^(1/3)/28 - (23*(24 - 233^(1/2))^(2/3))/588 + 1/6)


x_3=
exp(-t*((233^(1/2) + 24)^(1/3)/2 + (24 - 233^(1/2))^(1/3)/2 - 3))*(C2*cos((3^(1/2)*t*((233^(1/2) + 24)^(1/3) - (24 - 233^(1/2))^(1/3)))/2) + C4*exp(t*((3*(233^(1/2) + 24)^(1/3))/2 + (3*(24 - 233^(1/2))^(1/3))/2)) + C3*sin((3^(1/2)*t*((233^(1/2) + 24)^(1/3) - (24 - 233^(1/2))^(1/3)))/2))

Instabiel systeem betekent positieve eigenwaarden, maar dat verandert echt niks aan het feit dat beide oplossingen tot hetzelfde zouden moeten leiden. Ik zal later verder zoeken. De homogene is ook gemakkelijk analytisch met de hand op te lossen. Ik zal dat later hier posten.

[attachment=1]Mapleplot.png[/attachment][attachment=0]rootlocusplot.png[/attachment]
Het is een instabiel systeem. zou dat een verklaring voor het verschil kunnen zijn?

Als ik het juist begrijp, dan correspondeert mijn oplossing met de state space oplossing van Ukster in Maple, maar niet met de oplossing van Xilvo en Ukster als ze gewoon de DV oplossen.

A = [0, 1, 0; 0,0, 1; 12, -6, 9];
B = [0;0;3];
C = [0,1,0];
D = [0];
x0 = [0 ; 0; 0];
sys = ss(A,B,C,D)
initial(sys,x0)

t = 0:0.00001:1;  
u = heaviside(t);

lsimplot(sys,u,t)
grid on

Dit krijg ik met het Matlab State Space model.

[code]A = [0, 1, 0; 0,0, 1; 12, -6, 9];
B = [0;0;3];
C = [0,1,0];
D = [0];
x0 = [0 ; 0; 0];
sys = ss(A,B,C,D)
initial(sys,x0)

t = 0:0.00001:1;
u = heaviside(t);

lsimplot(sys,u,t)
grid on[/code]


[attachment=0]statespace.jpg[/attachment]

Nee,[attachment=1]1.png[/attachment][attachment=0]2.png[/attachment]

Voor de eerste afgeleide krijg ik 2409,7, voor x 208,5
Heb je niet toevallig de verkeerde variable geplot?
(ik heb me verder nog niet in het model verdiept).

met de toestandsvariabelen:
[attachment=0]toestandsvariabelen.png[/attachment]

ik blijf het vreemd vinden.[attachment=1]netwerk.png[/attachment]
ik gebruik voor dit netwerk geen initial conditions (volgens mij default nul)

de responsplot[attachment=2]stapresponsie.png[/attachment] vanuit het state space model[attachment=0]statespacemodel.png[/attachment] en is exact gelijk aan die van het professionele simulatiesoftwarepakket [b]Microcap[/b] (stapresponsie 1V)

Ik zal eens proberen in matlab.

ingevoerd: ic := [x1(0) = 0, x2(0) = 0, x3(0) = 0]
Helaas, nog steeds dezelfde melding!

[quote=ukster post_id=1177159 time=1686485740 user_id=76233]
ic := [y1(0) = 0, D(y1)(0) = 0]
opgenomen in de opdracht: ResponsePlot(sys, Step(1), output = y1(t), initialconditions = ic, duration = 1)
geeft [color=#FF0040]Error, (in DynamicSystems:-Simulate) initial conditions for a state-space model must be a Vector.[/color]

Nu nog uitzoeken hoe dat moet en of ic wel compleet is, (de toestandsvariabelen zitten er immers nog niet in)
bedoel je dit?
ic := [y1(0) = 0, D(y1)(0) = 0, x1(0) = 0, x2(0) = 0, x3(0) = 0]
[/quote]

Ik verwacht iets genre

ic := [x1(0) = 0, x2(0) = 0, x3(0) = 0]

De toestandsvariabelen specifiëren is voldoende. Probeer dan nog eens opnieuw...