M - File
1. Tính tích phần xác định theo công thức hình thang --- P1
2. Tính gần đúng tích phân xác định theo phương pháp Simpson O = Parabol(FUN,a,b,tol) ---P1
3. Richard.m --- P2
4. Cài đặt chương trình dùng pp MONTE – CARLO tính thể tích của phần hình cầu: (x-3)^2+(y-2)^2+(z+1)^2<=24 --- P2
5. Tính thể tích elipxoid --- P2
6. tính thể tích hình paraboloid-eliptic z=45-(x/2-3)^2-(y/3-2)^2 nằm phía trên mp oxy bằng pp montecarlo --- P2
Hàm nội trú --- P2
Bài 1: tính tích phân xác định bằng phương pháp Simson và Newton_cotes19
Bài 2: Tính tích phân kép- 19
Bài 3: Tính tích phân bội 3- 19
Bài 5: Tính xấp xỉ hàm f”(x) tại x = 1,5; 2,0; 3,2 bằng phương pháp nội suy
==============
A - M – file
1. Tính tích phần xác định theo công thức hình thang
%Ham HThang.m tinh tich phan xac dinh theo
%cong thuc hinh thang
function o = HThang(FUN,a,b,tol)
if nargin == 3
RelTol = 1e-3;
AbsTol = 1e-6;
else
AbsTol = tol(1);
RelTol = tol(2);
end
n = 200;
h = (b-a)/n;
It = 0;
I = 0;
I1 = 0;
for i = 1:n
t = feval(FUN,a + i*h);
It = It + t;
end
It = h*(It + (feval(FUN,a) + feval(FUN,b))/2);
for i = 1:n
y = feval(FUN,a + i*h);
y1 = feval(FUN,a + h*(i-1));
I = I + h*(y1 + y)/2;
o = I;
if abs(It - I) < RelTol
break;
end
y = feval(FUN,a + h/2 + i*h);
y1 = feval(FUN,a + h/2 + h*(i-1));
I1 = I1 + h*(y + y1)/2;
o = I1;
if abs(It - I1) < AbsTol
break;
end
end
2. Tính gần đúng tích phân xác định theo phương pháp Simpson O = Parabol(FUN,a,b,tol)
function o = Parabol(FUN,a,b,tol)
if nargin == 3
RelTol = 1e-3;
AbsTol = 1e-6;
else
AbsTol = tol(1);
RelTol = tol(2);
end
n = 200;
h = (b-a)/2*n;
It1 = 0;
It2 = 0;
I = 0;
I1 = 0;
for i = 1:2*n
t = feval(FUN,a + i*h);
if mod(i,2) == 0
It1 = It1 + t;
else
It2 = It2 + t;
end
end
It = h*(4*It2 +2*It1 + (feval(FUN,a) + feval(FUN,b)))/3;
disp(It);
for i = 1:2*n
y = feval(FUN,a + i*h);
y1 = feval(FUN,a + h*(i-1));
I = I + h*(y1 + y)/2;
o = I;
if abs(It - I) < RelTol
break;
end
y = feval(FUN,a + h/2 + i*h);
y1 = feval(FUN,a + h/2 + h*(i-1));
I1 = I1 + h*(y + y1)/2;
o = I1;
if abs(It - I1) < AbsTol
break;
end
end
3. Richard.m
%Ham Richar.m tinh tich phan bang phuong phap ngoai suy Richardson
function Tp = Richard(FUN,a,b,Tol)
if nargin == 3
Tol(1) =1e-3;Tol(2) = 1e-6;
else if length(Tol) == 1
Tol(2) = 1e-6;
h = (b-a)/20;
x = [a+h:h:b-h];
y = feval(FUN,x);
s = sum(y) + sum(feval(FUN,[a b]))/2;
I1 = h*s;
clear x y;
h =h/2;
x = [a+h:2*h:b-h];
y = feval(FUN,x);
s = s + sum(y);
I2 = h*s;
Tp = (4*I2-I1)/3;
Tp1 = realmax;
while abs(Tp - Tp1)>Tol(2) || abs((Tp - Tp1)/Tp) > Tol (1)
clear x y;
h = h/2;
x = [a+h:2*h:b-h];
y = feval(FUN,x);
s = s + sum(y);
I1 = I2;
Tp1 = Tp2;
I2 = h*s;
Tp = (h*I2 - I1)/3;
end
4. Cài đặt chương trình dùng pp MONTE – CARLO tính thể tích của phần hình cầu: (x-3)^2+(y-2)^2+(z+1)^2<=24
% tinh the tich cua elipxoid bang montocarlo
function V=monte()
n=input('nhap so diem ngau nhien');
m=0;
for i=1:n
x=18*rand(1);
y=24*rand(1);
z=10*rand(1);
if(z<=sqrt(30-(x/2-3)^2-(y/3-2)^2)+4)
m=m+1;
end
end
V=18*24*10*m/n;
5. Tính thể tích elipxoid
clear;
N=input('Give the ranom points');
M=0;
for k=1:N
x=18*rand(1);
y=24*rand(1);
z=10*rand(1);
if(x/2-3)^2 + (y/3-2)^2 + (z -4)^2<=30
M=M+1;
end
end
V=18*24*10*M/N;
fprint('the tich:%f',V);
6. tính thể tích hình paraboloid-eliptic z=45-(x/2-3)^2-(y/3-2)^2 nằm phía trên mp oxy bằng phương pháp montecarlo
%tinh the tich cua hinh cau theo monte
function V=monte_carlo()
n=input('nhap so diem ngau nhien');
m=0;
for i=1:n
x=10*rand(1)-2;
y=10*rand(1)-3;
z=4*rand(1);
if(z<=sqrt(24-(x-3)^2-(y-2)^2)-1)
m=m+1;
end
end
V=10*10*4*m/n;
B – Hàm nội trú
Bài 1: tính tích phân xác định bằng phương pháp Simson và Newton_cotes
>> F = inline('15*log10(exp(x)+cos(x))-2*(x.^2).*tan(1./x)');
>> I = quad8(F,2,5,[1e-4 1e-8])
>> F = inline('5*log10(cos(x))-2*x.*atan(1./x)');
>> I = quad8(F,1,10,[1e-4 1e-8])
>> F = inline('3*log(abs(sin(x)))-2*x.*exp(3*x)');
>> I = quad8(F,1,5,[1e-4 1e-8])
Bài 2: Tính tích phân kép
>> F = inline('y.*cos(x)-2*x.*tan(y./x)','x','y');
>> I = dblquad(F,-pi,pi,2,5)
>> F = inline('exp(y).*sin(x)-2*x.*atan(y)','x','y');
>> I = dblquad(F,0,pi,1,2)
Bài 3: Tính tích phân bội 3
>> F = inline('x*(y^2)*sin(z)+(x^3)*exp(y)*tan(z)','x','y','z');
>> I=triplequad(F,0,5,-1,2,-pi/3,pi/3)
Bài 4: tương tự bài 3
Bài 5: Tính xấp xỉ hàm f”(x) tại x = 1,5; 2,0; 3,2 bằng phương pháp nội suy
>> y=[-3.05001;-2.90250;-1.00020;2.75754;5.37002;4.83752];
>> x=[0 0.5 1 1.5 2 2.5];
>> p=vander(x)\y;// tính hệ số của đa thức
>> f1=polyder(p');// tính đạo hàm của đa thức
>> f2=polyder(f1);
>> X=[1.5 2 3.2];
>> Z=polyval(f2,X) tính giá trị hàm
Z =
-4.9147 -14.2810 77.7692
- 19
Bạn đang đọc truyện trên: Truyen2U.Pro