Asked by danish rafiq
on 23 Apr 2019

I am trying to run the following function via ode45 solver but not able to model the equations for it.

Problem is like this:

I have model in state-space representation form:

where x represents the state vector matrix .

Dimension of is , is , is wher are specified.

State equation are of the form:

and the initaila conditions are:

x1=2.3*ones(N,1);

x2=zeros(N,1);

x3=0.1013*ones(M*N,1);

How do I incorporate the state variables as vectors in the function?

Here is the code:

N=40;M=20;

n=(M+2)N;

utest=@(t) cos(4*pi*t);

tsim=[0:0.001:10];

x1=2.3*ones(N,1);

x2=zeros(N,1);

x3=0.1013*ones(M*N,1);

x0=[x1;x2;x3];

[consts, operators] = MemsConstsAndOperators(N,M); %function provided below

[tout,xout]=ode45(@(t,x)testfunc(t,x,consts,operators),[0 1],x0); % additional arguments passed to testfunc: (consts,operators)

B=zeros(size(x0,1),1);

B(N+1:2*N,:)=1;

C=zeros(1,M*N+2*N);

C(1,N/2)=1;

y=C*xout';

plot(y)

function xdot=testfunc(t,x,consts,operators)

M=consts(15);

N=consts(16);

S0 = consts(1);

E = consts(2);

In0 = consts(3);

rho0 = consts(4);

lambda = consts(5);

p0 = consts(6);

z0 = consts(7);

width=consts(8);

viscosity = consts(10);

height = consts(12);

del = consts(14);

B = operators{5};

rho = rho0*height*width;

In = In0*width*height^3;

INT = operators{1};

D4 = operators{10};

D2 = operators{3};

grady = operators{4};

gradx = operators{8};

LAP = operators{9};

gradx2 = operators{7};

C= operators(6);

%Here is the nonlinear function:

arg31 = (x(2)/(3*x(1)^3))*1/del;

f31 = -DiagHat((arg31),M,N) * x(3); %this is the line where error shows up "To RESHAPE the number of elements must not change".

arg32 = (x(1)+4*lambda)/(4*viscosity).*(gradx2*(x(1)-z0));

f32 = DiagHat((arg32),M,N) * (x(3)*(gradx*(x(3)-p0)));

arg33 = ((x(1)^2+6*lambda*x(1))/(12*viscosity));

f33 = DiagHat((arg33),M,N) * (x(3).*(LAP*(x(3)-p0))+(gradx*(x(3)-p0)).^2+(grady*(x(3)-p0))^2);

f3=f31+f32+f33;

%Here is the nonlinear function

xdot=[x(2)/(3*x(1)^2)*1/del;

(2*x(2)^2)/(3*x(1)^3)*1/del + (3*x(1)^2)/(rho)*(INT*(x(3)-p0)+(S0*height*width)*(D2*(x(1)-z0))-E*In*(D4*(x(1)-z0)))*del;

f3];

end

function Jtemp = DiagHat(x,M,N)

