How to solve nonlinear problem of magnetic bearing shaft system ?

function dx=ODE_magnetic(t,x)

dx=zeros(4,1); x=zeros(4,1);

m=5.65;

g=9.81;

Aa=1.44e-4;

N=80;

co=0.0002;

Io=4;

muo=pi*4e-7;

io=m*g*co^2/(muo*Aa*N^2*Io);

k=4.05e-4;

c=1;

w=10000;

e=0.001;

kxx=-(muo*Aa*N^2*Io^2)/co^3

kxi=(muo*Aa*N^2*Io)/co^2

kyi=(muo*Aa*N^2*Io)/co^2

kyy=-(muo*Aa*N^2*(Io^2+io^2))/co^3

kyip=-(2*muo*Aa*N^2*io)/co^3

kyyp=(3*muo*Aa*N^2*Io*io)/co^4

Kdx=1/kxi

Kpx=(k-kxx)/kxi

Kdy=1/kyi

Kpy=(k-kyy)/kyi

kx1=(kxi*Kdx*sqrt(c*g))/(m*g)

kx2=c*(kxi*Kpx+kxx)/(m*g)

ky1=(kyi*Kdy*sqrt(c*g))/(m*g)

ky2=(c*kyip*Kdy*sqrt(c*g))/(m*g)

ky3=(c*(kyi*Kpy+kyy))/(m*g)

ky4=(c^2*(kyyp+kyip*Kpy))/(m*g)

dx(1)=x(3)

dx(2)=x(4)

dx(3)=-(kx1/w)*x(3)-(kx2/w^2)*x(1)+e*cos(w*t)

dx(4)=-(ky1/w)*x(4)-(ky2/w)*x(4)*x(2)-(ky3/w^2)*x(2)-(ky4/w^2)*x(2).^2+e*sin(w*t)

[John D'Errico]

What does your question have to do with the mess of undocumented code that you show?

Answer 1

[Are Mjaavatten]

Your code will not work in an ODE solver because you reset the input value x on line 2. Get rid of that and make sure all lines end with a semicolon to avoid a hopeless amount of unneeded output. Then you may try solving your ODE system:

x0 = [1;1;1;1];  % or whatever
[t,x] = ode45(@ODE_magnetic,[0,5],x0);
plot(t,x)