MHB Implementation of numerical method

AI Thread Summary
The discussion focuses on implementing a numerical method to solve an initial value problem using a specific iterative approach. The method involves updating values of y using a formula that requires solving a nonlinear equation at each step, which is facilitated by Newton's method. The user is seeking clarification on how to handle initial values and whether a one-dimensional or two-dimensional array is necessary for storing results. They also share results from their function implementation, expressing confusion about the correctness of the outputs. Overall, the thread highlights the challenges of coding numerical methods and the need for proper data structure management.
evinda
Gold Member
MHB
Messages
3,741
Reaction score
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:
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
end
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)
 
Mathematics news on Phys.org
I have written a function and calling it with $a=0, b=1, N=5, y0=1, \rho=0$ and I got the following results:
[m]
t =

0 0.2000 0 0 0 0y =

1 1 0 0 0 0
1 1 0 0 0 0

[/m]

Then, changing the function, I got the following:[m]
t =

0
0.2000
0.4000
0.6000
0.8000
1.0000y =

1.0000
0.3529
15.1119
647.0365
228.3700
-133.2945

[/m]Could you tell me one of them is the right result? (Thinking)
 
Seemingly by some mathematical coincidence, a hexagon of sides 2,2,7,7, 11, and 11 can be inscribed in a circle of radius 7. The other day I saw a math problem on line, which they said came from a Polish Olympiad, where you compute the length x of the 3rd side which is the same as the radius, so that the sides of length 2,x, and 11 are inscribed on the arc of a semi-circle. The law of cosines applied twice gives the answer for x of exactly 7, but the arithmetic is so complex that the...
Back
Top