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)
 
Insights auto threads is broken atm, so I'm manually creating these for new Insight articles. In Dirac’s Principles of Quantum Mechanics published in 1930 he introduced a “convenient notation” he referred to as a “delta function” which he treated as a continuum analog to the discrete Kronecker delta. The Kronecker delta is simply the indexed components of the identity operator in matrix algebra Source: https://www.physicsforums.com/insights/what-exactly-is-diracs-delta-function/ by...
Suppose ,instead of the usual x,y coordinate system with an I basis vector along the x -axis and a corresponding j basis vector along the y-axis we instead have a different pair of basis vectors ,call them e and f along their respective axes. I have seen that this is an important subject in maths My question is what physical applications does such a model apply to? I am asking here because I have devoted quite a lot of time in the past to understanding convectors and the dual...
Thread 'Imaginary Pythagoras'
I posted this in the Lame Math thread, but it's got me thinking. Is there any validity to this? Or is it really just a mathematical trick? Naively, I see that i2 + plus 12 does equal zero2. But does this have a meaning? I know one can treat the imaginary number line as just another axis like the reals, but does that mean this does represent a triangle in the complex plane with a hypotenuse of length zero? Ibix offered a rendering of the diagram using what I assume is matrix* notation...
Back
Top