readme

来自「various matlab programs to slove various」· 代码 · 共 369 行

CHAPTER 6:

1) Example of usage of Program newthorn:

suppose to find the roots of:

f(x) = x^3 + 5x^2 -x + 1

Matlab instructions:

>> A=[1, 5, -1, 1];
>> [xn,iter,root]=newthorn(A,3,1e-5,0,100,0);  
>> iter

iter =

    11     9     1

>> root

root =

  Column 1 

     -5.227871411936592e+00                         

  Column 2 

      1.139357361976039e-01 - 4.222575255986668e-01i

  Column 3 

      1.139357375576224e-01 + 4.222575252177587e-01i

>> r=roots(A)

r =

     -5.227871411936587e+00                         
      1.139357059682952e-01 + 4.222571431798804e-01i
      1.139357059682952e-01 - 4.222571431798804e-01i

>> [xn,iter,rootN,itrefin]=newthorn(A,3,1e-5,0,100,1);
>> itrefin

itrefin =

     1     2     2

>> root

rootN =

  Column 1 

     -5.227871411936590e+00                         

  Column 2 

      1.139357059682953e-01 - 4.222571431798806e-01i

  Column 3 

      1.139357059682953e-01 + 4.222571431798806e-01i

(Compare the accuracy between r, root and rootN)

2) Example of usage of Program mulldefl:

(same polynomial as before)

>> [xn,iter,root]=mulldefl(A,3,1e-5,0,1,2,100,0);        
>> iter

iter =

     5     2     1

>> root

root =

  Column 1 

      1.139357060245700e-01 - 4.222571432421807e-01i

  Column 2 

      1.139357059225077e-01 + 4.222571432502485e-01i

  Column 3 

     -5.227871411947078e+00 - 8.067824186497319e-12i

>> [xn,iter,rootN,itrefin]=mulldefl(A,3,1e-5,0,1,2,100,1);
>> itrefin

itrefin =

     1     1     1

>> rootN

rootN =

  Column 1 

      1.139357059682953e-01 - 4.222571431798806e-01i

  Column 2 

      1.139357059682952e-01 + 4.222571431798806e-01i

  Column 3 

     -5.227871411936590e+00 + 2.465190328815662e-32i

(Check again the accuracy of the Newton refinement)


CHAPTER 7:

1) Example of usage of Program newtonsys:

suppose to compute the solution of the nonlinear
system 

x^2 + y^2 - 1 = 0
x - y         = 0

The Matlab instructions are the following

>> F=['x(1)^2+x(2)^2-1';...               
      'x(1)-x(2)      ']

F =

x(1)^2+x(2)^2-1
x(1)-x(2)      

>> J=['2*x(1)';...
      '2*x(2)';...
      '1     ';...
      '-1    ']

J =

2*x(1)
2*x(2)
1     
-1    

(be careful that each row of F and J have the same sizes!!)

>> [x, nit] = newtonsys(F, J, [0, -1]', 1e-10, 100, 0);      
>> x

x =

    -7.071067811742358e-01
    -7.071067811742358e-01

>> nit

nit =

    28

2) Example of usage of program broyden:

again, be careful to define matrix Q (nxn) and
the vector f (nx1). Using the same data as in the previous example
we get:

>> f=['x(1)^2+x(2)^2-1';...
      'x(1)-x(2)      '];         
>> Q=rand(2);
>> [x,it]=broyden(x,Q,nmax,toll,f);
>> x

x =

    0.7071
    0.7071

>> it

it =

     3

CHAPTER 9:

1) There is a small bug in the Program newtcot:

the line: 

for j=2:n+1,    x=x+h; y(j)=eval(f); end;    int=0; 

must be changed into:

for j=2:n+1,    x=x+h; y(j)=eval(fun); end;    int=0;

otherwise Matlab will echo the following error message:

>> int = newtcot(a,pi/2,1,fun)
??? Undefined function or variable 'f'.

Error in ==> /tmp_mnt/ipmma19/home3/calnum/PerSpringerNY/Programs/newtcot.m
On line 12  ==> for j=2:n+1,    x=x+h; y(j)=eval(f); end;    int=0; 

2) Here is an example of usage of Program trapmodc.
Suppose to compute the integral of 

f(x)=cos(x)

between -pi/2 and pi/2 (integral=2).
You must first define the function dfun.m

function y=dfun(x)
y=-sin(x);

and then the program can be executed as follows:

>> fun='cos(x)'; dfun='dfun'; a=-pi/2; b=pi/2; m=10;
>> int = trapmodc(a,b,20,fun,dfun)                  

