Euler method empty plot

I am trying to plot the vector h vs the x axis but I keep getting an empty plot when I use euler method but I do not get an error so I do not know where everything is failing. My H vector is full of NaN so I guess that is why I my plot is empty.

function h = forwardmethodprob2(t, I, A2, dt, dx, h_2, Qx, h0, k)
    for i = 1:length(t)
        g = dt.*Qx;
        g(1) = g(1) + (dt*h0)*k(1)/dx.^2;
        h_1 = (I + dt*A2)*h_2 + g;
        h_2 = h_1;
    end
        h = h_1;
end

Is there any issues with my Euler method?

Answer 1

[Walter Roberson]

%Variables
T = 1;
L = 1;
Nt = 1000;
dt = T/Nt;
Nx = 1000;
dx = L/Nx;
k1 = 10;
k2 = 0.1;
L12 = 0.1;
L1 = 0.7;
L2 = 0.2;
alpha = 200;
beta = 2;
theta = 20;
h0 = 100;
k = [];
c = 0.001;

%Vector t and x
t = 0:dt:T;
xvector = 0:dx:L;

%define k function
for i = 2:length(xvector)
   if xvector(i) <= L1
       k(i-1) = k1;
   elseif xvector(i) >= L1 + L12
       k(i-1) = k2;
   else 
       k(i-1) = k1+(k2-k1)*((xvector(i)-L1)/(L12));
   end
end

k = k';

%define Qx 
Qx = alpha.*exp(-beta.*xvector).*sin(theta.*xvector);
Qx(2) = Qx(2) + (k(2)/dx.^2)*100;
Qx = Qx(2:end);

A2 = diag(ones(Nx - 1, 1), 1) - 2 * eye(Nx) + diag(ones(Nx - 1, 1), -1); 
A2(end, end) = -1;
A2 = (k/dx^2).*A2 ;
A2 = A2 + diag(ones(Nx-1,1)*(-c/dx),1) + eye(Nx).*(c/dx);
I = eye(length(A2));

for i = 1:length(xvector) - 1 
    h_2(i) = 0;
end

h_2 = h_2';
Qx = Qx';

H = forwardmethodprob2(t, I, A2, dt, dx, h_2, Qx, h0, k);
figure;
size(xvector)

ans = 1×2

1 1001

size(H)

ans = 1×2

1000 1

nnz(isnan(xvector))

ans = 0

nnz(isnan(H))

ans = 1000

plot(xvector(2:end)', H);

function h = forwardmethodprob2(t, I, A2, dt, dx, h_2, Qx, h0, k)
    for i = 1:length(t)
        g = dt.*Qx;
        g(1) = g(1) + (dt*h0)*k(1)/dx.^2;
        h_1 = (I + dt*A2)*h_2 + g;
        h_2 = h_1;
    end
        h = h_1;
end

All of your H values are NaN. Nothing will be plotted.

Question: you loop i over length(t) but you never use t inside the loop. What is the point of passing in t if you are not going to use it?

You do not index anything by i inside your loop. Your dt does not change inside your loop, and your Qx does not change inside your loop, so the assignment to g is going to produce the same result each time. You do not change dt or h0 or k or dx inside your loop, so the modification go g(1) is going to produce the same result each iteration of the loop. You do not change I or dt or A2 or h_2 inside the loop, and we proved that g will be the same each iteration, so your h_1 is going to have the same result each iteration. Your h_1 is copied to h_2 so you are going to have the same h_2 for each iteration. You also never use h_2 afterwards.

So... your forwardmethodprob2() is going to produce the same output for each value of i in the for loop; what is the point of having the for loop ?