MATLAB Answers

0

How to solve nonlinear vectorized state space equations using ode solvers?

Asked by danish rafiq on 23 Apr 2019
Latest activity Commented on by danish rafiq on 27 Apr 2019 at 5:42
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}];

  4 Comments

Show 1 older comment
This function takes the imput vector x and makes a diagonal matrix of the form:
test.jpg
What is the size of arg31 variable?
arg31 = (x(2)/(3*x(1)^3))*1/del;
f31 = -DiagHat((arg31),M,N) * x(3);
Your function:
function Jtemp = DiagHat(x,M,N)
% x - you assume size(x) == (N,1) (column vector)
x1 = repmat(x',M,1); % size(x1) == (M,N) (matrix)
x2 = reshape(x1,N*M,1); % size(x2) == (N*M,1) (column vector)
Jtemp = diag(x2); % matrix (N*M, N*M)
end
What will happen if arg31 is (1,1) size?
x1 = repmat(x',M,1); % size(x1) == (M,1) (column vector)
% and here you are trying to make vector of N*M long
% from vector of long M
x2 = reshape(x1,N*M,1);
Yes that's the issue here.
Basically the code is written for dicretized state variables of dimensions N,N, and N*M respectively. and when try to run the function f=@(x) [f1;f2;f3] using an ode solver, it shows a matrix of dimension with each row as the first row containg the initial values only. There is no time stepping in the variables.

Sign in to comment.

0 Answers