User:Fuse809

From Wikipedia, the free encyclopedia
Jump to navigation Jump to search
Noia 64 apps karm.svg This user has been on Wikipedia for 7 years, 9 months and 3 days.
WP:RS This user believes in using Reliable Sources.
WikEd logo39x40 animated.gif
This user edits with the gadget wikEd.
About me
OEDThis user uses Oxford spelling.
INNThis user favours the use of International Nonproprietary Names (INNs), as designated by the World Health Organization, over all over generic names, in all circumstances.
Software
Tux.svg
Archlinux-icon-crystal-64.svgOpenSUSE Geeko button bling7.svg
This user presently contributes using theLinux distributions: Arch, NixOS and openSUSE Tumbleweed. Although, he has used other distributions in the past to contribute to Wikipedia.
Gnome footprint logo.This user contributes with GNOME.
KDEThis user contributes with KDE desktop environment.
An open box This user contributes using i3.
Firefox Logo, 2017.svgThis user contributes using Firefox.
The GIMP icon - gnome.svgThis user contributes using GIMP.
Inkscape Logo.svgThis user contributes using Inkscape.

My chemical structures[edit]

I create structural images to satisfy everyone, not just me. I try to use the structures provided by ChemSpider or PubChem just to ensure reliability. The software I use for chemical structure drawing is (beware User:Vaccinationist has been editing my files, contrary to my wishes, so some with the filenames described bellow may no longer be mine):

  • MarvinSketch/ChemSpider/PubChem (for the production of MOL files), ChemSketch (to create WMF files in ACS style) and then Scribus to convert the WMF to a SVG file (these are the images with 2DACS tags).
  • MarvinSketch (with settings according to WP:Chem/Structures) to create SVG files (2DCSD tag). Sometimes I do these images wrong the first time in that I do not crop them as well as I can and when I upload a better scaled cropped version I give it the tag 2DCSDS. Structures drawn in MarvinSketch often have whitespace in their atom labels that are due to MarvinSketch itself, e.g., see this forum topic I started about these spaces. The best solution for removing them is to use a text editor (e.g., Atom or Sublime Text) and search for xml:space="preserve" and replace it with nothing (in other words set the text editor up such that it deletes all these pieces of XML code). After this I use Inkscape to crop them. Sometimes I accidentally draw the opioids in a less conventional way, so when I fix this I give them the tag 2DCSDT.
  • MarvinSketch (to create SDF files) → OpenBabel (SDF->PDB; also adds hydrogens with this tool) → QuteMol (PDB->GIF). (3Dan tag). Unfortunately these images do not display double bonds in a way distinct from single bonds and hence I'm now favouring the 3DanJ images.
  • I also create 3D animated GIF files of drug molecules using Jmol (to create usually about 60 PNG images spaced 3 degrees apart) and GIMP to animate them into gif files. (3DanJ tag). This is the code I use to create these images:
background white;##
name = "./images/Work0000.png";
nFrames = 60;
nDegrees =6 ;
thisFrame = 0;
width = 539;
height = 260;

message loop;
thisFrame = thisFrame + 1;
fileName = name.replace("0000","" + ("0000" + thisFrame)[-4][0]);
rotate y @ndegrees;
frame next; # only use this if you have a multiframe file.
refresh;
write image @width @height @fileName;
if (thisFrame < nFrames);goto loop;endif;
  • I create 3D animated pictures of biomolecules, most often joined with their respective ligand(s) using the PDB files found here, Jmol and then PhotoScape or GIMP to convert these PNG images to GIF format. I usually colour shapely, unless doing so will make the ligand(s) difficult to see. I have now decided to use the name PDBID3DanJ for these images; for instance, https://en.wikipedia.org/wiki/File:3NT13DanJ.gif for this protein-ligand complex 3NT1.
Structures (Hidden as one of the animated images is insanely high resolution)
Image type Example's name Image example
2DACS Piritramide
Piritramide2DACSC.svg
2DCSD Dextropropoxyphene
Dextropropoxyphene2DCSD.svg
3Dan Fluoxetine
Fluoxetine3Dan3.gif
3DanJ Piritramide
Piritramide3DanJ2.gif
Biomolecule High resolution structure of naproxen:COX-2 complex. (PDB ID:3NT1)
3NT13DanJ2.gif

I upload all these images to commons now, as per the request by Leyo on my talk page. For a comprehensive list of images I've uploaded see Commons:Special:ListFiles/Fuse809.

Proof that a projectile will travel in a parabola[edit]

