how can i run this code need help?
古いコメントを表示
%its about Solving PDEs by Asymmetric MQ Collocation the above code is already given in the book which %i have
%pasted below i havn't much experience in MATALAB thats why i shared it here...
%Neumann Boundary Conditions
%The function f is set to f(x, y) = -2(2y3 - 3y2 + 1) + 6(1 - x2)(2y - 1) and Dirichlet boundary conditions
%of u(0, y) = 2y2 - 3y2 + 1 and u(1, y) = 0 are applied as well as Neumann boundary conditions of
%∂u/∂y = 0, along y = 0 and y = 1.
%A shape parameter of ε = 2.5 results in an evaluation matrix to be inverted that has a very high condition %number
%κ(H) ≈ 4.6e19.
clc
Nb = 22; % number of boundary point s
Np = 60; % t o t a l number of c ent e r , i n t e r i o r and boundary
centers = dlmread(centersUnitCircle60.txt ,'_');
x = centers(:,1);y = centers(:,2);
u = 65./(65+(x-0.2).^2+(y+0.1).^2 );% e x ac t s o l u t i o n
f = 130./(65+(x-0.2).^2+(y+0.1).^2 ).^3.*(2.*x-0.4).^2 - ...
260./(65+(x-0.2).^2+(y+0.1).^2).^2+...
130./(65+(x-0.2).^2+(y+0.1).^2).^3.*(2.*y+0.2).^2;
H = zeros (Np,Np ); rx = zeros (Np,Np ); ry = zeros (Np,Np); r = zeros (Np,Np );
f ( 1 :Nb) = u ( 1 :Nb ) ; % Di r i c h l e t Boundary c ond i t i ons
for i =1:Np
for j =1:Np
rx ( i , j ) = x ( i ) -x ( j ) ;
ry ( i , j ) = y ( i ) - y ( j ) ;
r ( i , j ) = sqrt ( rx ( i , j )^2 + ry ( i , j )^2 ) ;
end
end
index = 1 ;
for shape =1.0: -0.01: 0.1;
H( 1 :Nb , : ) = mq( r ( 1 :Nb , : ) , shape ) ; % enf or c e boudary c ond i t i ons
Hxx = mqDerivatives( r (Nb+1:Np , : ) , rx (Nb+1:Np , : ) , shape , 2 ) ; % enf or c e PDE
Hyy = mqDerivatives( r (Nb+1:Np , : ) , ry (Nb+1:Np , : ) , shape , 2 ) ;
H(Nb+1:Np , : ) = Hxx + Hyy ; % e v a l u a t i on matrix
kappa ( index ) = cond(H) ;
warning off
lambda = H\ f ; % expansion c o e f f i c i e n t s v i a Gaussian e l imi na t i on
warning on
B = mq( r , shape ) ; % system matrix
uh = B*lambda ; % approximate s o l u t i o n
e r ( index ) = norm(u?uh , i n f ) ;
sh ( index ) = shape ;
index = index + 1 ;
end
semilogy( sh , er , 'b' ) , xlabel ’ shape parameter ’ , ylabel ’max e r r o r ’
figure , semilogy( sh , kappa ) , xlabel ’ shape parameter ’ , ylabel ’ \kappa (H) ’
14 件のコメント
Walter Roberson
2022 年 3 月 9 日
what problems are you having? We do not have your file so we cannot test your code.
KSSV
2022 年 3 月 9 日
This line:
centers = dlmread(centersUnitCircle60.txt ,'_');
should be:
centers = dlmread('centersUnitCircle60.txt');
Akhtar Jan
2022 年 3 月 9 日
Walter Roberson
2022 年 3 月 9 日
We do not have your txt file, so we cannot run the code to test it. You will need to show us the error message instead of just telling us that it does not run.
@Akhtar Jan: Please do do explain, what is not happening, but explain, what you observe. "still not run" does not explain anything. You do have a useful error message on the screen already, so share it with the readers.
The code is full of typos. Did you copy&paste it from a PDF? E.g.
e r ( index ) = norm(u?uh , i n f ) ;
% should be:
er(index) = norm(uh , inf);
Akhtar Jan
2022 年 3 月 9 日
Jan
2022 年 3 月 9 日
Check what? It is your turn either to post some code and explain clearly, which problem you have with it, or to explain clearly, what you want to achieve and showing, what you have tried so far.
Currently we see a massively distorted text, which looks almost like Matlab, but the readers have to guess many details. Lines like:
e r ( index ) = norm(u?uh , i n f ) ;
are pure nonsense and you cannot expect, that readers in the internet can magically guess, what its intention is.
Walter Roberson
2022 年 3 月 9 日
編集済み: Walter Roberson
2022 年 3 月 9 日
clc
Nb = 22; % number of boundary point s
Np = 60; % t o t a l number of c ent e r , i n t e r i o r and boundary
centers = dlmread('centersUnitCircle60.txt' ,'_');
x = centers(:,1);y = centers(:,2);
u = 65./(65+(x-0.2).^2+(y+0.1).^2 );% e x ac t s o l u t i o n
f = 130./(65+(x-0.2).^2+(y+0.1).^2 ).^3.*(2.*x-0.4).^2 - ...
260./(65+(x-0.2).^2+(y+0.1).^2).^2+...
130./(65+(x-0.2).^2+(y+0.1).^2).^3.*(2.*y+0.2).^2;
H = zeros (Np,Np ); rx = zeros (Np,Np ); ry = zeros (Np,Np); r = zeros (Np,Np );
f ( 1 :Nb) = u ( 1 :Nb ) ; % Di r i c h l e t Boundary c ond i t i ons
for i =1:Np
for j =1:Np
rx ( i , j ) = x ( i ) -x ( j ) ;
ry ( i , j ) = y ( i ) - y ( j ) ;
r ( i , j ) = sqrt ( rx ( i , j )^2 + ry ( i , j )^2 ) ;
end
end
index = 1 ;
for shape =1.0: -0.01: 0.1;
H( 1 :Nb , : ) = mq( r ( 1 :Nb , : ) , shape ) ; % enf or c e boudary c ond i t i ons
Hxx = mqDerivatives( r (Nb+1:Np , : ) , rx (Nb+1:Np , : ) , shape , 2 ) ; % enf or c e PDE
Hyy = mqDerivatives( r (Nb+1:Np , : ) , ry (Nb+1:Np , : ) , shape , 2 ) ;
H(Nb+1:Np , : ) = Hxx + Hyy ; % e v a l u a t i on matrix
kappa ( index ) = cond(H) ;
warning off
lambda = H\f ; % expansion c o e f f i c i e n t s v i a Gaussian e l imi na t i on
warning on
B = mq( r , shape ) ; % system matrix
uh = B*lambda ; % approximate s o l u t i o n
er( index ) = norm(u - uh , inf ) ;
sh( index ) = shape ;
index = index + 1 ;
end
semilogy( sh , er , 'b' ) , xlabel 'shape parameter' , ylabel 'max error'
figure , semilogy( sh , kappa ) , xlabel 'shape parameter', ylabel '\kappa (H)'
Walter Roberson
2022 年 3 月 9 日
編集済み: Walter Roberson
2022 年 3 月 9 日
shape = 2.5;
N = 30;
Np = N^2;
[X,Y] = meshgrid(linspace(0,1,N),linspace(0,1,N));
x = reshape(X,N^2,1);
y = reshape(Y,N^2,1);
f = -2*(2*y.^3-3*y.^2+1)+6*(1-x.^2).*(2*y-1);
u = (1-x.^2).*(2*y.^3-3*y.^2+1); % exact solution
H = zeros(Np,Np);
rx = zeros(Np,Np);
ry = zeros(Np,Np);
r = zeros(Np,Np);
dirichletBCs = find(x==0 | x==1); % identify boundry and interior centers
neumannBCs = find ((y==0 | y==1)& ~(x==0 | x==1));
interior = find(x~=0 & x ~ = 1 & y~ = 0 & y~ = 1 );
f(dirichletBCs)=u(dirichletBCs);
f(neumannBCs)= 0;
for i =1:Np
for j =1:Np
rx( i , j ) = x( i )-x ( j ) ;
ry( i , j ) = y( i )-y ( j ) ;
r( i , j ) = sqrt( rx( i , j )^2 + ry( i , j )^2 );
end
end
H(dirichletBCs,:)=mq(r(dirichletBCs,:),shape); % e v a l u a t i on matrix
H(neumannBCs,:)= mqDerivatives(r(neumannBCs,:),ry(neumannBCs,:),shape,1);
Hxx = mqDerivatives(r(interior,:),rx(interior,:),shape,2);
Hyy = mqDerivatives(r(interior,:),ry(interior,:),shape,2);
H(interior,:)=Hxx+Hyy;
kappa=cond(H);
alpha=H\f;
B = mq( r , shape ) ; % system matrix
uh = B*alpha;
error = norm(uh-u,inf);
t = delaunay(x,y);trisurf(t,x,y,abs(uh-u));
xlabel 'x' , ylabel 'y' colormap('Summer')
Walter Roberson
2022 年 3 月 9 日
Neither of those are tested, as we do not have the data file and we do not have functions mq or mqDerivative
Akhtar Jan
2022 年 3 月 10 日
編集済み: Akhtar Jan
2022 年 3 月 10 日
Jan
2022 年 3 月 10 日
@Walter Roberson: Although the OP does not want to do it, I want to say: Thank you for spending your time for fixing the code.
回答 (1 件)
Image Analyst
2022 年 3 月 9 日
0 投票
The function is called 'linspace', not 'linespace'.
4 件のコメント
Image Analyst
2022 年 3 月 10 日
And many of your spaces are not allowed, like
close a l l
and
0 . 2
and
i n f
etc.
Akhtar Jan
2022 年 3 月 10 日
Akhtar Jan
2022 年 3 月 10 日
@Akhtar Jan: Your code contains to many spaces. Your description contains a lot of % characters. It is a good idea to spend the time to post clean messages, if you want to motivate others to help you.
You have been told already, that statements as "but not.....run" are not useful. Post the error message you see, or even better: read the message and try to fix the specific problem.
カテゴリ
ヘルプ センター および File Exchange で Partial Differential Equation Toolbox についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!