- #1
evinda
Gold Member
MHB
- 3,836
- 0
Hello! (Smile)
Consider the initial value problem
$$\left\{\begin{matrix}
y'(t)=f(t,y(t)) &, a \leq t \leq b \\
y(a)=y_0&
\end{matrix}\right. (1)$$
I want to write a program that implements the following numerical method to solve $(1)$
$\left\{\begin{matrix}
y^{n+1}=y^n+h[\rho f(t^n,y^n)+(1-\rho)f(t^{n+1},y^{n+1})] &, n=0,1, \dots, N-1 \\
y^0=y_0 &
\end{matrix}\right.$for a uniform partition of $[a,b]$ with step $h=\frac{b-a}{N}$.
For each $\rho \in [0,1)$ the method is implicit, that means that at each step we will have to solve a nonlinear equation for the computation of $y^{n+1}$. This can be done with the use of Newton's method, which is defined as follows:
Let $g(x)$ be a function and $x^{\star}$ a root of it, which we want to approach. Given an initial approximation of the root, $x_0$, we define the iterative procedure $x_{k+1}=x_k-\frac{g(x_k)}{g'(x_k)}, k=0,1,2, \dots$ that approaches our root. There are two termination criteria of the method. Either we define a maximum number of iterations $R_{max}$ or we require two consecutive approximations to differ less than an accuracy [m] tol [/m] that we define, i.e. we terminate the procedure when we have $|x_{k+1}-x_k|<tol$.For example, if we have $\rho=0$ then we have the implicit Euler method $y^{n+1}=y^n+hf(t^{n+1},y^{n+1})$.
In our case, $g(y^{n+1})=y^{n+1}-y^n-hf(t^{n+1},y^{n+1})$.
So in order to calculate an approximation of $y^{n+1}$, we will use Newton's method:
$$y_{k+1}^{n+1}=y_k^{n+1}- \frac{g(y_k^{n+1})}{g'(y_k^{n+1})}=y_k^{n+1}- \frac{y_k^{n+1}-y^n-hf(t^{n+1},y_k^{n+1})}{1-h \frac{\partial{f(t^{n+1},y_k^{n+1})}}{\partial{y_k^{n+1}}}}, k=1,2, \dots, R_{max}$$
given some initial approximation $y_0^{n+1}$ (let $y_0^{n+1}=y^n$).That's what I have tried so far:
Is is right so far? I don't know what to do with the initial value $y_0^{n+1}=y^n$ and in general with $y_k^{n+1}$.
We want to compute $y_1^{n+1},y_2^{n+1}, y_2^{n+1}, \dots$ for several values of $n$.
So does it suffice to consider a vector of length $N$?
Or should we have a two-dimensional array? (Worried)
Consider the initial value problem
$$\left\{\begin{matrix}
y'(t)=f(t,y(t)) &, a \leq t \leq b \\
y(a)=y_0&
\end{matrix}\right. (1)$$
I want to write a program that implements the following numerical method to solve $(1)$
$\left\{\begin{matrix}
y^{n+1}=y^n+h[\rho f(t^n,y^n)+(1-\rho)f(t^{n+1},y^{n+1})] &, n=0,1, \dots, N-1 \\
y^0=y_0 &
\end{matrix}\right.$for a uniform partition of $[a,b]$ with step $h=\frac{b-a}{N}$.
For each $\rho \in [0,1)$ the method is implicit, that means that at each step we will have to solve a nonlinear equation for the computation of $y^{n+1}$. This can be done with the use of Newton's method, which is defined as follows:
Let $g(x)$ be a function and $x^{\star}$ a root of it, which we want to approach. Given an initial approximation of the root, $x_0$, we define the iterative procedure $x_{k+1}=x_k-\frac{g(x_k)}{g'(x_k)}, k=0,1,2, \dots$ that approaches our root. There are two termination criteria of the method. Either we define a maximum number of iterations $R_{max}$ or we require two consecutive approximations to differ less than an accuracy [m] tol [/m] that we define, i.e. we terminate the procedure when we have $|x_{k+1}-x_k|<tol$.For example, if we have $\rho=0$ then we have the implicit Euler method $y^{n+1}=y^n+hf(t^{n+1},y^{n+1})$.
In our case, $g(y^{n+1})=y^{n+1}-y^n-hf(t^{n+1},y^{n+1})$.
So in order to calculate an approximation of $y^{n+1}$, we will use Newton's method:
$$y_{k+1}^{n+1}=y_k^{n+1}- \frac{g(y_k^{n+1})}{g'(y_k^{n+1})}=y_k^{n+1}- \frac{y_k^{n+1}-y^n-hf(t^{n+1},y_k^{n+1})}{1-h \frac{\partial{f(t^{n+1},y_k^{n+1})}}{\partial{y_k^{n+1}}}}, k=1,2, \dots, R_{max}$$
given some initial approximation $y_0^{n+1}$ (let $y_0^{n+1}=y^n$).That's what I have tried so far:
Code:
function [t,y]=pf(a,b,y0,N, rho)
RMAX=4;
tol=10^(-3);
h=(b-a)/N;
t=zeros(1,N+1);
y=zeros=(1,N+1);
y0=zeros(1,N+1);
t(1)=0;
y(1)=1;
for n=1:N
g=@(y(n+1)) y(n+1)-y(n)-h*(rho*f(t(n),y(n))-(1-rho)*func(t(n+1),y(n+1)));
dg=@(y(n+1)) 1-h*(1-rho)*dfunc(t(n+1));
y0=y(n);
for k=1:RMAX
y(n+1)=y(0)-g(y(0))/dg(y(0))
end
endWe want to compute $y_1^{n+1},y_2^{n+1}, y_2^{n+1}, \dots$ for several values of $n$.
So does it suffice to consider a vector of length $N$?
Or should we have a two-dimensional array? (Worried)