How to calculate u for every alpha?

Hi Guys,

In my project i have to Write a script that solves the heat conduction equation only for linear basis

functions with Robin boundary conditions. This time different values of α

are to be considered in a simulation. I have to determine for each α the function value

at the left and right end of the interval. Store these values in a matrix. For ¨

α = 2, 5, 10, 15, 20, for example, the matrix could look like the following.

I know how to calculate it for one alpha like alpha = 2, but i dont know how to do it with more alpha and build the matrix. With the for loop i only get the values for alpha=20... This is the code:

...
f= @(x)(0*zeros(size(x)));
alpha = [2,5,10,15,20];

%lineare Basisfunktion
[A]= assemDiffusion1D(nodes, epsilon);
b = assemRightSide1D(f,nodes, "exact");

%Robin RW
for alpha=[2,5,10,15,20]
n = length(A);
b(1)=b(1)+epsilon*alpha/400*271.15;
b(n)=b(n)+epsilon*alpha/400*ul;
A(1,1)=A(1,1)+epsilon*alpha/400;
A(n,n)=A(n,n)+epsilon*alpha/400;
ulinR= A\b;
end

T = [alpha ulinR(1) ulinR(n)]


end 

[Torsten]

For this harmless question you really need to change your Username ?

Answer 1

[Torsten]

...
f= @(x)(0*zeros(size(x)));
alpha = [2,5,10,15,20];
%lineare Basisfunktion
A = assemDiffusion1D(nodes, epsilon);
b = assemRightSide1D(f,nodes, "exact");
%Robin RW
T = zeros(numel(alpha),3);
for i = 1:numel(alpha)
    n = length(A);
    b(1)=b(1)+epsilon*alpha(i)/400*271.15;
    b(n)=b(n)+epsilon*alpha(i)/400*ul;
    A(1,1)=A(1,1)+epsilon*alpha(i)/400;
    A(n,n)=A(n,n)+epsilon*alpha(i)/400;
    ulinR= A\b;
    T(i,:) = [alpha(i),ulinR(1),ulinR(n)]
end
...

[Luluhihi2021]

Thank you so so much! It worked, but the values are not correct for example for alpha = 20 (not like in the example). With your code i get for alpha = 20 u(1) = 281.5421 and in the example it is u(1) = 281.9079... What could possibly be the reason?

[Torsten]

...
f= @(x)(0*zeros(size(x)));
alpha = [2,5,10,15,20];
%lineare Basisfunktion
AA = assemDiffusion1D(nodes, epsilon);
bb = assemRightSide1D(f,nodes, "exact");
%Robin RW
T = zeros(numel(alpha),3);
for i = 1:numel(alpha)
    A = AA;
    b = bb;
    n = length(A);
    b(1)=b(1)+epsilon*alpha(i)/400*271.15;
    b(n)=b(n)+epsilon*alpha(i)/400*ul;
    A(1,1)=A(1,1)+epsilon*alpha(i)/400;
    A(n,n)=A(n,n)+epsilon*alpha(i)/400;
    ulinR= A\b;
    T(i,:) = [alpha(i),ulinR(1),ulinR(n)]
end
...

[Luluhihi2021]

I appreciate it so much! Thank you !:)

[Luluhihi2021]

One last question for comprehension: What did you do by putting A= AA, and b = bb? Or why did you do that? Did u save it in another variable ?

[Torsten]

To save workspace, this might be a better solution:

...
f= @(x)(0*zeros(size(x)));
alpha = [2,5,10,15,20];
%lineare Basisfunktion
A = assemDiffusion1D(nodes, epsilon);
b = assemRightSide1D(f,nodes, "exact");
A11 = A(1,1);
Ann = A(n,n);
b1 = b(1);
bn = b(n);
%Robin RW
T = zeros(numel(alpha),3);
for i = 1:numel(alpha)
    n = length(A);
    b(1)=b1+epsilon*alpha(i)/400*271.15;
    b(n)=bn+epsilon*alpha(i)/400*ul;
    A(1,1)=A11+epsilon*alpha(i)/400;
    A(n,n)=Ann+epsilon*alpha(i)/400;
    ulinR= A\b;
    T(i,:) = [alpha(i),ulinR(1),ulinR(n)]
end
A(1,1) = A11;
A(n,n) = Ann;
b(1) = b1;
b(n) = bn;

...