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 Ansicht (letzte 30 Tage)
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')

Akzeptierte Antwort

Star Strider
Star Strider am 21 Aug. 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 Kommentare
Star Strider
Star Strider am 21 Aug. 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
.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu MATLAB finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by