MATLAB Answers

calculation of inverse matrix

3 views (last 30 days)
raha ahmadi
raha ahmadi on 5 Apr 2021
Edited: raha ahmadi on 5 Apr 2021
I need to calculate the inverse matrix of P2 but I dont know my matrix is sparse or not. also I knew that this problem was solved (I saw the results in a book)
but when I raise 'J', (2^(J+1) is the dimention of my space) I get the singular matrix error even for J=3 however the book has solved for J=6, I dont know exactly in this situation It whould be better to use matlab built-in function or there is more effective algorithm for invers matrix calculation.
This code indeed solve a PDE by wavelet method .
Thanks in advance
% non homogeneous Haar wavelet
close all
clear all
clc
%% wavelet Twaveoluton
co=1/24;
J=4;
M=pow2(J);
M2=2*M;
A=0;
B=1;
dx=(B-A)/M2;
x(M)=A+M*dx;
ksi1(1)=A; ksi2(1)=B; ksi3(1)=B;
ksi1(2)=A; ksi2(2)=x(M); ksi3(2)=B;
H=zeros(4);
P1=H;
P2=H;
P3=H;
P4=H;
for j=1:J
m=pow2(j);
mu=round(M/m);
x(mu)=A+mu*dx;
x(2*mu)=A+(2*mu)*dx;
ksi1(m+1)=A;
ksi2(m+1)=x(mu);
ksi3(m+1)=x(2*mu);
for k=1:m-1
i=m+k+1;
x(2*k*mu)=A+2*k*mu*dx;
x((2*k+1)*mu)=A+(2*k+1)*mu*dx;
x(2*(k+1)*mu)=A+(2*(k+1)*mu)*dx;
ksi1(i)=x(2*k*mu);
ksi2(i)=x((2*k+1)*mu);
ksi3(i)=x(2*(k+1)*mu);
end
end
xc(1)=0.5*(x(1)+A);
for l=2:M2
xc(l)=0.5*(x(l)+x(l-1));
c(l)=(ksi2(l)-ksi1(l))/(ksi3(l)-ksi2(l));
end
for i=1:M2
K(i)=0.5*(ksi2(i)-ksi1(i))*(ksi3(i)-ksi1(i));
for l=1:M2
X=xc(l);
if X<ksi1(i)
H(i,l)=0;
P1(i,l)=0; P2(i,l)=0;P3(i,l)=0;P4(i,l)=0;P6(i,l)=0;
elseif X<ksi2(i)
H(i,l)=1;
P1(i,l)=X-ksi1(i);
P2(i,l)=0.5*(X-ksi1(i)).^2;
P3(i,l)=(X-ksi1(i)).^3/6;
P4(i,l)=co*(X-ksi1(i)).^4;
P6(i,l)=1/(factorial(6))*(X-ksi1(i))^6;
elseif X<ksi3(i)
H(i,l)=-c(i);
P1(i,l)=c(i)*(ksi3(i)-X);
P2(i,l)=K(i)-0.5*c(i)*(ksi3(i)-X).^2;
P3(i,l)=K(i)*(X-ksi2(i))+(ksi3(i)-X).^3/6;
P4(i,l)=co*((X-ksi1(i)).^4-2*(X-ksi2(i)).^4);
P6(i,1)=1/(factorial(6))*((X-ksi1(i))^6-(1+c(i))*(X-ksi2(i))^6);
elseif X>=ksi3(i)
H(i,l)=0;
P1(i,l)=0; P2(i,l)=K(i);
P3(i,l)=K(i)*(X-ksi2(i));
P4(i,l)=co*((X-ksi1(i)).^4-2*(X-ksi2(i)).^4+(X-ksi3(i)).^4);
P6(i,1)=1/(factorial(6))*((X-ksi1(i))^6-(1+c(i))*(X-ksi2(i))^6+c(i)*(X-ksi3(i))^6);
else
end
F(l)=X-1-0.15*X.^2-X.^5./factorial(5)+sin(3/2.*X).*sin(X./2);
end
end
%% calculation the solution of the sine Gordon PDE by wavelet expanstion : 1/L^2V"-V''=sin v
L=10;
dt=1e-3;
tfin=30;
tin=0;
S=(tfin-tin)/dt;
%t=zeros(1,S);
L=20;
beta=0.025;
alpha=L/(sqrt(1-L^2*beta^2));
A=zeros(S+1,2*M);
Vzegx=A;
V=A;
Vdotzegx=A;
Vdot=A;
t=zeros(S+1,1);
B=A;
for d=1:S+2
t(d)=tin+(d-1)*dt;
end
si=(-4*alpha*beta*exp(-alpha*beta.*t))./(1+exp(-2*alpha*beta.*t));
sidot =(4*alpha^2*beta^2*exp(-alpha*beta.*t).*((1-exp(-2*alpha*beta.*t))))./...
((1+exp(-2*alpha*beta.*t)).^2);
si2dot=(4*alpha^3*beta^3*exp(-alpha*beta.*t).*(-1+5*exp(-4*alpha*beta.*t+...
5*exp(-2*alpha*beta.*t)-exp(-6*alpha*beta.*t))))./((1+exp(-2*alpha*beta.*t)).^4);
fi=4*atan(exp(alpha.*t));
fidot=(-4*alpha*beta*exp(-alpha*beta.*t))./(1+exp(-2*alpha*beta.*t));
fi2dot=(4*alpha^2*beta^2*exp(-alpha*beta.*t).*(1+exp(-2*alpha*beta.*t)))...
./((1+exp(-2*alpha*beta.*t)).^2);
for c=1:S+1
%for f=1:M2
P2inv=eye(2*M)/P2;
if c==1
V(1,:)=4*atan(exp(alpha.*xc));
Vzegx(1,:)=(4*alpha^2.*exp(alpha.*xc).*(1-exp(2*alpha.*xc)))/...
(1+exp(2*alpha.*xc)).^2;
B(1,:)=(1/L^2.*Vzegx(1,:)-sin(V(1,:))-fi2dot(2).*ones(1,2*M)-...
si2dot(2).*xc);
BB=B';
A(1,:)= B(1,:)/P2;
Vdot(1,:)=0;
V2dot(1,:)=0;
Vdotzegx(1,:)=0;
elseif c>=2
V(c,:)=0.5*dt^2*A(c-1,:)*P2+V(c-1,:)+dt*Vdot(c-1,:)+...
(fi(c)-fi(c-1)-dt*fidot(c-1)+(si(c)-si(c-1)-...
dt*sidot(c-1)).*xc);
Vdotzegx(c,:)=dt*A(c-1,:)*H+Vdotzegx(c-1,:);
Vzegx(c,:)=0.5*dt^2*A(c-1,:)*H+Vzegx(c-1,:)+dt*Vdotzegx(c-1,:);
B(c,:)=(1/L^2*Vzegx(c,:)-sin(V(c,:))-fi2dot(c+1).*ones(1,2*M)-...
- si2dot(c+1).*xc);
A(c,:)=B(c,:)/P2;
Vdot(c,:)=dt*A(c-1,:)*P2+Vdot(c-1,:)+(fidot(c)-fidot(c-1)).*ones(1,2*M)+(sidot(c)-sidot(c-1)).*xc;
end
end
% dd=[1,100]
% for d=1:length(dd)
%
% hold on
plot(xc,V(1000,:))
% end
% error=(B(1000,:)/P2)*P2-B(1000,:);
% cond(A(1000,:))

Answers (0)

Products


Release

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by