Jtemp = diag(reshape(repmat(x',M,1),N*M,1));

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%function that calculates the constants for the main function

function [consts, operators, xold] = MemsConstsAndOperators(N,M)

del = 1e-7;

height = 2.2;

long = 610;

width = 40;

z0=2.3;

dx=long/(N+1);

dy=width/(M+1);

% Material and air parameters

epsilon0=8.854e-6;

rho0=2330e-18;

E=149.0e3;

S0=-3.7;

In0=1/12;

atm_Pa_conversion=101.3e-3;

viscosity=1.82e-11;

lambda=0.064;

p0=1*atm_Pa_conversion;

% Input parameter

b0 = -(3*epsilon0*del/(2*rho0*height));

consts = [S0 E In0 rho0 lambda p0 z0 width dy viscosity dx height b0 del M N];

xold = zeros(N*(M+2),1); xold(1:N) = z0; xold((2*N+1):((M+2)*N)) = p0;

IntSeg = ones(1,M);

IntSeg = dy*IntSeg;

IntRow = zeros(1,M*N); IntRow(1:M) = IntSeg;

for i = 1:N

INT(i,:) = IntRow;

IntRow = circshift(IntRow,[0 M]);

end

INT=sparse(INT);

% Operator D4

D4 = zeros(N,N);

D4 = 6*diag(ones(N,1))-4*diag(ones(N-1,1),1)-4*diag(ones(N-1,1),-1)+diag(ones(N-2,1), 2)+diag(ones(N-2,1),-2);

D4(1,1)=D4(1,1)+1;

D4(N,N)=D4(N,N)+1;

D4=(1.0/dx^4)*D4;

D4 = sparse(D4);

% Operator D2

D2 = zeros(N,N);

D2 = -2*diag(ones(N,1))+diag(ones(N-1,1),1)+diag(ones(N-1,1),-1);

D2(1,1)=-1;

D2(N,N)=-1;

D2=(1.0/dx^2)*D2;

D2 = sparse(D2);

% grady

grady = zeros(M*N,M*N);

grady(1:M,1:M) = (1.0/(2.0*dy))*(diag(ones(M-1,1),1)-diag(ones(M-1,1), -1));

for index=2:N

grady([((index-1)*M+1):(index*M)], [((index-1)*M+1):(index*M)]) = ...

grady([1:M],[1:M]);

end;

grady=sparse(grady);

% gradx

gradx = zeros(M*N,M*N);

e = ones((N+1)*M,1)/(2.0*dx);

gradx = spdiags([e -1.0*e],[M -M],(M*N),(M*N));

for index = 1:M

gradx(index,index) = -1.0/(2*dx);

gradx(M*(N-1)+index, M*(N-1)+index) = 1.0/(2*dx);

end;

gradx = sparse(gradx);

e=1/dx^2; f1=-2*1/dx^2; f2=-2*1/dy^2; f=f1+f2; g=1/dy^2;

diag0=ones(N*M,1);

diag1=ones(N*M-1,1);

diag0(1:M)=1/2;

diag0((N-1)*M+1:N*M)=1/2;

for i=1:N-1

diag1(i*M)=0;

end

LAP= zeros(M*N, M*N);

LAP=f1*diag(diag0)+f2*diag(ones(N*M,1))+g*diag(diag1,1)+g*diag(diag1,-1)+e*diag(ones((N-1)*M,1),M)+e*diag(ones((N-1)*M,1),-M);

LAP=sparse(LAP);

% gradx2

gradx2 = zeros(N,N);

gradx2 = (1/(2*dx))*(diag(ones(N-1,1),1)-diag(ones(N-1,1),-1));

gradx2=sparse(gradx2);

% Operator B -- input vector

B = zeros((M+2)*N,1);

B((N+1):(2*N)) = 1;

B = B*b0;

B=sparse(B);

% Operator C -- Output vector

C1=zeros(1,(M+2)*N); C2 = C1; C3 = C1;

C1(round(N/2))=1; % Vertical displacement at center

C3(2*N+round(M*N/2)) = 1; % Pressure under beam center

C2(N+round(N/2)) = 1; % du^3/dt at center of beam

C = [C1];

C=sparse(C);

E = speye(length(B));

operators = [{INT} {E} {D2} {grady} {B} {C} {gradx2} {gradx} {LAP} {D4}];

Opportunities for recent engineering grads.

Apply Today
## 4 Comments

## darova (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/457927-how-to-solve-nonlinear-vectorized-state-space-equations-using-ode-solvers#comment_697200

## danish rafiq (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/457927-how-to-solve-nonlinear-vectorized-state-space-equations-using-ode-solvers#comment_697302

## darova (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/457927-how-to-solve-nonlinear-vectorized-state-space-equations-using-ode-solvers#comment_698708

## danish rafiq (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/457927-how-to-solve-nonlinear-vectorized-state-space-equations-using-ode-solvers#comment_698846

Sign in to comment.