Note: this is using some basic energy conservation laws of physics.
Assumptions: acceleration due to gravity () does not vary substantially with position along the trajectory. Air resistance is negligible and wind is non-existent.

With the Hamiltonian:

Figure 2: William Rowan Hamilton (1805-1865), the Irish physicist and mathematician after whom the Hamiltonian is named.

also known as the total energy of the system. Where T is kinetic energy (); V is potential energy (). where: and v is the velocity.

Re-expressing v as a function of y[edit]

Basic energy conservation laws state that the Hamiltonian should be static, assuming a lack of drag. Therefore, H at one time should equal H at another time. So let us take t=0 where y=0 and v=v0 and equate it to the variable form of the Hamiltonian. That is:

Expressing y in terms of x[edit]

Let, y be a function of x. That is let:

using the chain rule we get:

Where:

It follows with the substitution of this expression into our first expression for v2 we get:

Equating results for v[edit]

Equating this expression with what we got from the law of conservation of energy we get:

As there is no forces in the x direction we can safely assume that is constant, letting this constant be we can rearrange this equation to get:

with some rearrangement we get:

Identifying the ratio of the y and x velocities as where is the launch angle of the projectile we get:

First order differential equation for y

Taking the square root of both sides gives:

Solving our first-order differential equation for y[edit]

Dividing through by the RHS and integrating with respect to x gives:

Integrating by substitution with: , ,

Solving this integral gives:

Squaring both sides and rearranging, gives:

Back-substituting again gives:

Then with some rearranging:


Our solution for y[edit]

Solution for y

This equation is precisely the equation of a parabola, that is, .

Drag and projectile motion[edit]

Figure 3: An example trajectory, the green line denotes the closest fit a parabola can give to the trajectory.

When drag is taken into consideration there are two major possible corrections to the force acting on an object, Stokes' and Newton. Stokes' gives the correction whereas Newton gives .[1] Therefore Newton's second law reduces to:

Letting:

Working in two dimensions (x and y), yields two equations:

   (1)

   (2)

Where and are the velocities in the x and y directions, respectively. These equations, being nonlinear coupled ordinary differential equations are impossible to solve in an exact (analytical) way, hence only numerical approximations are available. Letting, (refer to the [1] for the reasons for these choices) where the units are given here and and applying quadratic regression to the numerical solution gives the following MATLAB code and an estimate with a root mean square error of 0.068 metres.[1]

%dragcalc.m
clear all
options = odeset('RelTol',1e-10,'AbsTol',[1e-10 1e-10 1e-10 1e-10]);

%Using the ODE45 function http://www.mathworks.com.au/help/matlab/ref/ode45.html. 

[t,x]=ode45(@drag,[0 1.92],[0 10 0 10],options);

