Advanced_Math_Final_Part

1.1. Ví dụ 1:

Ví dụ 1: Vẽ đồ thị nghiệm riêng của phương trình vi phân: \[y''-(1-y^{2})y'+y=0\] Với điều kiện ban đầu là: \(y(0)=2, y'(0)=0\) trên miền \([0,100]\)

Đặt \(y[1]=y, y[2]=y'\), ta đưa phương trình vi phân về hệ: \[y[1]'=y[2]+0.y[1]\] \[y[2]''=(1-y[1]^{2})y[2]+y[1]\] Điều kiện ban đầu \(y1=2, y2=0\)

library(deSolve)

## Warning: package 'deSolve' was built under R version 4.0.5

hamso <- function(x, y, parms){list(c(y[2], (1-y[1]^2)*y[2]-y[1]))}
dieukien <- c(y1=2, y2=0)
nghiem <- ode(y=dieukien, func= hamso, times=0:100, parms=1)
plot(nghiem, type="l", which="y1", lwd=2, ylab="Truc Y", main="Phuong trinh vi phan cap 2")

1.2. Ví dụ 2:

Vẽ đồ thị nghiệm riêng của phương trình vi phân: \[y''-(1-y^{2}y')+y=e^{x}sinx\], với điều kiện ban đầu là \(y(0)=2, y'(0)=0\) trên miền \([0,100]\)

ham <- function (x, y,parms) {list(c(y[2], (1 - y[1]^2) * y[2] - y[1] + exp(1)^x*sin(x)))}
dieukien <- c(y1 = 2, y2 = 0)
nghiem <- ode(y = dieukien, func = ham,times = 0:100, parms = 1)

## DLSODA-  At T(=R1) and step size H(=R2), the error
##       test failed repeatedly or with ABS(H) = HMIN
## In above message, R1 = 47.1239, R2 = 1.06413e-11
## 

## Warning in lsoda(y, times, func, parms, ...): repeated error test failures on a
## step, but integration was successful - singularity ?

## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go

plot(nghiem, type = "l", which = "y1",lwd = 2, ylab = "Truc y",main = "PTVP cấp 2")

2.1. Phương trình sai phân cấp 1

Xét phương trình sai phân cấp 1: y(n+1)=f(n,y(n)). Giả sử ta cần giải phương trình sai phân y(n+1)+2y(n)=0. Phương trình này có nghiệm tổng quát là \(y(n)= C(-2)^{n}\) với điều kiện ban đầu là y(2)=3 ta được nghiệm riêng là \(y(n)=3/4(-2)^{n}\). Ta sẽ vẽ đồ thị nghiệm riêng này: Vẽ đồ thị nghiệm riêng của phương trình sai phân y(n+1) + 2y(n) = 0 với điều kiện ban đầu y(2) = 3 trên miền [2, 20]

ham <- function (n, y,parms) {list(c(-2*y))}
dieukien <- c(y=3)
ketqua <- ode(y = dieukien, func = ham,times = 2:20, parms = 0, method = "iteration")
plot(ketqua, type = "l", which = "y",lwd = 2, xlab = "Truc n", ylab = "Truc y",main = "PTSP")

# Để thực hiện hiển thị các kết quả của hàm y, thực hiện lệnh:
ketqua

##    time       y
## 1     2       3
## 2     3      -6
## 3     4      12
## 4     5     -24
## 5     6      48
## 6     7     -96
## 7     8     192
## 8     9    -384
## 9    10     768
## 10   11   -1536
## 11   12    3072
## 12   13   -6144
## 13   14   12288
## 14   15  -24576
## 15   16   49152
## 16   17  -98304
## 17   18  196608
## 18   19 -393216
## 19   20  786432

Ghi chú: times = 2:20 thể hiện miền nghiệm là [2, 20], vì trong PTSP biến n thường là biến thời gian nên có quy ước chuẩn là times. Phần mềm sẽ xuất ra đồ thị (trên 1 cửa sổ khác – R Graphic) như hình trên (bạn có thể vào File -> Copy to the clipboard (CTRL+C) và paste vào word) Giải PTSP chỉ khác PTVP ở: method = “iteration”

2.2. Ví dụ bổ sung

Vẽ đồ thị nghiệm riêng của phương trình sai phân y(n+1) = (n+1)y(n) + (n+1)!n; với điều kiện ban đầu y(4) = 2 trên miền [4, 10]. Cần bổ sung thêm gói \(pracma\)

library(pracma)

## Warning: package 'pracma' was built under R version 4.0.5

## 
## Attaching package: 'pracma'

## The following object is masked from 'package:deSolve':
## 
##     rk4

ham <- function (n, y,parms) {list(c((n+1)*y+(n+1)*fact(n)))}
dieukien <- c(y = 2)
nghiem <- ode(y = dieukien, func = ham, times = 4:10, parms = 0, method = "iteration")
plot(nghiem, type = "l", which = "y",lwd = 2, xlab = "Truc n", ylab = "Truc y",main = "PTSP")