==================Tham khảo giải bài tập thực hành ==========================

Bài 1:
theta=0:10:90;
a=asind((1+3*cosd(theta))./sqrt((1+3*cosd(theta)).^2+(3-3*sind(theta)).^2));
F=300*4.5*sind(theta)./(3*cosd(a-theta));
T=[theta', F'];
fprintf(' Theta(deg) F(lb) \n')
disp(T)

bài 2:

clear all; clc

format shortg
fd=fopen('result_file.dat', 'w');
v=350*1000/3600; h=1500; % 1km/h doi sang m/s
t=0:0.5:6;
H=h+v*t*sind(10);
X=v*t*cosd(10);
r=sqrt(H.^2+X.^2);
th=acosd(X./r);
Table=[t', th', r'];
disp(' t(s) theta (degree) r(m)')
disp(Table)
fprintf(fd,' t(s) theta (degree) r(m) \n');
fprintf(fd,'%10g %10g %10g \n', Table);

%===================================================
Bai 3: GPT bằng các phương pháp Euler, Rk2, Rk3, Rk4
clc, clear
h=0.02; a=0; b=2; % h is the step size, t=[a,b] is the domain
y0=1; % initial condition; in MATLAB indices start at 1
%F= @(t,y) 3.*exp(-t)-0.4*y; % solving equation dy/dt=3.*exp(-t)-0.4*y
%F=@(t, y) t.^3
F= @(t,y) y*t.^2-1.1*y;
y4= RK4(F, a, b, y0, h)';
plot(y4(1,:), y4(2,:), '*-b', 'Linewidth', 3);
y3= RK3(F, a, b, y0, h)'; hold on
plot(y3(1,:), y3(2,:), 'o-r', 'Linewidth', 1);
y2= RK2(F, a, b, y0, h)'; hold on
plot(y2(1,:), y2(2,:), '+-g', 'Linewidth', 1);
YY=Euler_(F, a, b, y0, h); hold on
plot(YY(1,:), YY(2,:), '+-g', 'Linewidth', 1);
legend('RK4', 'RK3', 'RK2', 'Euler')
title('dy/dt=y*t.^2-1.1*y'); xlabel('t'); ylabel('y(t)')
%------------------RK4---------------------------
function yy= RK4(F, a, b, y0, h)
t = a:h:b; % Computes t-array from a to b,
y(1) = y0;
for i=1:1:length(t)-1 % the number of elements of t
k1 = F(t(i),y(i));
k2 = F(t(i)+0.5*h,y(i)+0.5*h*k1);
k3 = F((t(i)+0.5*h),(y(i)+0.5*h*k2));
k4 = F((t(i)+h),(y(i)+h*k3));
y(i+1) = y(i) + (h/6)*(k1+2*k2+2*k3+k4); % main equation
end
t=t'; y=y';yy=[t, y];
return
end
%-------------RK3-----------------------------
function yy= RK3(F, a, b, y0, h)
t = a:h:b; % Computes t-array from a to b,
y(1) = y0;
for i=1:1:length(t)-1 % the number of elements of t
k1 = F(t(i),y(i));
k2 = F(t(i)+0.5*h,y(i)+0.5*h*k1);
k3 = F((t(i)+h),(y(i)-h*k1+2*k2*h));
y(i+1) = y(i) + (h/6)*(k1+4*k2+k3); % main equation
end
t=t'; y=y';yy=[t, y];
return
end
%------------------RK2------------------------
function yy= RK2(F, a, b, y0, h) % Heun method
t = a:h:b; % Computes t-array from a to b,
y(1) = y0;
for i=1:1:length(t)-1 % the number of elements of t
k1 = F(t(i),y(i));
k2 = F(t(i)+h,y(i)+h*k1);
y(i+1) = y(i) + (h/2)*(k1+k2); % main equation
end
t=t'; y=y';yy=[t, y];
return
end
%------------------Euler-------------------------
function yy= Euler_(F, a, b, y0, h) % Heun method
t = a:h:b; % Computes t-array from a to b,
y(1) = y0;
for i=1:length(t)-1 % the number of elements of t
k1 = F(t(i),y(i));
y(i+1) = y(i) + h*k1; % main equation
end
t=t'; y=y'; yy=[t, y]';
return
end
%=================================================
Bài 4 shooting method

clc, clear all% program shooting method for fireworks problem y''=-g+a*y^2
global a te ye h ss;
a = 0.01; % coeff of nonlinear acceleration
te = 5; % end time [sec]
ye = 40; % end height [m]
h = 0.1; % time step
clf;
hold on;
ss=[];
for dy0=19:4:23
[tv,yv] = euler2_(h,0,te,0,dy0);
plot(tv,yv,'o-', 'LineWidth',2);
ss=[ss; dy0, yv(te/h+1)-ye];
end
dy = secant1_(20,20.7,1e-3);
% draw last curve
[tv,yv] = euler2_(h,0,te,0,dy);
plot(tv,yv,'o-', 'LineWidth',2);
% invariant: tv(te/h+1)==te
y = yv(te/h+1); % returns y at t=te
text(te+.2,y,sprintf('y\047(0)=%g',dy), 'FontSize',15);
text(te+.2,y,sprintf('y\047(0)=%g',dy), 'Color','r', 'FontSize',15);

set(gca, 'FontSize', 16); % for tick marks

line([te te],[0 60], 'Color','k');
line([te te],[ye ye], 'Color','r', 'Marker','o', 'LineWidth',3);
xlabel('t(s)', 'FontSize',20);
ylabel('y(m)', 'FontSize',20);
title(sprintf('Shooting Method on y\047\047=-g+0.01y^2'),
'FontSize',20);
figure(2), plot(ss(:,1), ss(:,2), '*-k')
%ylabel('Error(m))'); xlabel('t(s)')
%=============================================================
function x = secant1_(x1,x2,tol)
% secant method for one-dimensional root finding
global ye h te;
[tv,yv] = euler2_(h,0,te,0,x1); % x1 is the dy0, that is unknown and
assume It is x1)
y = yv(te/h+1);% returns y at t=te; invariant: tv(te/h+1)==te
y1=y-ye;
[tv,yv] = euler2_(h,0,te,0,x2); %x1 is the dy0, that is unknown and
assume It is x2)
y = yv(te/h+1); % returns y at t=te
y2=y-ye;
while abs(x2-x1)>tol
disp(sprintf('(%g,%g) (%g,%g)', x1, y1, x2, y2));
x3 = x2-y2*(x2-x1)/(y2-y1); % interporation value of dy0 is x3
%================
[tv,yv] = euler2_(h,0,te,0,x3);
 plot(tv,yv,'*b', 'LineWidth',2);