%Quadratic regression
N=length(t);
P=[ones(N,1) x(:,1) x(:,1).^2];
a=P\x(:,3); Pr=P*a; 
err=(Pr-x(:,3)); rms=sqrt(err'*err/N); 

%Plotting
plot(x(:,1),x(:,3),'r',x(:,1),Pr,'g')

Where the drag function is given by the following m file:

%drag.m
function dx=drag(t,x)
  g=9.8; alf=0.016; bet=1.024e-5; %alf=alpha; bet=beta
  dx(1,1)=x(2,1); %dx/dt
  dx(2,1)=-alf*sqrt(x(2,1)^2+x(4,1)^2)*x(2,1)-bet*x(2,1); %d^2x/dt^2
  dx(3,1)=x(4,1); %dy/dt
  dx(4,1)=-g-alf*sqrt(x(2,1)^2+x(4,1)^2)*x(4,1)-bet*x(4,1); %d^2y/dt^2
end

Proof that a satellite will travel in an elliptical path around a central mass[edit]

Figure 4: The elliptical orbit of a planet in orbit around the sun. The sun is at the focus of this ellipse.

Conservation of energy[edit]

Let the Hamiltonian (which is the total energy of the system) be:

Where T and V are the kinetic and potential energies, respectively. Assuming we have the central mass at r=0 (the focus of the trajectory, although we do not know this quite yet as we are still to prove that the trajectory is elliptical in the first place), then:

 

 

 

 

(1)

Where G is the gravitational constant, which is approximately, .[2] is the mass of the larger object, m is mass of the smaller object (like the particle in orbit). As always the dot above r or theta, indicates the derivative with respect to time, that is, the rate at which said coordinate changes with respect time. According to the law of the conservation of energy, this quantity should be static (i.e. constant) throughout the projectile's path.

Substitutions galore to make our job easier[edit]

Now it is time for some substitutions to make our lives easier. Let: , . According to the law of the conservation of angular momentum . Let: , which is also a constant, provided m is constant, then .

Manipulating the Hamiltonian to get r dot[edit]

Multiplying (1) by , gives:

Extensive substitutions[edit]

To simplify this further, let: then:

 

 

 

 

(2)

Let us call this expression on the right-hand side . We shall rearrange f so that it is easier for us to use later.

Let the expression in brackets be called . Making some further substitutions: . Then:

Let: then:

Two final substitutions are required, let: and .

Inverting this, for future reference, gives:

 

 

 

 

(3)

Then:

This, in turn gives:


 

 

 

 

(4)

Expressing theta in terms of rho[edit]

 

 

 

 

(5)

Integrating with respect to rho gives:

 

 

 

 

(6)

Substituting, into (6), (3) and (4) gives:

 

 

 

 

(7)

These constants cancel out, giving:

 

 

 

 

(8)

Back-substituting z gives:

hence:

Solution for θ

Taking the cosine of both sides, multiplying by omega over beta and reversing sides gives:

re-expressing in terms of r:

Solution for r

Where:

Mathematics/Physics (Old)[edit]

Mathematics/Physics (Old)

Schrodinger's Equation for a Linear Potential in 1d
I have not found the following elsewhere on the internet which I have found fustrating before when I needed to know this or would have liked to have known this. So I thought out of 'do unto others, what you would like others to do to you' I better put it somewhere and where better than here.

If:

Where is the Hamiltonian operator then if working on the semi-infinite interval the eigenfunctions (which in fact have been normalized) and eigenvalues are:


Where are the roots of the airy function

Schrodinger's Equation for a Spherically Symmetric Linear Potential


Perturbation theory yields using our solution in 1d with as the perturbation to the unperturbed Hamiltonian .


Which is first order accurate (in the perturbation) where is the one dimensional wave function (from the above section) with the input argument of r instead of x and:

Where is the energies of the unperturbed problem (that is the 1d problem with a linear potential).

Using what's given above and the following one can construct the three dimensional wave function.

Perturbation Theory for Nonlinear Differential Equations
Suppose:

Where is a parameter that can be varied from 0 to 1 where at the solution is known and is the value of alpha for which a solution is desired.
Subject to:


If denotes the solution for when then we can expand as:

In turn we may substitute this into our orignal equation:

Here we use a derivative approximation to linearise the above equation

Where the zero index is to denote that it is evaluated with and the u and in the index indicates partial differention with respect to said variables. This leads to a series of differential equations in terms of the corrections (by equating the powers of alpha) and now we should point out that the boundary conditions on each of these differential equations are zero except for the zeroth order equation which is idenitcal to the equation for . It should be noted however that in F we require that should appear in some finite order polynomial for perturbation theory to work.
If the polynomial inside F in terms of alpha that appears is linear then the following should hold:

Where in terms without alpha derivatives is set to 0.



1D PDE similar to the Wave Equation


such that
such that

Where:


Where are the zeros of the airy function and

From the two initial conditions stated earlier we can deduce that the coefficients and are respectively:


2D Sturm-Liouville Problem of the form of the 2D Simple Harmonic Oscillator problem from quantum mechanics
MATLAB code and figures for solving/viewing the solutions of the problem:

With k set to one the following MATLAB code may be constructed and executed to give the following graphs and eigenvalues.

clear all
N=50; NN=90;

k=1;
%Truncated domain method is being used to Solve the above problem the computational domain [-inf,inf] is 

%being truncated to [-L,L]
L=6; a=-L; b=L;

%T_mn =T_n(x_m) where T_n(x) denotes the Chebyshev Poly. of the second kind with the input argument of the 

%extrema of T_N (x) on the interval [-1,1],  D1 is the first order chebyshev differentiation matrix 

%transformed to the domain [a,b], D2 is just the 2nd order variant of D1, E1 and E2 are D1 and D2 

%respectively with the first and last rows and columns removed. x is the chebyshev extrema points mapped
%onto the domain [a,b]. xsub is x with endpoints removed.

%[T D1 D2 E1 E2 x xsub]=cheb(N,a,b); %(function code given below)

%Linearly spaced points xc are used at the end to improve the resolution of the solution.
xc=linspace(-1,1,NN+1)'; TT=cos(acos(xc)*(0:N));
%Solving the above problem and ordering the eigenvalues from biggest (in this case least negative) to smallest (or most negative)
I=eye(N-1); [xx yy]=meshgrid(xsub,xsub); X=diag(xx(:)); Y=diag(yy(:));
H=kron(E2,I)+kron(I,E2)-k*(X.^2+Y.^2);
[uu lam]=eig(H);
lam=diag(lam); [lam IX]=sort(real(lam),'descend');
u=zeros(N+1,N+1,(N-1)^2); uc=zeros(NN+1,NN+1,(N-1)^2);
uuc=zeros((NN+1)^2,(N-1)^2);
uue=zeros((N+1)^2,(N-1)^2);
for i=1:(N-1)^2
  u(2:N,2:N,i)=reshape(uu(:,IX(i)),N-1,N-1); u(1:N+1,1:N+1,i)=[zeros(1,N+1); zeros(N-1,1) u(2:N,2:N,i) zeros(N-1,1); zeros(1,N+1)];
  
end
%2d fast chebyshev transform used to determine the chebyshev series expansion coefficients. (code given below)
a=fct2d(u);

  
clear uu
%Getting the solution onto the equally spaced points xc
uuc=(kron(TT,TT))*a;
uu=reshape(uuc,(NN+1),NN+1,(N-1)^2);

  
[xx yy]=meshgrid(x,x); xc=L*linspace(-1,1,NN+1)';
[xxc yyc]=meshgrid(xc,xc);

function [T,D,D2,E1,E2,x,xsub] = cheb(N,a,b)
  if N==0, D=0; x=1; return, end
    x=-cos(pi*(0:N)/N)'; xsub=x(2:N); T=cos(acos(x)*(0:N));
%This part was borrowed from Nick Trevethan's book Spectral Methods in MATLAB (2001)
    c=[2; ones(N-1,1); 2].*(-1).^(0:N)';
    X=repmat(x,1,N+1);
    dX=X-X';
    D=(c*(1./c)')./(dX+eye(N+1));
    D=D-diag(sum(D'));
%------end of section borrowed--------
    D=2/(b-a)*D;
    E1=D(2:N,2:N);
    D2=D^2;
    E2=D2(2:N,2:N);
    x=(b-a)/2*x+(a+b)/2; xsub=(b-a)/2*xsub+(a+b)/2;
end

function a=fct2d(u)
[N , M, L]=size(u);

for i=1:L
ak=fct(fct(u(:,:,i))'); %fct and ifct's files are shown below
a(:,i)=ak(:);
end

end

function u=ifct2d(a)
[N2 L]=size(a);
N=sqrt(N2);
%it is assumed that the a is in the form of a N^2xL object
a=reshape(a,N,N,L);
%preallocating for speed

u=zeros(N,N,L);
for i=1:L
u(:,:,i)=ifct(ifct(a(:,:,i))');

end

end


function a=fct(f,l)
%Input of l=1 gives the Fast chebyshev transform at the points
%x=-cos(pi/N*((0:N-1)'+1/2)); if doesn't equal this value it is at the
%points x=-cos(pi/N*0:N);
if ~exist('l','var')
    l=0;
end
f=f(end:-1:1,:);
A=size(f); N=A(1); 
if exist('A(3)','var') && A(3)~=1
    for i=1:A(3)
        if l==1
            a(:,:,i)=sqrt(2/N)*dct(f(:,:,i));
            a(1,:,i)=a(1,:,i)/sqrt(2);
        else
            F=ifft([f(1:N,:,i);f(N-1:-1:2,:,i)]);
            a(:,:,i)=([F(1,:); 2*F(2:(N-1),:); F(N,:)]); %This code and the code below that is the function 

%ifct is in part derived from a file on MATLAB file exchange by Greg von Winckel (although not in its entirety)
        end
    end
else
    if l==1
            a=sqrt(2/N)*dct(f(:,:,i));
            a(1,:)=a(1,:)/sqrt(2);
    else
            F=ifft([f(1:N,:);f(N-1:-1:2,:)]);
            a=([F(1,:); 2*F(2:(N-1),:); F(N,:)]);
    end
a=real(a);
end


function f=ifct(a,l)
%l=1 x=-cos(pi/N*((0:N-1)'+1/2)) otherwise x=-cos(pi/N*(0:N)')
k=size(a); N=k(1); 
if ~exist('l','var')
    l=0;
end

if l==1

    f=idct(sqrt(N/2)*[a(1,:)*sqrt(2); a(2:end,:)]);

else
    A=a;
   F=fft([A(1,:); [A(2:N,:);A(N-1:-1:2,:)]/2]); 
 
    B=(F(1:N,:));
    f=B;
end
f=real(f(N:-1:1,:));

end

(Not sure If I'm required to do this but just in case there's any legal issues with omitting this I'm gonna place it in this page)
Copyright (c) 2009, Greg von Winckel All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

   * Redistributions of source code must retain the above copyright 
     notice, this list of conditions and the following disclaimer.
   * Redistributions in binary form must reproduce the above copyright 
     notice, this list of conditions and the following disclaimer in 
     the documentation and/or other materials provided with the distribution
     

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

1st eigenfunction of the above MATLAB Code
For this eigenfunction the corresponding eigenvalue is:
1st Eigenfunction of the 2D Simple Harmonic Oscillator Plan view.jpeg1st Eigenfunction of the 2D Simple Harmonic Oscillator 2nd perspective view.jpeg 2nd eigenfunction of the above MATLAB Code

2nd eigenfunction of the 2d simple harmonic oscillator.jpeg 3rd eigenfunction of the above MATLAB Code
(Note this is the same as the 2nd eigenvalue, because as indicated by the only 1 subscript for the eigenvalues and two for the eigenfunctions this problem is degenerate.)
3rd Eigenfunction of the 2D Simple Harmonic Oscillator plan view.jpeg 3rd Eigenfunction of the 2D Simple Harmonic Oscillator 2nd perspective view.jpeg

4th eigenfunction of the above MATLAB Code

4th Eigenfunction of the 2D Simple Harmonic Oscillator plan view.jpeg 4th Eigenfunction of the 2D Simple Harmonic Oscillator 2nd perspective view.jpeg 5th eigenfunction of the above MATLAB Code

5th Eigenfunction of the 2D Simple Harmonic Oscillator plan view.jpeg 5th Eigenfunction of the 2D Simple Harmonic Oscillator 2nd perspective view.jpeg 6th eigenfunction of the above MATLAB Code

6th Eigenfunction of the 2D Simple Harmonic Oscillator plan view.jpeg 6th Eigenfunction of the 2D Simple Harmonic Oscillator 2nd perspective view.jpeg

FEM and Domain Decomposition Methods


For this problem the following MATLAB code may be constructed (This code uses the finite element method with linear shape functions):

function [x u]=FEMH(p,q,rd,a,b,alpha,beta,N)
  %d/dx(p(x) du/dx)+q(x)u=r(x), u(a)=alpha; u(b)=beta

  x=linspace(a,b,N+1)';
function L=q1f(x)
  xi=(x-x(1))/(x(end)-x(1));
    H1=(1-xi);
    L=q(x).*H1.^2;
end
function l=pp(x)
    l=p(x);
end
function L=q2f(x)
  xi=(x-x(1))/(x(end)-x(1));
    H2=xi;
    L=q(x).*H2.^2;
end
function L=q3f(x)
  xi=(x-x(1))/(x(end)-x(1));
    H1=1-xi; H2=xi;
    L=q(x).*H1.*H2;
end
function L=r1f(x)
  xi=(x-x(1))/(x(end)-x(1));
    H1=1-xi;
    [~, dl]=p(x);
    r=rd(x)-(beta-alpha)/(b-a)*(dl+x-a)+alpha;
    L=r.*H1;
  end
function L=r2f(x)
  xi=(x-x(1))/(x(end)-x(1));
  H2=xi;
  L=r(x).*H2;
end
p1=zeros(N,1); q1=zeros(N,1); q2=zeros(N,1); q3=zeros(N,1); r1=zeros(N,1); r2=zeros(N,1); F=zeros(N-1,1);
  for i=1:N
  p1(i)=N^2/(b-a)^2*quad(@pp,x(i),x(i+1));
  q1(i)=quad(@q1f,x(i),x(i+1));
  q2(i)=quad(@q2f,x(i),x(i+1));
  q3(i)=quad(@q3f,x(i),x(i+1));
  r1(i)=quad(@r1f,x(i),x(i+1));
  r2(i)=quad(@r2f,x(i),x(i+1));
  if i~=N
  F(i)=-(r1(i)+r2(i));
  end
  end
  p11=p1(2:N);
  p1=p1(1:N-1); q2=q2(1:N-1);
  q11=q1(2:N);
  q31=q3(2:N);
  
    K=diag(p1+p11-q11-q2);
   
    K=K+diag(q31(1:N-2)-p11(1:N-2),1)+diag(q31(1:N-2)-p11(1:N-2),-1);
 
  
  uu=K\F;
 v=[0; uu; 0];
 [~, dl]=p(x);
 
 u=v+(beta-alpha)/(b-a)*(dl+x-a)+alpha; 
end
function [x u]=SubDomCheb(N,M,G,bd,A)
  %This function solves the problem:
  %d/dx (p(x) du/dx)+q(x)u=r(x)
  %Where the function G is defined such that:
  %[p q r]=G(x,A);
  %alpha1*u(a)+alpha2*u'(a)=alpha3, beta1*u(b)+beta2*u'(b)=beta3;
  %using domain decomposition method with chebyshev cardinal functions.

[T, D1]=cheb(N,-1,1);
x=-cos(pi*(0:N)'/N);
a=bd(1); alpha1=bd(2); alpha2=bd(3); alpha3=bd(4); b=bd(5); beta1=bd(6); beta2=bd(7); beta3=bd(8);


D1=2*M/(b-a)*D1; 
xx=zeros(N+1,M);
for i=1:M
  xx(:,i)=(b-a)/(2*M)*(x+1)+(b-a)/M*(i-1)+a;
  xt((N+1)*(i-1)+1:i*(N+1),1)=xx(:,i);
end
[p, q, r]=G(xt,A);
F=r;
clear p q

for i=1:M
  if i==1
    clear H
    [p q]=G(xx(:,i),A);
    H=D1*diag(p)*D1+diag(q);
    H(1,1:N+1)=alpha1*[1 zeros(1,N)]+alpha2*D1(1,:); H(N+1,1:N+1)=D1(N+1,:);
    H(N+1,(N+2):2*(N+1))=-D1(1,:);
    K(1:N+1,1:2*(N+1))=H;
    F(N+1,1)=0; F(1)=alpha3;
  elseif i==M
    clear H
    [p q]=G(xx(:,i),A);
    H(:,(N+2):2*(N+1))=D1*diag(p)*D1+diag(q);
    H(1,1:2*(N+1))=[zeros(1,N) 1 -1 zeros(1,N)];
    H(N+1,(N+2):2*(N+1))=beta1*[zeros(1,N) 1]+beta2*D1(N+1,:);
    K((i-1)*(N+1)+1:i*(N+1),(i-2)*(N+1)+1:i*(N+1))=H;
    F((i-1)*(N+1)+1)=0; F(i*(N+1))=beta3;
  else
    clear H
    [p q]=G(xx(:,i),A);
    H(1:(N+1),(N+2):2*(N+1))=D1*diag(p)*D1+diag(q);
    H(1,1:2*(N+1))=[zeros(1,N) 1 -1 zeros(1,N)];
    H((N+1),(N+2):3*(N+1))=[D1(N+1,:) -D1(1,:)];
    K(((i-1)*(N+1)+1):(i*(N+1)),((i-2)*(N+1)+1):((i+1)*(N+1)))=H;
    F((i-1)*(N+1)+1)=0; F(i*(N+1))=0;
  end
end

uu=K\F;
xx=linspace(-1,1,N+1)'; TT=cos(acos(xx)*(0:N));
for i=1:M
  aa(:,i)=fct(uu((i-1)*(N+1)+1:i*(N+1)));
  us(:,i)=TT*aa(:,i);
  u((i-1)*(N+1)+1:i*(N+1),1)=us(:,i);
  xs(:,i)=(b-a)/(2*M)*(xx+1)+(b-a)/M*(i-1)+a;
  x((i-1)*(N+1)+1:i*(N+1),1)=xs(:,i);
end
end

PDEs


And Sturm-Liouville problem boundary conditions with respect to the x coordinate.
Where is a Sturm-Liouville operator. That is it satisfies:

Subject to standard SL boundary conditions. And of course the orthogonality condition, however for convenience an additional normality condition will be imposed on the norms of each eigenfunction, .






SLE

the solution is:

Reference list[edit]

  1. ^ a b c Taylor, JR (2005). "Chapter 2: Projectiles and Charged Particles". Classical Mechanics (PDF)|format= requires |url= (help). Sausalito, USA: University Science Books. ISBN 978-1891389221.
  2. ^ Mohr, PJ; Taylor, BN; Newell, DB (2011). Baker, J; Douma, M; Kotochigova, S (ed.). "The 2010 CODATA Recommended Values of the Fundamental Physical Constants" (Web Version 6.0)". Gaithersburg, MD: National Institute of Standards and Technology.CS1 maint: multiple names: authors list (link)