making pde matlab code

1 回表示 (過去 30 日間)
Mr.DDWW
Mr.DDWW 2022 年 3 月 31 日
回答済み: Torsten 2022 年 4 月 1 日
I want to show high values of "Nu", the solution will converge, from the equtions in image
clc;clear all;close all;
L = 1;
x = linspace(0,L,100);
t = linspace(0,1,100);
% Because symmetry
m = 1;
sol = pdepe(m,@heatpde,@heatic,@heatbc,x,t);
sol1=sol;
figure(1)
figure
hold all % the use of all makes the colors cycle with each plot
for i=1:10:numel(t)
plot(x,sol1(i,:),'LineWidth',1.5)
end
% we use -DynamicLegend to pick up the display names we set above.
%legend('-DynamicLegend','location','bestoutside');
%surf(x,t,sol1);
xlabel('y/b');
zlabel('Phsi');
title('Fig 12.1-1'); grid on;
function [c,f,s] = heatpde(x,t,u,dudx)
c = 1;
f = dudx;
s = 0;
end
function u0 = heatic(~)
u0 = 1;
end
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t )
Nu = 1 ;
% left Bc = ul
pl = 0 ;
ql = 1 ;
% right BC= ur
pr = Nu*ur ;
qr = 1 ;
end

回答 (2 件)

Torsten
Torsten 2022 年 3 月 31 日
Setting
m = 1
does not implement your equation in the image.
Setting m = 1 would mean solving
d(theta)/d(tau) = 1/eta * d/d(eta) (eta * d(theta)/d(eta))
I want to show high values of "Nu", the solution will converge, from the equtions in image
And the code does not converge if you set Nu to a high value ?
  5 件のコメント
Torsten
Torsten 2022 年 3 月 31 日
What error message do you get from MATLAB ?
Mr.DDWW
Mr.DDWW 2022 年 3 月 31 日
I am not sure how to implement (show high values of "Nu", the solution will converge)

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


Torsten
Torsten 2022 年 4 月 1 日
Here is a simple code to check your results from pdepe:
% Spatial interval
L = 1;
nL = 101;
dL = L/(nL-1);
eta = linspace(0,L,nL);
% temporal interval
Tau = 1;
nTau = 101;
tau = linspace(0,Tau,nTau);
% Nusselt number
Nu = 10;
% Call solver
theta0 = ones(nL,1);
tauspan = tau;
[t, theta] = ode15s(@(tau,theta)fun(tau,theta,nL,dL,Nu),tauspan,theta0);
% Plot results Here: theta-profile over eta at tau = Tau
plot(eta,theta(end,:))
function dthetadtau = fun(tau,theta,nL,dL,Nu)
dthetadtau = zeros(nL,1);
% @eta = 0 : d(Theta)/d(eta) = 0 @ eta = 0
dthetadtau(1) = 2*(theta(2)-theta(1))/dL^2;
dthetadtau(2:nL-1) = (theta(3:nL) - 2*theta(2:nL-1) + theta(1:nL-2))/dL^2;
% @eta = 1: d(Theta)/d(eta) + Nu*Theta = 0
dthetadtau(nL) = ((theta(nL-1) - 2*dL*Nu*theta(nL)) - 2*theta(nL) + theta(nL-1))/dL^2;
end

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by