# Sagar an Payne stress-strain relationships and Boussinesq aproximation in Matlab

I have to do 2 problems in Matlab, and the Math course is not my favourite one. However, I have tried to resolve the first problem, based on another exercise, but I'm pretty sure it's wrong. Can you please tell me what's wrong with it and correct me or lead me to the correct result? On the second problem, I know I should resolve it with some kind of interpolation or the small squares method (or something like this). But as far as I know, these methods are used with polynomials, and I can't see one in the problem.

Please, any explanation would be great, hopefully I'll understand something.

What I've tried:

dt=0.0013; dx=0.05;
Nx=1/(dx)+1;
x=0:dx:1;
tf=0.05;
Nt=tf/dt;
Uo=zeros(Nx,1);
Un=zeros(Nx,1);
K=0.1;

%initialization
for i=1:round(Nx/2)
Uo(i)=(i-1)*dx;
end
for i=round(Nx/2):Nx
Uo(i)=1-(i-1)*dx;
end

n=0;
while (n<Nt)
n=n+1;
for i=2:Nx-1
%Un(i)=Uo(i)+dt/(dx*dx)*(Uo(i-1)-2*Uo(i)+Uo(i+1));
Un(i) = Uo(i) + (4*K * dt * (x * (Uo(i-1)-2*Uo(i)+Uo(i+1)) + dx * (Uo(i+1) - Uo(i-1))))/(x * dx * dx);
end
Uo=Un;
end
plot(x,Uo,'.-r')
hold on
%analytic solution
U=HeatAnalytic(x,tf);
plot(x,U,'b')


The HeatAnalytic function:

function rez=HeatAnalytic(x,t)
rez=0;
for k=1:100
rez=rez+4/(k*pi)^2*sin(k*pi/2)*sin(k*pi*x)*exp(-k^2*pi^2*t);
end


The second problem:

Given the table with the water density variation with temperature: T 0 4 10 15 20 22 25 30 40
60 80 100 ρ 999.84 999.97 999.7 999.1 998.21 997.77 997.05 995.65 992.2 983.2 971.8 958.4 T -> Celsius degrees and ρ -> kg/m^3

And the Boussinesq approximation: ρ = ρ0 * (1 - β*(T −T0)), where the temperature is expressed in Kelvin (T(C) = T(K) + 273.15), and ρ0 and T0 represents the values at 0 degrees. Find the β coefficient.

-

Regarding the first problem:

The initialization is wrong, it should contain the initial and boundary conditions from the problem statement. The for loop is correct, but it does not have an analytic solution, so all the HeatAnalytic part shouldn't be there.

About the second problem, I managed to solve it, using the method of least squares. This is the Matlab code:

The least squares function:

function rez=least_squares(xtrimis,ytrimis,N);

N= N+1;
x=xtrimis;
y=ytrimis;
for i=1:N
for j=1:N
A(i,j) = sum(x.^(i+j -2));
end
C(i)= dot(y,x.^(i-1));
end

A;
D=C';
B = A\D;

B=flipud(B);
rez=B;


And the main part:

clc;
clear;

T0 = 273.15;
ro0 = 999.84;

T = [0 4 10 15 20 22 25 30 40 60 80 100] ;
ro = [999.84 999.97 999.7 999.1 998.21 997.77 997.05 995.65 992.2 983.2 971.8 958.4];

x = T + T0;
y = ro;

polin = E28P2_fPatratele(x,y,1);

a = polin(1)
b = polin(2)

bet = (b - ro0)/T0

plot(T + T0,ro, 'o')
hold on;

t1= 273.15:373.15;

ro1 = (ro0 + bet * T0) - (ro0 * bet * t1);
plot(t1,ro1, '-r');

-

## protected by T. Bongers May 26 '14 at 8:13

Thank you for your interest in this question. Because it has attracted low-quality or spam answers that had to be removed, posting an answer now requires 10 reputation on this site.