Bài báo cáo môn Phương pháp tính

ĐẠI HỌC QUỐC GIA THAØNH PHỐ HỒ CHÍ MINH

TRÖÔØNG ÑAÏI HOÏC BAÙCH KHOA

KHOA KHOA HOÏC ÖÙNG DUÏNG

Boä moân : Phöông phaùp tính



Baøi baùo caùo

GVHD : Ths Leâ Thaùi Thanh

Lôùp : XD09BXD1

Nhoùm 7

Nguyeãn Quang Thaéng

Tröông Khaùnh Duy

Toâ Döông Nhaät Ñaêng

Traàn Phöôùc Vónh

Voõ Trung Hieáu

80902548

80900428

Nguyeãn Anh Döông

Buøi Danh Toaïi

80900472

80902845

Baøi 1 :

Vieát chöông trình xaây döng ña thöùc noäi suy Newton (SGK, trang 72-75). Aùp duïng ñoái

vôùi baûng soá :

X

Y

6.0

7.0

2.3

8.0

9.2

10.5

1.4

11.3

0.9

11.6

0.7

12.0

0.6

12.6

0.5

13.0

0.4

2.25

1.95

Veõ ñoà thò ña thöùc thu ñöôïc.

Baøi laøm

Söû duïng PHP :

Code :

<h1> Noi suy Newton - Design by Khanh Duy</h1>

<?php