y3 = yv(te/h+1)-ye; %yv(te/h+1)returns y at t=te; solve equation
y=ye
%===================
x1 = x2; y1 = y2;
x2 = x3; y2 = y3;
end
x = x2;
return;
end
%=======================================================
function [tv,yv] = euler2_(h,t0,tmax,y0,dy0)
% use Euler's method to solve 2nd order ODE y''=-g+a*y^2
% return tvec and yvec sampled at t=(t0:h:tmax) as col. vectors
% y(t0)=y0, y'(t0)=dy0
global a te ye ss; % coeff of nonlinear acceleration
g = 9.8; % accel. of gravity, [m/sec^2]
y = y0; % position
dy = dy0; % velocity
tv = [t0];
yv = [y0];
for t = t0:h:tmax
y = y+dy*h; % this and following line are Euler's method
dy = dy+(-g+a*y^2)*h;% dy=vy (velocity); dvy=(-g+a*y^2), accelerator
tv = [tv; t+h];
yv = [yv; y];
end
%figure();
% plot(tv,yv,'*', 'LineWidth',2);
% text(te+.2,y,sprintf('y\047(0)=%g',dy0), 'FontSize',15);
return;
end

%==============================================================
Bai 5:-Iteractive and Matrix Method-----
clear all; clc
global Ta heso dx n T1 x
Ta=20; heso=0.1; dx=0.05; L=10; n=L/dx;
T0=zeros(n+1,1);
T0=T0+2+heso*dx^2;
T0(1)=200;T0(n+1)=40;
tam=1e-3;
T1=T0;ss=10;
while ss>1e-6
for i=2:n
T1(i)=(T0(i-1)+T0(i+1)+dx*dx*Ta*heso)/(2+heso*dx^2);
ss=max(abs(T1-T0));
end
T0=T1;
end
x=linspace(0, 10, n+1);
plot(x, T1, 'O-b'); hold on;
xlabel('x (m)', 'FontSize',20);
ylabel('T (^oC)', 'FontSize',20);
title(sprintf('Matrix and Iterative methods for BVP: T\047\047=0.01*(T-
Ta)'), 'FontSize',20);
Matrix_()
legend('Matrix methods', 'Iterative methods')
%============================
function Matrix_()
global Ta n dx heso T1 T2 x
for i=1:n+1
for j=1:n+1
if i==j
A(i, j)=2+heso*dx^2;
elseif(i==j+1)|(j==i+1)
A(i,j)=-1;
else A(i,j)=0;
end
end
end

for i=1:n
B(i)=dx*dx*Ta*heso;
end
B(1)=200; B(n+1)=40; B=B';
T2=A\B;
T12=[T1, T2]
plot(x,T2, '*-r'); xlabel('x(m)'); ylabel('T(^oC)');
end
%==================================================

%==============================================================
Bài 6: Con lac lo xo
%--------------------main function------------------------
% k=2.5; m=1.5; g=9.8;
% state0 = [-2*m*g/k, 0.0];
state0 = [0.0, 0.0];
t = 0:.1:20;
[t,state] = ode45('MassSpring', t, state0);
plot(t, state,'Linewidth', 3)
xlabel('TIME (sec)')
ylabel('STATES')
title('Mass-Spring System')
legend({'$x$ (m)', '$\dot{x}$ (m/sec)'},'interpreter','latex')

function stated = MassSpring(t, state)

% unpack the state vector
x = state(1);
xd = state(2);
% these are our constants
k = 2.5; % Newtons per metre
m = 1.5; % Kilograms
g = 9.8; % metres per second
% compute acceleration xdd
xdd = ((-k*x)/m) + g;
stated = [xd; xdd];
end

%================================================================
Bài 7: -------------Vật rơi-------------------------
clear all; clc
[t, v] = ode45(@v_falling, [0,30], 0);
plot(t,v, 'r', 'Linewidth', 3)
title("falling with air resistance")
function vp = v_falling(t,v)
% function vp = v_falling(t,v)% compute dv/dt
m = 70;
g = 9.8;
b = .25; vp = -g +(b/m)*v.^2;
end

%==================================================================
Bai 8:
a/ Chương trình tạo mạng tinh thể SC và tính số phối trí, hàm phân bố xuyên tâm:
clear all, clc
%====================create a simple cubic lattice=====
nn=2;a=1;l=nn*a; l2=l/2;
n=1;
for i=1:nn
for j=1:nn
for k=1:nn
x(n)=(i-1)*a; y(n)=(j-1)*a; z(n)=(k-1)*a; c(n)=0; n=n+1;
end
end
end
%=====================Calculate coordination number================
rmin=1.001;% distance of nearest neighbor
SPT(1,8)=0;N=nn^3;
for i=1:N
k=0;
for j=1:N
if(i~=j)
r=khoangcach(x(i), y(i), z(i), x(j), y(j), z(j), l, l2);;
if (r<rmin)
k=k+1;
end
end
end
PT(i)=k;
SPT(k)=SPT(k)+1;
end
No=1:8;
Tabl=[SPT', No'];
fprintf(' No atoms CN \n');
disp(Tabl)

%===========Calculate radial distribution function=====

SizeG=1000; g=zeros(SizeG, 1);

dg=0.01;
for i=1:SizeG
R(i)=i*dg;
end
L=dg*(SizeG-1);
if L>l
L=l;
end
for i=1:N
for j=(i+1):N
r=khoangcach(x(i), y(i), z(i), x(j), y(j), z(j), l, l2);
if(r<L) k=round(1+r/dg);
g(k)=g(k)+1.0;
end
end
end

for i=1:SizeG
g(i)=g(i)/N/(2*pi*dg*dg*i*i*dg+1.e-20)/(N/(l^3));% tinh theo cong thuc
end
plot(R, g)
%axis([0 10 0 1])
%==========================function khoangcach================
function kc=khoangcach(x1, y1, z1, x2, y2, z2, l, l2)
rx=abs(x1-x2);if(rx>l2)rx=rx-l;end
ry=abs(y1-y2);if(ry>l2)ry=ry-l;end
rz=abs(z1-z2);if(rz>l2)rz=rz-l;end
kc=sqrt(rx^2+ry^2+rz^2);
end
%==========================================================
b/ Chương trình tạo mạng tinh thể BCC và tính số phối trí, hàm phân bố xuyên tâm:
clear all, clc
%====================create a simple cubic lattice=====
nn=2;a=1;l=nn*a; l2=l/2; N=2*nn^3;
n=1;
for i=1:nn
for j=1:nn
for k=1:nn
x(n)=(i-1)*a; y(n)=(j-1)*a; z(n)=(k-1)*a; c(n)=0; n=n+1;
x(n)=(i-1)+0.5; y(n)=(j-1)+0.5; z(n)=(k-1)+0.5; c(n)=1;n=n+1;
end
end
end
%=====================Calculate coordination number================
rmin=1.001;% distance of nearest neighbor
SPT(1,20)=0;
for i=1:N
k=0;
for j=1:N
if(i~=j)
r=khoangcach(x(i), y(i), z(i), x(j), y(j), z(j), l, l2);;
if (r<rmin)
k=k+1;
end
end
end
PT(i)=k;
SPT(k)=SPT(k)+1;
end
No=1:20;
Tabl=[SPT', No'];
fprintf(' No atoms CN \n');
disp(Tabl)

%===========Calculate radial distribution function=====

SizeG=1000; g=zeros(SizeG, 1);

dg=0.01;
for i=1:SizeG
R(i)=i*dg;
end
L=dg*(SizeG-1);
if L>l
L=l;
end
for i=1:N
for j=(i+1):N
r=khoangcach(x(i), y(i), z(i), x(j), y(j), z(j), l, l2);
if(r<L) k=round(1+r/dg);
 g(k)=g(k)+1.0;
end
end
end

for i=1:SizeG
g(i)=g(i)/N/(2*pi*dg*dg*i*i*dg+1.e-20)/(N/(l^3));% tinh theo cong thuc
end
plot(R, g)

c/ Chương trình tạo mạng tinh thể FCC và tính số phối trí, hàm phân bố xuyên tâm:
clear all, clc
%====================create a simple cubic lattice
nn=2;a=1;l=nn*a; l2=l/2;N=4*nn^3;
n=1;
for i=1:nn
for j=1:nn
for k=1:nn

x(n)=(i-1)*a; y(n)=(j-1)*a; z(n)=(k-1)*a; n=n+1;c(n)=0;

x(n)=(i-1)*a+0.5; y(n)=(j-1)*a+0.5; z(n)=(k-1)*a;n=n+1; c(n)=1;
x(n)=(i-1)*a+0.5; y(n)=(j-1)*a; z(n)=(k-1)*a+0.5; n=n+1;c(n)=1;
x(n)=(i-1)*a; y(n)=(j-1)*a+0.5; z(n)=(k-1)*a+0.5; n=n+1; c(n)=1;
end
end
end
%=====================Calculate coordination number================
rmin=0.8;% distance of nearest neighbor
SPT(1,N)=0;
for i=1:N
k=0;
for j=1:N
if(i~=j)
r=khoangcach(x(i), y(i), z(i), x(j), y(j), z(j), l, l2);
if (r<rmin)
k=k+1;
end
end
end
PT(i)=k;
SPT(k)=SPT(k)+1;
end
No=1:N;
Tabl=[SPT', No'];
fprintf(' No atoms CN \n');
disp(Tabl)

%===========Calculate radial distribution function=====

SizeG=1000; g=zeros(SizeG, 1);

dg=0.01;
for i=1:SizeG
R(i)=i*dg;
end
L=dg*(SizeG-1);
if L>l
L=l;
end
for i=1:N
for j=(i+1):N
r=khoangcach(x(i), y(i), z(i), x(j), y(j), z(j), l, l2);
if(r<L) k=round(1+r/dg);
g(k)=g(k)+1.0;
end
end
end

for i=1:SizeG
g(i)=g(i)/N/(2*pi*dg*dg*i*i*dg+1.e-20)/(N/(l^3));% tinh theo cong thuc
end
plot(R, g)

%==========================function khoangcach================
function kc=khoangcach(x1, y1, z1, x2, y2, z2, l, l2)
rx=abs(x1-x2);if(rx>l2)rx=rx-l;end
ry=abs(y1-y2);if(ry>l2)ry=ry-l;end
rz=abs(z1-z2);if(rz>l2)rz=rz-l;end
kc=sqrt(rx^2+ry^2+rz^2);
end
%=============================================================

Bài 9: Xây dựng chương trình mô phỏng động lực học phân tử có 1000 nguyên tử, sử dụng thế Lennard-Jones
với  = 1; L = 10*a (a là hằng số mạng, lấy a=2).
------------------------------%Chương trình chính: Main.m -------------------------------------------------------1
% main.m
clear all; clc
initial; M=300;
for i1=1:M
newcoordinate;
if rem(i1,10) == 1
energy;
fprintf('%10.3f %10.3f %10.3f\n', K, U, K+U);
end
updatecoordinate;
end
----------------initial.m-----------------------------2
%initial.m
n0=10; a=2; l=n0*a; l2=l/2; h=0.001; n=0;
for i=1:n0
 for j=1:n0
for k=1:n0
n=n+1; x1(n)=(i-1)*a;
y1(n)=(j-1)*a; z1(n)=(k-1)*a; x2(n)=x1(n)+0.0001*(-1)^n;
y2(n)=y1(n)+0.0001*(-1)^n; z2(n)=z1(n)+0.0001*(-1)^n;
end
end
end
--------------- energy.m------------------------------3
% energy.m
U=0; K=0;
for i=1:n-1
for j=i+1:n
rx=x2(i)-x2(j); ry=y2(i)-y2(j); rz=z2(i)-z2(j);
rx=dd(rx, l, l2); ry=dd(ry, l, l2); rz=dd(rz, l, l2);
r=sqrt(rx^2+ry^2+rz^2); U=U+4*pot(r, r0);
end
end
for i=1:n
rx=x(i)-x1(i); ry=y(i)-y1(i); rz=z(i)-z1(i);
rx=dd(rx, l, l2); ry=dd(ry, l, l2); rz=dd(rz, l, l2);
K=K+(rx^2+ry^2+rz^2)/r0/r0;
end
K=K/8/h/h;
--------------------- pot.m------------------------------4
function u=pot(r, r0)
r=r0/r; u=r^12-r^6;
end

-------------------force.m--------------------------------5

% force.m
function f=force(r, r0)
r=r0/r; f=(-12*r^13+6*r^7)/r0;
end
-------------------dd.m------------------------------------6
% dd.m
function d2=dd(d1, l, l2)
d2=d1;
if d1> l2
d2=d1-l ;
end
if d1<-l2
d2=d1+l ;
end
end

---------------------- updatecoordinate.m---------------------7
%updatecoordinate.m
x1=x2; x2=x; y1=y2; y2=y; z1=z2; z2=z;

-------------------- newcoordinate.m------------------------8
% newcoordinate.m
fx=zeros(1,n); fy=fx; fz=fx;r0=1;
for i=1:n-1
for j=i+1:n
rx=x2(i)-x2(j); ry=y2(i)-y2(j); rz=z2(i)-z2(j);
rx=dd(rx, l, l2); ry=dd(ry, l, l2); rz=dd(rz, l, l2);
r=sqrt(rx^2+ry^2+rz^2); f=force(r, r0)/r;
fx(i)=fx(i)-rx*f; fx(j)=fx(j)+rx*f;
fy(i)=fy(i)-ry*f; fy(j)=fy(j)+ry*f;
fz(i)=fz(i)-rz*f; fz(j)=fz(j)+rz*f;
end
end

for i=1:n
rx=x2(i)-x1(i); ry=y2(i)-y1(i); rz=z2(i)-z1(i);
rx=dd(rx, l, l2); ry=dd(ry, l, l2); rz=dd(rz, l, l2);
x(i)=x2(i)+rx+4*fx(i)*h*h*r0;
if x(i)< 0 x(i)=x(i)+l; end
if x(i)> l x(i)=x(i)-l; end
y(i)=y2(i)+ry+4*fy(i)*h*h*r0;
if y(i)< 0 y(i)=y(i)+l; end
if y(i)> l y(i)=y(i)-l; end
z(i)=z2(i)+rz+4*fz(i)*h*h*r0;
if z(i)< 0 z(i)=z(i)+l; end
if z(i)> l z(i)=z(i)-l; end
end

%==========================Bài 10=============================
Bai 10: Tính toán phân bố nhiệt trên một tấm kim loại với một nguồn tạo nhiệt ở giữa của tấm phẳng như
hình phía dưới. Giả sử tấm phẳng bằng hợp kim đồng có độ dẫn nhiệt là C=0.3 W.K/cm. Nguồn tạo nhiệt có công
suất là P=100W/cm3. Chiều dày của tấm phẳng rất nhỏ so với các chiều khác. Bỏ qua dòng nhiệt dọc theo hướng
z (chiều dày của tấm). Tấm phẳng hình vuông có cạnh là L2=10 cm, Nguồn tạo nhiệt hình vuông có cạnh là L1=3
cm
Phân bố nhiệt độ tuân theo phương trình sau

  2T ( x , y )  2T ( x , y ) P
    L1
L2
 x 2 y 2 C
 T ( x , y )  3 0; C tại biên của tấm phẳng
o

%--------------------------------------------------------------
clear all
P=200; C=0.3;
Q1=P/C;L2=10; L1=3;
Lx=L2;Ly=L2;dx=0.5;dy=dx;
nx=round(Lx/dx);
ny=round(Ly/dy);
n0=round((L2-L1)/2/dx); n1=round((L2+L1)/2/dx);
Q=zeros(ny+1, nx+1);
Q(n0:n1,n0:n1)=Q1;
%w=0.523647138;
for j=1:nx+1
for i=1:ny+1
T(i,j)=30;
TNew(i,j)=T(i,j);
end
end
%repeat
repeat=1;
while (repeat>1.e-9)
for j=2:nx
for i=2:ny
dT=(T(i+1,j)+T(i,j+1)+T(i-1,j)+T(i,j-1)-4*T(i,j)+dx*dx*Q(i,j))/4;
TNew(i,j)=T(i,j)+dT;
end
end
repeat= max(max(abs(T-TNew)));
for j=2:nx
for i=2:ny
T(i,j)=TNew(i,j);
end
end
end
[x,y]=meshgrid(0:dx:Lx, 0:dy:Ly);
subplot(2,1,1); surface(x,y,T);
axis([0 Lx 0 Ly 0 max(max(T))])
subplot(2,1,2); contour(T, 10);
[C,h] = contour(T, 10);
set(h,'ShowText','on','TextStep',get(h,'LevelStep')*1);
f=fopen('secondoder.txt','w')
fprintf(f,' ');
for i=1:nx+1
fprintf(f,' x=%3.2f',(i-1)*dx);
end
fprintf(f,'\n')
for j=1:ny+1
fprintf(f,'y=%3.2f',(j-1)*dy);
for i=1:nx+1
fprintf(f,'%11.2f',T(j,i));
 end
fprintf(f,'\n');
end
fclose(f);

%================================Bai 11==============
Bài 11: Xây dựng chương trình tạo các số ngẫu nhiên bằng phương pháp khác nhau.
a/ Phương pháp chuyển đổi, với hàm phân bố như sau:
𝑐. 𝑥 1≤𝑥≤2
𝑃 (𝑥) =
0 𝑥 < 1 ℎ𝑜ặ𝑐 𝑥 > 2

b/ Phương pháp loại trừ với hàm phân bố như sau

𝑐. 𝑠𝑖𝑛𝑥 0.3 ≤ 𝑥 ≤ 2
𝑃 (𝑥) =
0 𝑥 < 0.3 ℎ𝑜ặ𝑐 𝑥 > 2

-------------------a/Phương pháp chuyển đổi -----------------

clear all; clc
m=1000;
y=rand(m, 1);
x=((2.^6-1).*y+1).^(1/6); %he so chuan hoa la C=(6/2^6-1)=>
%integral(C.x'^6, 1, x)=y=>x=((2.^6-1).*y+1).^(1/6)
title('Distribution of random points according to C*x^5'); hold on;
figure(1)
hist(x)
figure(2)
x1=1:0.1:2;y1=x1.^5;
plot(x1, y1);
figure(3)
index=1:m;
plot(x, index, 'b.');
xlim([1 2]);

--------------------------Phương pháp loại trừ---------------

clear all, clc

n=10000;
a=0.3, b=2;
c=1/(cos(a)-cos(b));% he so chuan hoa
cmax=c; % camx=max(c*sin(x))=c
x_ace = zeros(1,n); % save accept variables
fx_ace = zeros(1,n);% save value of function c*sin(x)
irv = 0;
while irv<n
y=a+rand*(b-a);
u=rand;
fx=c*sin(y);
 if u <= fx/cmax;
irv = irv+1;
x_ace(irv) = y; % save in x
fx_ace(irv) = fx;
end
end
figure(1); hist(x_ace); title('pdf=c*sin(x)');
axis([a, b, 0, inf])
figure(2); plot(x_ace, fx_ace, 'o')
title('fx=c.sin(x)');
axis([a, b, 0, 1])

%===============================Bai 12=========
Bài 12: Tính tích phân sau bằng phương pháp monte-carlo
.
𝐼= (sin (𝑥) + 𝑥 )𝑑𝑥
.

%------------------loi giai----------------------------

clear all; clc

%f=sin(x)+x^3;
a=0.1; b=1.2;
I=0; m=5000000;
x=a+(b-a)*rand(m,1);
f=@(x) x.^3+sin(x);
I=(b-a)*sum(f(x))/m;
format long
fprintf('ket qua la:%.4f\n\n', I)

%----------------------------------------------------------------------------------------------------

%===================Bài 13====================================

Bài 13: Xâydựng mô hình Izing hai chiều kích thước 10x10 với điều kiện biên tuần hoàn và không
có từ trường ngoài. Vẽ sự phụ thuộc moment từ M, nhiệt dung C V và độ cảm từ  vào nhiệt độ.

---------------------------Giải-----------------------------
clear all
clc
dT=0.05;
Temp0=1;
Temp=Temp0;
Tmax=4;
MM=[];k=0;CV=[];CHI=[];
while(Temp<Tmax)
k=k+1;
Temp(k)=Temp0+dT*k;
N_exp=10000;
n=10;
N_step=n*n;m=round(rand(n))*2-1;
M=sum(sum(m));
Emed=0; Emed2=0;Mmed=0; Mmed2=0;
c0=1.0/Temp(k);
for i=1:N_exp
Mg=sum(sum(m));
Ener=0;
for j=1:N_step
i1=1+round(rand*(n-1));
j1=1+round(rand*(n-1));

spin=m(i1,j1);
L=i1-1; if L<1, L=n;end
R=i1+1; if R>n, R=1; end
B=j1-1; if B<1, B=n;end
T=j1+1; if T>n, T=1; end
DE=2*m(i1,j1)*( m(L,j1)+m(R,j1)+m(i1,B)+m(i1,T) );

if DE>0
p=exp(-DE*c0);
else
p=1;
end

if rand<p
m(i1,j1)=-spin;
end
Mg=Mg-2*spin;
Ener=Ener+DE;
end
% Mg=sum(sum(m));
Ener=Ener/2/N_step;
Mg=Mg/N_step;
M=M+abs(Mg);
Emed=Emed+Ener;
Emed2=Emed2+Ener*Ener;
Mmed=Mmed+Mg;
Mmed2=Mmed2+Mg*Mg;

end
MM=[MM, M/N_exp];
Emed=Emed/N_exp;
Emed2=Emed2/N_exp;
Mmed=Mmed/N_exp;
Mmed2=Mmed2/N_exp;
CV=[CV, (Emed2-Emed*Emed)/Temp(k)/Temp(k)];
CHI=[CHI, (Mmed2-Mmed*Mmed)/Temp(k)];
end

figure, plot(Temp, MM, 'o')

figure, plot(Temp, CV, 'o')
figure, plot(Temp, CHI, 'o')

%===========================Bài 14=======================================
14. Số liệu đo thực nghiệm của hai đại lượng x và y như sau:
x 2.00 2.50 3.00 3.50 4.00 4.50 5.00

y -0.09 -0.82 -1.43 -1.93 -2.37 -2.76 -3.11

Biết x và y liên hệ với nhau bởi phương trình:

y  a ln( x )  b
Tính hệ số a, b bằng phương pháp bình phương tối thiểu
%-----------------Linear least squares Method via function polyfit-----------------------------------------
x=[2.00 2.50 3.00 3.50 4.00 4.50 5.00];
y=[-0.09 -0.82 -1.43 -1.93 -2.37 -2.76 -3.11];
X=log(x);
heso_ab=polyfit(X, y, 1);
a=heso_ab(1); b=heso_ab(2);
f=@(x1) a*log(x1)+b;
x1=2:0.01:5;
y1=f(x1);
plot(x, y,'o', x1, y1,'-b', 'linewidth', 2)
title('Function Fitting');
legend('experimental data','Function Fitting');
axis([1.8,5.1, -4, 0.1])

%=======================Bài 15====================================
15. Tốc độ phát triển của vi khuẩn k (1/ngày) phụ thuộc vào nồng độ ôxy c (mg/L) theo hàm sau:

k m ax c 2
k 
cs  c 2
Biết số liệu thực nghiệm đo sự phụ thuộc của k vào c được cho trong bảng sau
c 0.5 0.8 1.5 2.5 4
k 1.1 2.4 5.3 7.6 8.9
Tìm các hệ số kmax và cs, tính giá trị của k tại c=3 mg/L
%==============================giải==========================================
clc, clear all
c=[0.5 0.8 1.5 2.5 4];
k=[1.1 2.4 5.3 7.6 8.9];
x=1./c.^2;
y=1./k;
heso_ab=polyfit(x, y, 1);
a=heso_ab(1); b=heso_ab(2); %a=cs/kmax; b=1/kmax;
kmax=1/b; cs=kmax*a;
k1=@(c1) kmax*c1.^2./(cs+c1.^2);
x1=0.5:0.01:4;
y1=k1(x1);
plot(c, k,'o', x1, y1,'-b', 'linewidth', 2)
title('Function Fitting');
legend('experimental data','Function Fitting');
axis([0.2,4.3, 0, 9]);
fprintf('cs=%8.3f \nkmax=%8.3f \n\n', cs, kmax);
%======================================