I am having error in this programme. i gave the call function also.

28 ビュー (過去 30 日間)
GEETHA NANIAPPAN
GEETHA NANIAPPAN 2021 年 9 月 15 日
コメント済み: Walter Roberson 2021 年 9 月 15 日
Nx=41;XL=0;XR=2*pi;
dx=(XR-XL)/(Nx-1);
for i=1:Nx
XX(i)=XL+(i-1)*dx;
end
cfl=0.125;dt=cfl*dx^3;tlast=1;
Nstep=tlast/dt;
cc=1;
IRK=4;
[u0ex]=feval(@JSC1,XX,0,cc);
uold=u0ex;
subplot(2,1,1),plot(XX,u0ex,'k-');
hold on
for k=1:Nstep
if IRK==4
u0x=reconuxxxp(uold,Nx,dx);
k0=dt*resJSC1(uold,u0x,cc);
u1=uold+k0/2;
u1X=reconuxxxp(u1,Nx,dx);
k1=dt*resJSC1(uold,u1x,cc);
u2=uold+k1/2;
u2x=reconuxxxp(u2,Nx,dx);
k2=dt*resJSC1(uold,u2x,cc);
u3=uold+k2;
u3x=reconuxxxp(u3,Nx,dx);
k3=dt*resJSC1(uold,u3x,cc);
unew=uold+(k0+2*k1+2*k2+k3)/6;
elseif IRK==2
u0x=reconuxp(uold,Nx,dx);
u0xx=reconuxp(u0x,Nx,dx);
u0xxx=reconuxp(u0xx,Nx,dx);
u2x=reconuxp(uold.^2,Nx,dx);
ko=dt*resJSC1(u2x,u0xxx,cc);
u1=uold+k0;
u1x=reconuxp(u1,Nx,dx);
uxx=reconuxp(u1x,Nx,dx);
uxxx=reconuxp(uxx,Nx,dx);
u2x=reconuxp(u1.^2,Nx,dx);
k1=dt*resJSC1(u2x,uxxx,cc);
u2=u1+k1;
unew=(uold+u2)/2;
end
uold=unew;
eps=0.3*dt;
if abs(k*dt-0.25)<eps|abs(k*dt-0.5)<eps|abs(k*dt-0.75)<eps|abs(k*dt-1)<eps disp(k),
[u0ex]=feval(@JSC1,XX,k*dt,cc);
subplot(2,1,1),plot(XX,unew,'b-',XX,uoex,'g.');
xlabel('x');ylabel('u(x,t)');
hold on
subplot(2,1,2),
if abs(k*dt-0.25)<eps,plot(XX,abs(u0ex-unew),'r-');end
if abs(k*dt-0.5)<eps,plot(XX,abs(u0ex-unew),'r-.');end
if abs(k*dt-0.75)<eps,plot(XX,abs(u0ex-unew),'r--');end
if abs(k*dt-1)<eps,plot(xx,abs(u0ex-unew),'r:');end
xlabel('x');ylabel('inf error');
hold on
end
end
function uu=JSC1(x,ap,cc)
uu=sin(cc*(x+ap));
function uu=resJSC1(u1,u2,cc)
uu=-cc*u2;
function ux=reconuxxxp(u,N,h)
IEX=0;
if IEX==4
a=2;b=-1;
for i=1:N
if i==1
tmp=b*(u(i+3)-3*u(i+1)+3*u(N-1)-u(N-3))/4+a*(u(i+2)-2*u(i+1)+2*u(N-1)-u(N-2));
else if i==2
tmp=b*(u(i+3)-3*u(i+1)+3*u(i-1)-u(N-2))/4+a*(u(i+2)-2*u(i+1)+2*u(i-1)-u(N-1));
else if i==3
tmp=b*(u(i+3)-3*u(i+1)+3*u(i-1)-u(N-1))/4+a*(u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2));
else if i==(N-2)
tmp=b*(u(2)-3*u(i+1)+3*u(i-1)-u(i-3))/4+a*(u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2));
else if i==(N-1)
tmp = b*(u(3)-3*u(i+1)+3*u(i-1)-u(i-3))/4+a*(u(2)-2*u(i+1)+2*u(i-1)-u(i-2));
else if i==N
tmp=b*(u(4)-3*u(2)+3*u(i-1)-u(i-3))/4+a*(u(3)-2*u(2)+2*u(i-1)-u(i-2));
else
tmp=b*(u(i+3)-3*u(i+1)+3*u(i-1)-u(i-3))/4+a*(u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2));
end
ux(i)=tmp/(2*h^3);
end
return
end
if IEX==6
a=-488/240;b=338/240;c=-72/240;d=7/240;
for i=1:N
if i==1
tmp=a*(u(i+1)-u(N-1))+b*(u(i+2)-u(N-2))+c*(u(i+3)-u(N-3))+d*(u(i+4)-u(N-4));
else if i==2
tmp=a*(u(i+1)-u(i-1))+b*(u(i+2)-u(N-1))+c*(u(i+3)-u(N-2))+d*(u(i+4)-u(N-3));
else if i==3
tmp=a*(u(i+1)-u(i-1))+b*(u(i+2)-u(i-2))+c*(u(i+3)-u(N-1))+d(u(i+4)-u(N-2));
else if i==4
tmp=a*(u(i+1)-u(i-1))+b*(u(i+2)-u(i-2))+c*(u(i+3)-u(i-3))+d*(u(i+4)-u(N-1));
else if i==(N-3)
tmp=a*(u(i+1)-u(i-1))+b*(u(i+2)-u(i-2))+c*(u(i+3)-u(i-3))+d(u(2)-u(i-4));
else if i==(N-2)
tmp=a*(u(i+1)-u(i-1))+b*(u(i+2)-u(i-2))+c*(u(2)-u(i-3))+d*(u(3)-u(i-4));
else if i==(N-1)
tmp=a*(u(i+1)-u(i-1))+b*(u(2)-u(i-2))+c*(u(3)-u(i-3))+d*(u(4)-u(i-4));
else if i==N
tmp=a*(u(2)-u(i-1))+b*(u(3)-u(i-2))+c*(u(4)-u(i-3))+d*(u(5)-u(i-4));
else
tmp=a*(u(i+1)-u(i-1))+b*(u(i+2)-u(i-2))+c*(u(i+3)-u(i-3))+d*(u(i+4)-u(i-4));
end
ux(i)=tmp/(h^3);
end
return
end
ISC=6;
if ISC==4
alfa=15/32;aa=2;bb=2*alfa-1;
else if ISC==6
alfa=7/16;aa=2;bb=-1.0/8;
end
amat=zeros(N,N);
B=zeros(N,1);
for i=1:N
if i==1
amat(1,1)=1;amat(1,2)=alfa;amat(1,N-1)=alfa;
B(i)=bb*(u(i+3)-3*u(i+1)+3*u(N-1)-u(N-3))/4+aa*(u(i+2)-2*u(i+1)+2*u(N-1)-u(N-2));
else if i==2
amat(2,1)=alfa;amat(2,2)=1;amat(2,3)=alfa;
B(i)=bb*(u(i+3)-3*u(i+1)+3*u(i-1)-u(N-2))/4+aa*(u(i+2)-2*u(i+1)+2*u(i-1)-u(N-1));
else if i==3
amat(i,i-1)=alfa;amat(i,i)=1;amat(i,i+1)=alfa;
B(i)=bb*(u(i+3)-3*u(i+1)+3*u(i-1)-u(N-1))/4+aa*(u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2));
else if i==N-2
amat(i,i-1)=alfa;amat(i,i)=1;amat(i,i+1)=alfa;
B(i)=bb*(u(2)-3*u(i+1)+3*u(i-1)-u(i-3))/4+aa*(u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2));
else if i==N-1
amat(i,i-1)=alfa;amat(i,i)=1;amat(i,i+1)=alfa;
B(i)=bb*(u(3)-3*u(i+1)+3*u(i-1)-u(i-3))/4+aa*(u(2)-2*u(i+1)+2*u(i-1)-u(i-2));
else if i==N
amat(N,1)=-1;amat(N,N)=1;
B(i)=0
else
amat(i,i-1)=alfa;amat(i,i)=1;amat(i,i+1)=alfa;
B(i)=bb*(u(i+3)-3*u(i+1)+3*u(i-1)-u(i-3))/4+aa*(u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2));
end
B(i)=B(i)/(2*h^3);
ux=(amat\B)';
end
end
end
end
end
end
end
end
end
end
end
end
end
end
end
end
end
end
end
end
end
end
  1 件のコメント
KSSV
KSSV 2021 年 9 月 15 日
What error you are getting? Your ocde has got lot of loops. I think there is a lot of room to make your code effective.

サインインしてコメントする。

サインインしてこの質問に回答する。

回答 (1 件)

Walter Roberson
Walter Roberson 2021 年 9 月 15 日
The function you are missing is on PDF page 108 of https://www.researchgate.net/file.PostFileLoader.html?id=5696da346143257e508b4567&assetKey=AS%3A317578013020160%401452727906601
  2 件のコメント
Star Strider
Star Strider 2021 年 9 月 15 日
+1
... because I have absolutely no idea how you figured that out!
.
Walter Roberson
Walter Roberson 2021 年 9 月 15 日
just googled reconuxp and looked at the first few results.

サインインしてコメントする。

カテゴリ

Help Center および File Exchange で Loops and Conditional Statements についてさらに検索

タグ

製品

Community Treasure Hunt

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

Start Hunting!

Translated by