Please someone should help me debug my error. The approximate solution is supposed to converge to the exact solution. I think i have problem with the code.

1 visualización (últimos 30 días)
i want to plot a graph of the numerical and exact solution. however, the errore is too much. can someone help me to trace my error. the solution are supposed to be ALMOST THE SAME. Here is my code:
alpha=0.25;
r=5;
t=0.002;
x=[-5:0.2:5];
k = reshape(0:100, 1, 1, []);
A=exp(x);
B = (t.^(k.*alpha))./factorial(k).*alpha.^k;
C = sum(B,3);
Numerical=A.*C;
Exact=A.*exp((r-4).*(t.^alpha)./(alpha));
Z1=Numerical;
Z2=Exact;
hZ1 = plot(x,Z1,'-r');
hold on
hZ2 = plot(x,Z2,'-g');
hold off
grid on
xlabel('x')
ylabel('u')
legend([hZ1(1),hZ2(1)], 'CFRDTM','EXACT', 'Location','NE')

Respuesta aceptada

Star Strider
Star Strider el 21 de Ag. de 2021
There are several typographical errors, at lelast with respect to the code representing the symbolic equations in the image.
With those corrections, ‘Numerical’ and ‘Exact’ are the same (within floating-point approximation error), at least for the values provided.
alpha=0.25;
r=5;
t=0.002;
x=[-5:0.2:5];
k = reshape(0:100, 1, 1, []);
A=exp(x);
% B =(t.^(k.*alpha))./(factorial(k).*alpha.^k);
B = (t.^(k.*alpha)).*(r-4).^k ./ (factorial(k).*alpha.^k);
C = sum(B,3);
Numerical=A.*C;
% Exact=A.*exp((r-4).*(t.^alpha)./(alpha));
Exact = exp(x+((r-4).*(t.^alpha)./(alpha)));
Z1=Numerical
Z1 = 1×51
0.0157 0.0192 0.0234 0.0286 0.0349 0.0427 0.0521 0.0637 0.0778 0.0950 0.1160 0.1417 0.1731 0.2114 0.2582 0.3153 0.3852 0.4704 0.5746 0.7018 0.8572 1.0470 1.2788 1.5619 1.9077 2.3301 2.8460 3.4761 4.2457 5.1857
Z2=Exact
Z2 = 1×51
0.0157 0.0192 0.0234 0.0286 0.0349 0.0427 0.0521 0.0637 0.0778 0.0950 0.1160 0.1417 0.1731 0.2114 0.2582 0.3153 0.3852 0.4704 0.5746 0.7018 0.8572 1.0470 1.2788 1.5619 1.9077 2.3301 2.8460 3.4761 4.2457 5.1857
Error = Numerical - Exact;
ErrorRMS = sqrt(mean(Error.^2))
ErrorRMS = 1.1945e-14
figure
hZ1 = plot(x,Z1,'-r');
hold on
hZ2 = plot(x,Z2,'--g'); % Change To Dashed Line
hold off
grid on
xlabel('x')
ylabel('u')
legend([hZ1(1),hZ2(1)], 'CFRDTM','EXACT', 'Location','NE')
It would also help if the code was parsed into different lines (as I had to do here), in order for it to run.
.
  5 comentarios
Star Strider
Star Strider el 21 de Ag. de 2021
The function is part of the code I posted.
Please copy the entire code, then run it, as I do here —
x = [-4,-3,-2,-1,2,3,4];
t = [0.002,0.004,0.006,0.008];
[Xm,Tm] = ndgrid(x,t);
[Exact,Numerical] = runCode(Xm(:),Tm(:));
absErr = abs(Exact-Numerical);
Results = table(Xm(:), Tm(:), Exact, Numerical, absErr, 'VariableNames',{'x','t','Exact','Numerical','AbsoluteError'})
Results = 28×5 table
x t Exact Numerical AbsoluteError __ _____ ________ _________ _____________ -4 0.002 0.042677 0.042677 1.3878e-17 -3 0.002 0.11601 0.11601 2.7756e-17 -2 0.002 0.31534 0.31534 5.5511e-17 -1 0.002 0.85718 0.85718 1.1102e-16 2 0.002 17.217 17.217 3.5527e-15 3 0.002 46.801 46.801 1.4211e-14 4 0.002 127.22 127.22 4.2633e-14 -4 0.004 0.050084 0.050084 1.3878e-17 -3 0.004 0.13614 0.13614 2.7756e-17 -2 0.004 0.37007 0.37007 0 -1 0.004 1.006 1.006 0 2 0.004 20.205 20.205 0 3 0.004 54.924 54.924 1.4211e-14 4 0.004 149.3 149.3 2.8422e-14 -4 0.006 0.055758 0.055758 0 -3 0.006 0.15157 0.15157 0
figure
scatter3(Results.x, Results.t, Results.Exact, 'xb')
hold on
scatter3(Results.x, Results.t, Results.Numerical, '+r')
hold off
grid on
xlabel('x')
ylabel('t')
set(gca, 'ZScale','log')
legend('Exact', 'Numerical', 'Location','bestoutside')
function [Exact,Numerical] = runCode(x,t)
alpha=0.25;
r=5;
% t=0.002;
% x=[-5:0.2:5];
k = reshape(0:100, 1, 1, []);
A=exp(x);
% B =(t.^(k.*alpha))./(factorial(k).*alpha.^k);
B = (t.^(k.*alpha)).*(r-4).^k ./ (factorial(k).*alpha.^k);
C = sum(B,3);
Numerical=A.*C;
% Exact=A.*exp((r-4).*(t.^alpha)./(alpha));
Exact = exp(x+((r-4).*(t.^alpha)./(alpha)));
end
.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Resizing and Reshaping Matrices en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by