function noisuynewton($xcankim) {

//$xcantim = 13.5;

//so phan tu

$a[1][0] = 10;

$x[1] = 6;

$x[2] = 7;

$x[3] = 8;

$a[1][1] = 2.25;

$a[1][2] = 2.3;

$a[1][3] = 2.25;

$x[4] = 9.2; $a[1][4] = 1.95;

$x[5] = 10.5; $a[1][5] = 1.4;

$x[6] = 11.3; $a[1][6] = 0.9;

$x[7] = 11.6; $a[1][7] = 0.7;

$x[8] = 12; $a[1][8] = 0.6;

$x[9] = 12.6; $a[1][9] = 0.5;

$x[10] = 13; $a[1][10] = 0.4;

// khoi tao cot 0 - so phan tu

for ($i = 2 ; $i <= $a[1][0] ; $i++) {

$a[$i][0] = $a[1][0] - $i + 1 ;

}

// sinh mang 2 chieu

for ($i = 2; $i <= $a[1][0] ; $i++) {

for ($j = 1; $j <= $a[$i][0] ; $j++) {

$a[$i][$j] = ($a[$i-1][$j+1] - $a[$i-1][$j]) / ($x[$i] -

$x[1]);

//debug

//echo $a[$i][$j] . "<br>";

}

// xuat ra man hinh

// tinh ket qua cuoi cung voi xcankim

$s = 0;

for ($i=1; $i <= $a[1][0] ; $i++) {

$stam = 1;

if ($i >1) {

for ($j=1 ; $j <= $i-1; $j++) {

$stam = $stam * ($xcankim - $x[$j]);

//debug

//echo $xcankim . " ";

//echo $x[$j];

//echo "<br>";

}

// debug

//echo "<br>";

//echo $stam . " ";

//echo $a[$i][1] . " ";

//echo "<br>";

$s = $s + $a[$i][1] * $stam;

}

return $s;

}

?>

<div id="dothi" style="position:relative;border:1px solid black; width:1000px;

height:1000px;">

<?php

for ($i = 0; $i <= 500 ; $i = $i+0.05) {

$top = 1000 - noisuynewton($i);

$left = $i;

echo "<div class='diem'

style='position:absolute;width:10px;height:10px;background-color:black;top:$top;left:$left;'>

</div>" ;

} ?>

</div>

Keát quaû :

Baøi 2 :

Cho phöông trình vi phaân caáp hai :

Söû duïng phöông phaùp Runge – Kutta baäc 4, haõy xaáp xæ y(x) trong [0.1] vôùi böôùc h = 0.1. Veõ ñoà thò

haøm y(x)

Baøi laøm

Söû duïng matlab

Code

function Bai2RungeKutta

% Voi quy uoc dat y1= y(x) va y2= y'(x)

% => y'(x)= y2 va y"(x)=f la ham so theo bien (y2, y1, x)

f= input('Nhap ham so y"(y2,y1,x)= ');

f= inline(f);

% Nhap f= '(x^2+1)*y2 + y^2 - (x^2+1)^1/3'

% Doi chuoi ki tu thanh ham so

% Nhap x0= 0

x0= input('Nhap gia tri diem dau x0= ');

y10= input('Nhap dieu kien ban dau y1(x0)= '); % Nhap y(x0)= 1

y20= input('Nhap dieu kien ban dau y2(x0)= '); % Nhap y'(x0)= 0

h= input('Nhap gia tri buoc nhay h = ');

x1= input('Nhap gia tri diem cuoi x1= ');

% Nhap h= 0.1

% Nhap x1= 1

[x,y1,y2]= RungeKuta(f,x0,y10,y20,h,x1); % Goi ham con RungeKuta va xuat ra man hinh ket qua tra

ve cua ham

disp(' ');

% Bo trong 1 dong

clf;

% Xoa cac do thi da ve. Dung dong nay khi ve do thi lien tiep cho cac ham khac

nhau.

figure(2);

% Tao figure so 2

subplot(111);

plot(x, y1, 'r-');

hold on;

% Tao vung do hoa so 1 o hang 1 cot 1

% Ve do thi tron x va y voi net lien mau do

% Giu nguyen do thi vua ve

plot(x, y2, 'g--');

grid on;

% Ve do thi tron x va y voi net dut mau xanh

% Hien mat luoi

title('\it{Do thi bai toan phuong trinh vi phan cap 2}', 'Fontsize', 16); % Ghi ten do thi

xlabel('Truc x');

ylabel('Truc y');

% Ghi ten truc hoanh

% Ghi ten truc tung

legend('Duong bien dien y1(x)', 'Duong bieu dien y2(x)', 2); % Ghi chu thich net ve cho duong bieu

dien

global n

% Bien n duoc dung nhu bien toan cuc

x= x(1:n+1);

% Trich vector x tu phan tu thu I den phan tu thu n+1

% Trich vector y1 tu phan tu thu I den phan tu thu n+1

% Trich vector y2 tu phan tu thu I den phan tu thu n+1

% Tao vector hang SoLanIn co n+1 phan tu

y1= y1(1:n+1);

y2= y2(1:n+1);

SoLanIn= 1:n+1;

s=num2str(n+1);

% Chuyen bien n thanh chuoi ki tu va dat ten la s

disp(['Ket qua duoc in ra ', s, ' cap gia tri x va y.']);

disp([' ', 'So lan in Gia tri x Gia tri y1(x) Gia tri y2(x)']);

disp([SoLanIn' x' y1' y2']); % Chuyen thanh vector cot va in ra man hinh

disp(' ');

% Bo trong 1 dong

function [x,y1,y2]= RungeKuta(f,x0,y10,y20,h,x1)

global n

% Bien n duoc dung nhu bien toan cuc

% Xoa tat ca phan tu cua vector x va tao vector x moi voi phan tu dau tien

x= []; x= [x,x0];

x(1) la x0

n= (x1-x0)/h;

for j=1:n, x(j+1)= x(j)+h; end;% Tinh gia tri x(j) ke tiep

y1= []; y1=[y1,y10];

tien y1(1) la y10

% Xoa tat ca phan tu cua vetcor y1 va tao vector x moi voi phan tu dau

y2= []; y2=[y2,y20];

tien y2(1) la y20

% Xoa tat ca phan tu cua vector y2 va tao vector x moi voi phan tu dau

for j = 1:n

k11= h*y2(j);

k12= h*f(x(j),y1(j),y2(j));

k21= h*(y2(j)+0.5*k12);

k22= h*f(x(j)+h/2,y1(j)+0.5*k11,y2(j)+0.5*k12);

k31= h*(y2(j)+0.5*k22);

k32= h*f(x(j)+h/2,y1(j)+0.5*k21,y2(j)+0.5*k22);

k41= h*(y2(j)+k32);

k42= h*f(x(j)+h,y1(j)+k31,y2(j)+k32);

y1(j+1) = y1(j) + (k11+2*k21+2*k31+k41)/6;% Tinh gia tri y1(j) ke tiep - tuc la y(x(j))

y2(j+1) = y2(j) + (k12+2*k22+2*k32+k42)/6;% Tinh gia tri y2(j) ke tiep - tuc la y'(x(j))

end;

Keát quaû :

Ñoà thò

Baøi 3 :

Haøm y(x) laø nghieäm cuûa baøi toaùn bieân tuyeán tính caáp hai :

1. Söû duïng phöông phaùo sai phaân höõu haïn, tìm nghieäm gaàn ñuùng haøm y(x) trong ñoaïn

[0.100] vôùi böôùc h = 0.1

2. Vôùi keát quaû vöøa tìm ñöôïc, haõy veõ ñoà thò cuûa haøm.

Baøi laøm

Söû duïng matlab

Code

function Bai3BaiToanBien

p= input('Nhap vao ham p= ');

p= inline(p);

% Nhap p= '(1)'

% Chuyen chuoi ki tu thanh ham so

% Nhap q= 'x* cos(x)'

q= input('Nhap vao ham q= ');

q= inline(q);

r= input('Nhap vao ham r= ');

r= inline(r);

% Nhap r= '2* sin(x)'

% Nhap f= '(x-1)/(x+1)'

% Nhap a= 0

f= input('Nhap vao ham f= ');

f= inline(f);

a= input('Nhap vao gia tri diem dau a= ');

b= input('Nhap vao gia tri diem cuoi b= ');

% Nhap b= 100

h= input('Nhap vao gia tri buoc nhay h= '); % Nhap h= 0.1

ya= input('Nhap vao gia tri diem dau y(a)= '); % Nhap y(0)= 0

yb= input('Nhap vao gia tri diem cuoi y(b)= '); % Nhap y(100)=0

s= input('Nhap vao so lan in ket qua s= ');

[y,x]= BaiToanBien(p,q,r,f,a,b,h,ya,yb);%Goi ham con BaiToanBien va Dau ; cho phep ko in

ra man hinh ra ket qua tra ve

x;

% Co the bo dong nay vi dau ; cho phep ko in ra man hinh tat ca gia tri cua x

% Co the bo dong nay vi dau ; cho phep ko in ra man hinh tat ca gia tri cua y

% Bo trong 1 dong

y;

disp(' ');

figure (3);

% Tao figure so 3

clf;

% Xoa cac do thi da ve. Dung dong nay khi ve do thi lien tiep cho cac ham

khac nhau.

subplot(121);

plot(x, y, 'r-');

hold on;

% Tao vung do hoa so 1 o hang 1 cot 2

% Ve do thi tron x va y voi net lien mau do

% Giu nguyen do thi vua ve

grid on;

% Hien mat luoi

title('\sl{Do thi bai toan bien tuyen tinh cap 2}', 'Fontsize', 16); % Ghi ten do thi voi co chu 16

va in chu nghieng

xlabel('Truc x')

ylabel('Truc y')

x= x(1:s);

% Ghi ten truc hoanh

% Ghi ten truc tung

% Trich vector x tu phan tu thu I den phan tu thu s

% Trich vector y tu phan tu thu I den phan tu thu s

y= y(1:s);

SoLanIn= 1:s;

s=num2str(s);

% Tao vector hang SoLanIn co 2 phan tu

% Chuyen bien s thanh chuoi ki tu

disp(['Ket qua duoc in ra ', s, ' cap gia tri x va y.']);

disp([' ', 'So lan in Gia tri x Gia tri y']);

disp([SoLanIn' x' y']); % Chuyen thanh vector cot va in ra man hinh

subplot(122);

plot(x, y, 'b--');

hold on;

% Tao vung do hoa so 2 o hang 1 cot 2

% Ve do thi tron x va y voi net dut mau xanh

grid on;title('\it{Do thi bai toan bien tuyen tinh cap 2}', 'Fontsize', 16);

xlabel('Truc x');

ylabel('Truc y');

disp(' ');

return;

function [y,x]= BaiToanBien(p,q,r,f,a,b,h,ya,yb)

x= [];x= [x,a]; % Xoa tat ca phan tu cua vector x va tao vector x moi voi phan tu dau tien x(1)

la a

n= (b-a)/h;

for k= 1:n+1, x(k)= a+(k-1)*h; end; % Tinh gia tri x(i) ke tiep

y(1)= ya; y(n+1)= yb;

for k= 1:n-1, A(k)= r(x(k+1))-2*p(x(k+1))/h/h; end; % Tinh cac phan tu cua he so Ak

for k= 1:n-2, B(k)= p(x(k+1))/h/h+q(x(k+1))/2/h; end; % Tinh cac phan tu cua he so Bk

for k= 1:n-2, C(k)=p(x(k+2))/h/h-q(x(k+2))/2/h; end; % Tinh cac phan tu cua he so Ck

for k= 1:n-1, D(k)=f(x(k+1)); end;

% Tinh cac phan tu cua he so Dk

D(1)= D(1)-(p(x(2))/h/h-q(x(2))/2/h)*ya;

D(n-1)= D(n-1)-(p(x(n))/h/h+q(x(n))/2/h)*yb;

Y= MaTranBaDuongCheo(n-1,A,B,C,D);% Goi ham con MaTranBaDuongCheo va Dau ; cho

phep ko in ra man hinh ket qua tra ve

for k= 1:n-1, y(k+1)= Y(k); end;

function [x]= MaTranBaDuongCheo(N,a,b,c,d) % a, b, c la 3 duong cheo chinh va d la vector

bo sung cua ma tran mo rong

u(1)= a(1); v(1)= b(1); l(1)= c(1)/u(1); y(1)= d(1);

for i=2:N-1

u(i)= a(i)-l(i-1)*v(i-1);

v(i)= b(i);

l(i)= c(i)/u(i);

y(i)= d(i)-l(i-1)*y(i-1);

end;

u(N)= a(N)-l(N-1)*v(N-1);

y(N)= d(N)-l(N-1)*y(N-1);

x(N)= y(N)/u(N);

for i=N-1:-1:1

x(i)= (y(i)-v(i)*x(i+1))/u(i); % Ket qua cua Bai Toan Bien la tim tat ca cac gia tri y(x)

nghiem dung Ma Tran Ba Duong Cheo

end;

Keát quaû :

Ñoà thò :

Bài báo cáo môn Phương pháp tính - Nguyễn Quang Thắng