How to start the finite difference method on this equation?

97 Views Asked by At

I have to find the numerical solutions for the boundary value problem through matlab,


$2yy\prime\prime -(y\prime)^{2} + 4y^{2} = 0, \quad y(\pi/6) = 1/4, \quad y(\pi/2) = 1. $


It took so much time to how to start, I tried $dy/dx = u$ but failed...

Can anybody help??

1

There are 1 best solutions below

4
On

If the words "difference method" are part of your task description, it would seem that you are expected to directly insert $$ y'(t)=\frac{y(t+h)-y(t-h)}{2h},\qquad y''=\frac{y(t+h)-2y(t)+y(t-h)}{h^2} $$ to get $$ 2y(t)(y(t+h)-2y(t)+y(t-h))-\frac14(y(t+h)-y(t-h))^2+4h^2y(t)^2=0 $$ which is a quadratic equation for $y(t+h)$.


Otherwise

function dz=prime(t,z)
     y=z(1), dy=z(2)
     dz = [ dy; (dy^2-4*y^2)/(2*y) ]

is a variant of the function to put into the BVP solver.


If you want to go Newton on your equation, compute the directional derivative in direction $v$ as $$ 2yv''+2y''v-2y'v'+8yv=0 $$

This works in Scilab, a close cousin of Matlab

function dz=prime(t,z)
     y=z(1), dy=z(2), v=z(3), dv=z(4)
     d2y = (dy^2-4*y^2)/(2*y);
     d2v = (dy*dv-4*y*v-d2y*v)/(y)
     dz = [ dy; d2y; dv; d2v ]
endfunction

Pi=4*atan(1)
t0 = Pi/6
y0 = 1/4
tf = Pi/2
yf = 1
a=1
dy0=1
while abs(a)>1e-14
    z0 = [y0; dy0; 0; 1]
    T = [tf]
    Z = ode(z0, t0, T, prime )
    n=length(T); n=n(1)
    a = (yf-Z(1,n))/Z(3,n)
    dy0=dy0+a
    disp([dy0,a])
end

T = linspace(t0,tf,301)
Z = ode(z0, t0, T, prime )
clf()
plot(T,Z(1,:),T,sin(T).^2)