int =

     1.999998307875836e+00
    
3) Here is an example of usage of Programs midptr2d and
traptr2d (two-dimensional integration).
Suppose to compute the integral of

f(x,y)=y/(1+xy)

over the square (0,1)^2 (the exact integral is log(4)-1 = 0.3863)

The Matlab instructions can be as follows:

>> x(1)=0; x(2)=1/2; x(3)=1;
>> x(4)=0; x(5)=1/2; x(6)=1;
>> x(7)=0; x(8)=1/2; x(9)=1;
>> y(1)=0; y(2)=0; y(3)=0; 
>> y(4)=0.5; y(5)=0.5; y(6)=0.5;
>> y(7)=1; y(8)=1; y(9)=1;
(Coordinates of the integration points)

>> lel(1,1)=1; lel(1,2)=2; lel(1,3)=5;
>> lel(2,1)=1; lel(2,2)=5; lel(2,3)=4;
>> lel(3,1)=2; lel(3,2)=3; lel(3,3)=6;
>> lel(4,1)=2; lel(4,2)=6; lel(4,3)=5;
>> lel(5,1)=4; lel(5,2)=5; lel(5,3)=8;
>> lel(6,1)=4; lel(6,2)=8; lel(6,3)=7;
>> lel(7,1)=5; lel(7,2)=6; lel(7,3)=9;
>> lel(8,1)=5; lel(8,2)=9; lel(8,3)=8;
 
(local numbering of the mesh nodes)

>> fun='y./(1+x.*y)';
>> inte=0;
>> for ie=1:8, i=lel(ie,1); j=lel(ie,2); k=lel(ie,3);
   xv(1)=x(i); xv(2)=x(j); xv(3)=x(k);
   yv(1)=y(i); yv(2)=y(j); yv(3)=y(k);
   inte=inte+midptr2d(xv,yv,fun);
   end
>> inte

inte =

    0.3918

>> inte=0;
>> for ie=1:8, i=lel(ie,1); j=lel(ie,2); k=lel(ie,3);
   xv(1)=x(i); xv(2)=x(j); xv(3)=x(k);
   yv(1)=y(i); yv(2)=y(j); yv(3)=y(k);
   inte=inte+traptr2d(xv,yv,fun);
   end
>> inte

inte =

    0.3708

CHAPTER 10:

There is a bug in Remark 10.3. 
Instead of 

x=phi(xi)=(a+b)/2*xi + (b-a)/2

one should read:

x=phi(xi)=(b-a)/2*xi + (a+b)/2. 

Therefore, the formulae yielding the 
integration nodes and weights
after the affine change of coordinate are:

x_j = (b-a)/2*xi_j + (a+b)/2
alpha_j = (b-a)/2 * beta_j

An example of this is the compitation of 
the integral on [-pi/2, pi/2] of f(x)=cos(x) 
(exact integral=2). Using the Program zplege we get

>> n=5;            
>> [t,w]=zplege(n);
>> x=pi/2*t;
>> w=pi/2*w;
>> int

int =

     2.000000110284471e+00

>> 

CHAPTER 11:
	
1) Example of usage of Program multistep.
We consider the Cauchy problem 

y'(t)=-y(t)      0 < t < 10
y(0)=1

and employ the Crank-Nicolson method to solve it.
The Matlab instructions are as follows:

>> a=[1]; b=[1/2; 1/2]; tf=10; t0=0; u0=1;               
>> h=0.1; dfun='-1+0.*y'; fun='-y'; tol=1e-10; itmax=100;
>> [t,u] = multistep (a,b,tf,t0,u0,h,fun,dfun,tol,itmax);

Executing the instruction 

>> plot(t,u) 

will show that the scheme has worked correctly. 


2) Example of usage of Program predcor.
Solve the same problem as above using the Heun method
(i.e.: Forward Euler as predictor, Crank-Nicolson as corrector, 
with m=1, i.e., only one correction per step).
Be quite careful in the definition of the coefficient 
vectors at and bt (predictor) and a and b (for the corrector).
Precisely, the size of b should always be greater than the
size of the corresponding vector bt.

>> a = [1];
>> b = [1/2 1/2]';
>> at = a;
>> bt = [1];
>> h = 0.1;
>> f = '-y';
>> t0=0;
>> u0=1;
>> tf = 10;
>> pece='n';
>> m = 1;
>> [u,t]=predcor(a,b,at,bt,h,f,t0,u0,tf,pece,m);

Again, executing the instruction

>> plot(t,u)

will evidence the right behaviour of the method.