calculation of inverse matrix
3 views (last 30 days)
Show older comments
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,:))
0 Comments
Answers (0)
See Also
Categories
Find more on Discrete Multiresolution Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!