solve pde with variable coefficent
1 view (last 30 days)
Show older comments
How can i discretize my PDE if it contains a coefficient that include the func. u(x,t) ?
PDE
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/535164/image.png)
where e.g.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/534664/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/534669/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/534674/image.png)
so i tried to use explicit methode until here; assuming ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/534679/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/534679/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/535169/image.png)
this lee equation from this
but how can i code the a, since i dont know the u(x,t) and solve the pde numerically ?
%% FTCS- ITeration
M=40;
K=100;
a=0;
b=1;
h=(b-a)/M;
r=0.49;
k=r*h^2;
x=linspace(a,b,M+1);
t=0:k:K*k; % timestep
U=zeros(K+1,M+1);
uL=0;uR=0; % BC
uIC=sin(4*pi*x); %IC
U(1,:)=uIC; U(:,1)=uL; U(:,M+1)=uL;
a_u=ones(K+1,M+1); %a is a constant of 1
for kk=2:K+1
for jj=2:M
a_p=a_u(kk-1,jj);a_e=a_u(kk-1,jj-1);a_w=a_u(kk-1,jj+1);
U_W=U(kk-1,jj-1);U_E=U(kk-1,jj+1);U_P=U(kk-1,jj);
U(kk,jj)= r/4*(a_w-a_e)*(U_W-U_E)+r*a_p*U_W+r*a_p*U_E+(1-2*r*a_p)*U_P;
end
U(kk,M+1)=r/4*(a_w-a_e)*(U_W-U_E)+r*a_p*U(kk-1,M)+(1-2*r*a_p)*U(kk-1,M+1);
U(kk,1)=r/4*(a_w-a_e)*(U_W-U_E)+r*a_p*U(kk-1,2)+(1-2*r*a_p)*U(kk-1,1);
end
% plot
figure
plot(x,U(1:end,1:end),'linewidth',2);
a1 = ylabel('Y');
set(a1,'Fontsize',14);
a1 = xlabel('x');
set(a1,'Fontsize',14);
a1=title(['FTCS Method - r =' num2str(r)]);
legend('explicite')
set(a1,'Fontsize',16);
xlim([0 1]);
grid;
figure
[X, T] = meshgrid(x,t);
s2=surf(X,T,U,'FaceAlpha',0.5,'FaceColor','interp');hold on;
title(['3-D plot FTCS - r = ' num2str(r)])
%set(s2,'FaceColor',[1 0 1],'edgecolor',[0 0 0],'LineStyle','--');
xlabel('x= Length[mm]');ylabel('t= time [s]');zlabel('T = Temperature [°C]');
2 Comments
Answers (0)
See Also
Categories
Find more on Boundary Conditions 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!