h=input('输入空间步长 h=');tau=input('输入时间步长 tau=');
r=tau/h^2;
m=1/h;n=1/tau;
i=1:m-1;k=1:n-1;
xi=0+1/m:1/m:1-1/m;tk=1/n:1/n:1;
ui0=sin(pi/4+xi*pi/2)';u00=sqrt(2)/2;
um0=sqrt(2)/2;
u0k=sqrt(2)/2*exp(-tk)';
umk=sqrt(2)/2*exp(-tk)';
fik=zeros(m-1,n);
for j=1:n
fik(:,j)=(exp(-tk(j))*sin(pi/2+...
xi*pi/2).*(pi^2/4*exp(-4*tk(j)).*sin(pi/4+xi*pi/2).^4-1))';
end
uik=zeros(m-1,n);uik2=ones(m-1,n-1);
uik(:,1)=(sin(pi/4+xi*pi/2)+tau*(sin(pi/4+...
xi*pi/2).^4).*(-pi^2/4*sin(pi/4+xi*pi/2))...
+tau*sin(pi/4+xi*pi/2).*(pi^2/4*sin(pi/4+xi*pi/2).^4-1))';
%first k=1%
A(1,1)=1+2*r*(uik(1,1)^4);A(1,2)=-r*(uik(1,1)^4);
B(1,1)=1-2*r*(uik(1,1)^4);B(1,2)=r*(uik(1,1)^4);
A(m-1,m-2)=-r*(uik(m-1,1)^4);
A(m-1,m-1)=1+2*r*(uik(m-1,1)^4);
B(m-1,m-2)=r*(uik(m-1,1)^4);
B(m-1,m-1)=1-2*r*(uik(m-1,1)^4);
for i=2:m-2
A(i,i-1)=-r*(uik(i,1)^4);
A(i,i)=1+2*r*(uik(i,1)^4);
A(i,i+1)=-r*(uik(i,1)^4);
B(i,i-1)=r*(uik(i,1)^4);
B(i,i)=1-2*r*(uik(i,1)^4);
B(i,i+1)=r*(uik(i,1)^4);
end
C=B*ui0+2*tau*fik(:,1)+[r*(uik(1,1)^4*(u0k(2)+u00));...
zeros(m-3,1);r*(uik(m-1,1))^4*(umk(2)+um0)];
%uik(:,2)=tridiag(A,C,m-1);
uik(:,2)=inv(A)*C;
%%%%%%%%%%%%%%%%
for k=2:n-1
A=zeros(m-1,m-1);B=zeros(m-1,m-1);
A(1,1)=1+2*r*(uik(1,k)^4);A(1,2)=-r*(uik(1,k)^4);
B(1,1)=1-2*r*(uik(1,k)^4);B(1,2)=r*(uik(1,k)^4);
A(m-1,m-2)=-r*(uik(m-1,k)^4);A(m-1,m-1)=1+2*r*(uik(m-1,k)^4);
B(m-1,m-2)=r*(uik(m-1,k)^4);B(m-1,m-1)=1-2*r*(uik(m-1,k)^4);
for i=2:m-2
A(i,i-1)=-r*(uik(i,k)^4);
A(i,i)=1+2*r*(uik(i,k)^4);
A(i,i+1)=-r*(uik(i,k)^4);
B(i,i-1)=r*(uik(i,k)^4);
B(i,i)=1-2*r*(uik(i,k)^4);
B(i,i+1)=r*(uik(i,k)^4);
end
C=B*uik(:,k-1)+2*tau*fik(:,k)+...
[r*(uik(1,k)^4)*(u0k(k+1)+u0k(k-1));zeros(m-3,1);...
r*(uik(m-1,k)^4)*(umk(k+1)+umk(k-1))];
uik(:,k+1)=inv(A)*C;
end
true=zeros(m-1,n);
for i=1:m-1
true(i,:)=exp(-tk).*(sin(pi/4+xi(i)*pi/2));
end
% e=max(true-uik);
% h=max(e)
% plot(e)
% plot(uik(:,90));
% hold on;
plot(true(:,end)-uik(:,end))
%surf(abs(uik-true))
fprintf('h=%.3f tau=%.3f\n',h,tau);
fprintf('k (x,t) 数值解 精确解 误差 \n');
for k=1:10
p=int16(k*0.1/tau);
if p<100
fprintf('%d (0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),...
true(m/2,p),abs(uik(m/2,p)-true(m/2,p)));
else
fprintf('%d (0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),...
true(m/2,p),abs(uik(m/2,p)-true(m/2,p)));
end
end
a=max(abs(uik-true));
E=max(a);
fprintf('Emax=%e',E);
MATLAB.rar_Jacobi迭代法_grass5tr_jacobi迭代_matlab 五点迭代_椭圆型方程
版权申诉
5星 · 超过95%的资源 75 浏览量
2022-09-21
01:41:09
上传
评论 1
收藏 4KB RAR 举报
周楷雯
- 粉丝: 74
- 资源: 1万+
最新资源
- STM32单片机FPGA毕设电路原理论文报告汽车电动助力转向单片机控制系统设计与试验研究
- STM32单片机FPGA毕设电路原理论文报告气压传感器神经网络算法及单片机实现
- STM32单片机FPGA毕设电路原理论文报告频率的测量在单片机设计中的应用
- 音频转码工具(用于将微信语音 amr 格式转换为 mp3 格式以便在 html5 的 audio 标签中进行播放).zip
- RDK-Web-Performance-Node
- STM32单片机FPGA毕设电路原理论文报告片式电容器浪涌及老化测试系统的设计与实现
- 一个简易的贪吃蛇小游戏(C语言作业).zip
- 国家游戏防沉迷系统接口(php)
- STM32单片机FPGA毕设电路原理论文报告片剂硬度测试仪的液晶显示界面设计
- STM32单片机FPGA毕设电路原理论文报告偏磁式消弧线圈的动态调谐装置
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈
评论1