Exerciţii programate în Scilab

1. Să se adune primii 20 de termeni din seria

1

2

!1k k

k

//ex 1 din lectia 1s=1; k=1; t=1;for k=2:20

t=t*k^2/(k-1)^3;s=s+t;

end

disp(s)

//gata

2. Să se adune termenii mai mari decât 710 din seria

1

2

!1k k

k .

//ex 2 din lectia 1s=0; k=1; t=1; eps=10^(-7);while (t>eps)

s=s+t;k=k+1;t=t*k^2/(k-1)^3;

ends,k

//gata

3. Un sistem de 100 ecuaţii cu 100 are matricea egală cu 4 pe diagonalaprinciupală iar diagonalele de deasupra şi dedesuptul diagonalei principale suntegale cu 1. Celelalte elemente sunt egale cu zero. Termenii liberi sunt egali cu 1. Săse rezolve sistemul.

//un sistemclear;n=100;a=zeros(100,100);for i=1:n

a(i,i)=4;if(i<n) then

a(i,i+1)=1;endif (i>1) then

a(i,i-1)=1;

endb(i)=1;

end

x=a\b;disp(x);

//gata

4. Să se reprezinte pe aceeaşi figură dar pe grafice diferite funcţiile f(x)=sin(kx)pentru k=1,2,3,4 şi ],0[ x .

//grafice de functii planeclear;x=linspace(0,%pi,100);y1=sin(x);subplot(221);plot(x,y1,'r:');xgrid;xtitle('SINUSURI')y2=sin(2*x);subplot(222);plot(x,y2,'b-');xgrid;y3=sin(3*x);subplot(223)plot(x,y3,'k-');xgrid

5. Să se arate că ecuaţia 01sin3 xax are pentru orice 0a o singură soluţie înintervalul [0, 1]. Cum se poate determina această soluţie în Scilab? Să se reprezintegrafic aceste soluţii în funcţie de a pentru ]1,0[a

//rezolvam o ecuatiex0=0.5;function y=g(a)

function z=f(x)z=x^3+a*sin(x)-1

endfunctiony=fsolve(x0,f)

endfunction

a=linspace(0,1,20);for i=1:20

x(i)=g(a(i));end

plot(a,x)xgrid;xtitle('radacina in functie de a')

6. Fie x un vector cu n=200 componente. Să se calculeze:

a. 2

3

ixix

b. 0sin ix

ix

c. 11

21ix i

i

x

x

// niste sumeclear;n=200000;x=10*(rand(n,1)-0.5);tic();//metoda 1s=0;for i=1:n

if(x(i)>2)s=s+x(i)^3;

endendtoc()disp(s)

//metoda 2tic()s1=sum(x(x>2).^3)toc()

7. Să se facă graficul suprafeţei:

22

22

vuz

uvy

vux

,11

11

v

u

function [x, y, z]=sup(u, v)x=u.^2+v.^3;y=u+v;z=u.^2-v.^2;

endfunction

u=linspace(-1,1,20);v=linspace(-1,1,20);[xx,yy,zz]=eval3dp(sup, u, v);cc=(zz-min(zz))*128/max(zz-min(zz));//se determina culoarea in functie de zplot3d(xx,yy,list(zz,cc));figura=gcf();figura.color_map = rainbowcolormap(128);

8. Să se determine epsilon maşină pentru calculatorul pe care lucrăm.Rezolvare.

-->format("e",22)

-->nearfloat("succ",1) - 1 ans =

2.220446049250313D-16Prin comenzile următoare vedem care este distanţa între două numere realeconsecutive reprezentabile în calculador:

-->nearfloat("succ",10^10)-10^10 ans =

1.907348632812500D-06

-->nearfloat("succ",10^15)-10^5ans =

9.999999999000001D+14

10. Se dă tabelul:

x 1 3 4 5y -2 4 1 6

Să se scrie polinomul de interpolare sub forma Lagrange. Să se programezecalculul polinomului de interpolare sub forma Lagrange.Rezolvare.

//interpolare Lagrangex=[1;3;4;5];y=[-2;4;1;6];

function v=lag(t, x, y)v=0;n=length(x);for i=1:n

p=y(i);for j=1:n

if(j~=i)p=p*(t-x(j))/(x(i)-x(j));

endendv=v+p;

endendfunction

//facem graficul puctelorplot(x,y,'o')//graficul polinomuluivx=linspace(x(1),x(length(x)),100);for i=1:length(vx)

vy(i)=lag(vx(i),x,y);endplot(vx,vy,'r-')xgrid

Rezultatul rulării este graficul următor:

11. Să se studieze cât de bine este aproximată funcţia cos(x) pe intervalul ],0[ depolinomul Lagrange cu grad din ce in ce mai mare.,Rezolvare.Definim

n=10;x=linspace(0,%pi,n);y=cos(x);

In rest procedăm ca în programul anterior.Pentru n=10 graficul valorilor exacte şi al celor înterpolate practic se suprapun:

Pentru n=70 vedem că apar mari diferenţe la capete datorită erorilor de3 rotunjire:

12. Să se programeze calculul valorii polinomului de interpolare prin metodaNewton şi să se compare cu valorile obţinute prin metoda lui Lagrange.

//interpolare Newtonx=[1;3;4;5];y=[-2;4;1;6];

function d=dd(x, y)n=length(x);d=y;for k=1:n-1

for i=0:n-k-1d(n-i)=(d(n-i)-d(n-i-1))/(x(n-i)-x(n-i-k));

endend

endfunction

13. Să se facă graficul funcţiei spline naturală ce interpolează tabelul

x -2 -1 1 3 4 5y 2 3 5 4 0 -1

//splinex=[-2,-1,1,3,4,5];y=[2, 3, 5, 4, 0, -1];plot(x,y,'o');xgridd=splin(x,y,"natural");xx=linspace(x(1),x(length(x)),100);yy=interp(xx,x,y,d);plot(xx,yy,'r')

14. Metoda celor mai mici pătrate

15. Derivare numerică

16. Să se arate că 6

2sin2

xxxf este contracţie pe [-1, 1]. Să se programeze

determinarea punctului fix.//contractie

function y=f(x);y=(x^2+sin(x)-2)/6;

endfunction

xi=2;n=20;x=zeros(n,1);x(1)=xi;for i=2:n

x(i)=f(x(i-1))end

disp(x);

17. Să se programeze rezolvarea sistemului

05

0444

33

yxx

yxyx prin iteraţii Newton.

//rezolvarea sistemelor neliniaren=2; //nr. de necunoscutenmax=200;

function y=f(x) //sistemul de ecuatiiy=zeros(n,1)y(1)=x(1)^3+x(1)*x(2)+x(2)^3-4y(2)=x(1)^4+x(1)+x(2)^4-5endfunction

eps=0.0000001;//criteriul de oprire er=||xi-xu||<epsh=0.0001 ; //pasul pentru calculul derivatelorxi=[1;1] ; //x initialer=eps+1 ; // ne asiguram ca er>eps la inceputa=zeros(n,n); //a va contine derivatele partialecontor=1;while ((er>eps) | (contor>nmax))

contor=contor+1;disp(xi);val=f(xi);for i=1:nxi(i)=xi(i)+h;a(:,i)=(f(xi)-val)/h;xi(i)=xi(i)-h;endxu=xi-inv(a)*f(xi);er=norm(xi-xu);xi=xuenddisp(xi)

18. Să se rezolve sistemul

05

0444

33

yxx

yxyx prin procedura fsolve din Scilab.

//rezolvarea sistemelor neliniare

function y=f(x) //sistemul de ecuatiiy=zeros(2,1)y(1)=x(1)^3+x(1)*x(2)+x(2)^3-4y(2)=x(1)^4+x(1)+x(2)^4-5endfunction

xi=[1;1] ;sol=fsolve(xi,f)

19. Să se programeze rezolvarea sistemului

052

34

24

zyx

zyx

zyx

prin metoda Iacobi.

//metoda Iacobia=[4, -1, 1; 1,4,-1; -1,2,5]b=[2;3;0]D=diag(diag(a));A=-inv(D)*(a-D)B=inv(D)*b;

function y=f(x)y=A*x+B

endfunction

x1=[1;1;1]for i=1:20

x1=f(x1)end

20. Să se programeze calculul integralei

5

3

2

21

1sindx

x

x prin metoda trapezelor,

metoda Simpson şi procedura integrate din Scilab.a.

//integrare metoda trapezelor

function y=f(x)y=sin(x^2+1)/(1+sqrt(2+x))

endfunction

a=3b=5n=10;x=linspace(a,b,n+1);s=0;for i=2:n

s=s+f(x(i))end

rez=(f(a)+f(b)+2*s)*(b-a)/(2*n)disp(rez)

b.//metoda lui Simpson

function y=f(x)y=sin(x^2+1)/(1+sqrt(2+x))

endfunction

a=3b=5n=10;x=linspace(a,b,n+1);

s1=0;for i=2:n

s1=s1+f(x(i))end

s2=0;for i=1:n

s2=s2+f((x(i)+x(i+1))/2)

end

rez=(f(a)+f(b)+2*s1+4*s2)*(b-a)/(6*n)disp(rez)

c.//utilizam procedura de integrare din scilab

function y=f(x)y=sin(x^2+1)/(1+sqrt(2+x))

endfunction

rez=integrate('f(x)','x',3,5)disp(rez)

21. Să se calculeze

D

dxdyyx

xyx22

2

1cos , 11,20,, 2 xyxxyxD prin

procedura integrate din Scilab.

//integrale duble

function z=f(x, y)z=(x*cos(y)+x^2)/(1+x^2+y^2)

endfunction

function y=jos(x)y=x-1

endfunction

function y=sus(x)y=x^2+1

endfunction

function val=asta(x)val=integrate('f(x,y)','y',jos(x),sus(x))

endfunction

a=0; b=2;rez=integrate('asta(x)','x',a,b)

disp(rez)

22. Să se calculeze integrala triplă

V

dxdydzzyx

zxy

1, prin procedura integrate

din Scilab, unde 11,11,1,,, 2222 xxyxzyxzyxV

//integrala tripla

function u=f(x, y, z)u=(x*y+z)/(1+x+y+z)

endfunction

function z=zjos(x, y)z=x^2+y^2

endfunction

function z=zsus(x, y)z=1

endfunction

function y=jos(x)y=-sqrt(1-x^2)

endfunction

function y=sus(x)y=sqrt(1-x^2)

endfunction

a=-1b=1

function v=intz(x, y)v=integrate('f(x,y,z)','z',zjos(x,y),zsus(x,y))

endfunction

function v=inty(x)v=integrate('intz(x,y)','y',jos(x),sus(x))

endfunction

tic()rez=integrate('inty(x)','x',a,b)

disp(rez)disp(toc())

23. Să se programeze problema Cauchy xeyy ' , 00 y prin metodele Euler,Euler modificată, Euler îmbunătăţită, Euler predictor-corector, Runge-Kutta deordinal 4.

//metode numerice pentru ecuatii diferentialeclearfunction z=f(x, y)

z=y+exp(x)endfunction

//conditiile initialexi=0; yi=0;//intervalulL=1;n=10;h=L/n;x=linspace(xi,xi+L,n+1);

//metoda Euleryeuler=zeros(1,n+1);yeuler(1)=yi;for i=1:n

yeuler(i+1)=yeuler(i)+h*f(x(i),yeuler(i));endplot(x,yeuler,'r-');

//metoda Euler modificatayemodif=zeros(1,n+1);yemodif(1)=yi;for i=1:n

yemodif(i+1)=yemodif(i)+h*f(x(i)+h/2,yemodif(i)+h/2*f(x(i),yemodif(i)));endplot(x,yemodif,'g-.')

//Euler imbunatatitayeimb=zeros(1,n+1);yeimb(1)=yi;for i=1:n

yeimb(i+1)=yeimb(i)+h/2*(f(x(i),yeimb(i))+f(x(i)+h,yeimb(i)+h*f(x(i),yeimb(i))));endplot(x,yeimb,'y:')

//Euler predictor-corectoreps=0.0000000001;yepc=zeros(1,n+1);yepc(1)=yi;

for i=1:nyp=yepc(i)+h*f(x(i)+h/2,yepc(i)+h/2*f(x(i),yepc(i)));yc=yepc(i)+h/2*(f(x(i),yepc(i))+f(x(i)+h,yp));while (abs(yc-yp)/(1+abs(yp))>eps)

yp=yc;yc=yepc(i)+h/2*(f(x(i),yepc(i))+f(x(i)+h,yp));

endyepc(i+1)=yc;

endplot(x,yepc,'k:');

//metoda Runge-Kuttayrk4=zeros(1,n+1);yrk4(1)=yi;for i=1:n

k1=h*f(x(i),yrk4(i));k2=h*f(x(i)+h/2,yrk4(i)+k1/2);k3=h*f(x(i)+h/2,yrk4(i)+k2/2);k4=h*f(x(i)+h,yrk4(i)+k3);yrk4(i+1)=yrk4(i)+(k1+2*k2+2*k3+k4)/6;

endplot(x,yrk4,'c-')

//y exactyexact=zeros(1,n+1);for i=1:n+1

yexact(i)=x(i)*exp(x(i));endplot(x,yexact,'b:')

xgrid

23. Să se rezolve problema Cauchy

2122

2111

02.03.0'

1.01.0'

yyyy

yyyy, 200,200 21 yy pe un

interval de lungime L=200. Programele se vor face pentru metodele: Euler, Eulermodificată, Euler îmbunătăţită, Runge-Kutta de ordinal 4. Să se facă graficul curbei 2000, 21 xxyxy .//metoda numerice pentru sisteme de ecuatii diferentialeclearnec=2;function z=f(x, y)

z=zeros(nec,1);z(1)=0.1*y(1)-0.01*y(1)*y(2)z(2)=-0.3*y(2)+0.02*y(1)*y(2)

endfunction

//conditiile initialexi=0; yi=[20;20];//intervalulL=100;n=200;h=L/n;

x=linspace(xi,xi+L,n+1);

//metoda Euleryeuler=zeros(nec,n+1);yeuler(:,1)=yi;for i=1:n

yeuler(:,i+1)=yeuler(:,i)+h*f(x(i),yeuler(:,i));end

//metoda Euler modificatayemodif=zeros(nec,n+1);yemodif(:,1)=yi;for i=1:n

yemodif(:,i+1)=yemodif(:,i)+h*f(x(i)+h/2,yemodif(:,i)+h/2*f(x(i),yemodif(:,i)));end

//Euler imbunatatitayeimb=zeros(nec,n+1);yeimb(:,1)=yi;for i=1:n

yeimb(:,i+1)=yeimb(:,i)+h/2*(f(x(i),yeimb(:,i))+f(x(i)+h,yeimb(:,i)+h*f(x(i),yeimb(:,i))));end

//metoda Runge-Kuttayrk4=zeros(nec,n+1);yrk4(:,1)=yi;for i=1:n

k1=h*f(x(i),yrk4(:,i));k2=h*f(x(i)+h/2,yrk4(:,i)+k1/2);k3=h*f(x(i)+h/2,yrk4(:,i)+k2/2);k4=h*f(x(i)+h,yrk4(:,i)+k3);yrk4(:,i+1)=yrk4(:,i)+(k1+2*k2+2*k3+k4)/6;

end

//plot(x,yrk4(1,:),'r');//plot(x,yrk4(2,:),'b');//plot(x,yrk4(3,:),'k')plot(yrk4(1,:),yrk4(2,:))

xgrid

Se obţine următorul grafic:

24. Să se rezolve problema Cauchy

2122

2111

02.03.0'

1.01.0'

yyyy

yyyy, 200,200 21 yy pe un

interval de lungime L=200. Programele se vor face pentru metodele Adams-Bashfort şi Adams-Moulton. Paşii iniţiali se vor calcula prin metoda RK4. Să sefacă graficul curbei 2000, 21 xxyxy .

//metoda multipas pentru ecuatii diferentialeclearnec=2;function z=f(x, y)

z=zeros(nec,1);z(1)=0.1*y(1)-0.01*y(1)*y(2)z(2)=-0.3*y(2)+0.02*y(1)*y(2)

endfunction

//conditiile initialexi=0; yi=[20;20];//intervalulL=1000;n=20000;h=L/n;x=linspace(xi,xi+L,n+1);

//primii 3 pasiymp=zeros(nec,n+1);ymp(:,1)=yi;for i=1:3

k1=h*f(x(i),ymp(:,i));k2=h*f(x(i)+h/2,ymp(:,i)+k1/2);k3=h*f(x(i)+h/2,ymp(:,i)+k2/2);k4=h*f(x(i)+h,ymp(:,i)+k3);ymp(:,i+1)=ymp(:,i)+(k1+2*k2+2*k3+k4)/6;

end

//Adams Bashfort de ordinul 2for i=4:n

ymp(:,i+1)=ymp(:,i)+h/2*(3*f(x(i),ymp(:,i))-f(x(i-1),ymp(:,i-1)))

endyab2=ymp;//plot(yab2(1,:),yab2(2,:));//xgrid

//Adams Bashfort de ordinul 3for i=4:n

ymp(:,i+1)=ymp(:,i)+h/12*(23*f(x(i),ymp(:,i))-16*f(x(i-1),ymp(:,i-1))+5*f(x(i-2),ymp(:,i-2)));endyab3=ymp;//plot(yab3(1,:),yab3(2,:));xgrid;

//Adams Bashfort de ordinul 4for i=4:n

ymp(:,i+1)=ymp(:,i)+h/24*(55*f(x(i),ymp(:,i))-59*f(x(i-1),ymp(:,i-1))+37*f(x(i-2),ymp(:,i-2))-9*f(x(i-3),ymp(:,i-3)));endyab4=ymp;//plot(yab4(1,:),yab4(2,:));//xgrid;

//Adams Moulton de ordinul 3eps=0.0000000001;for i=4:n

yp=ymp(:,i)+h/24*(55*f(x(i),ymp(:,i))-59*f(x(i-1),ymp(:,i-1))+37*f(x(i-2),ymp(:,i-2))-9*f(x(i-3),ymp(:,i-3)));yc=ymp(:,i)+h/12*(5*f(x(i+1),yp)+8*f(x(i),ymp(:,i))-f(x(i-1),ymp(:,i-1)));while (abs(yc-yp)/(1+abs(yp))>eps)

yp=yc;yc=ymp(:,i)+h/12*(5*f(x(i+1),yp)+8*f(x(i),ymp(:,i))-f(x(i-1),ymp(:,i-1)));

endymp(:,i+1)=yc;

endyam3=ymp;//plot(yam3(1,:),yam3(2,:));xgrid

//Adams Moulton de ordinul 4eps=0.0000000001;for i=4:n

yp=ymp(:,i)+h/24*(55*f(x(i),ymp(:,i))-59*f(x(i-1),ymp(:,i-1))+37*f(x(i-2),ymp(:,i-2))-9*f(x(i-3),ymp(:,i-3)));yc=ymp(:,i)+h/24*(9*f(x(i+1),yp)+19*f(x(i),ymp(:,i))-5*f(x(i-1),ymp(:,i-1))+f(x(i-2),ymp(:,i-2)));while (abs(yc-yp)/(1+abs(yp))>eps)

yp=yc;yc=ymp(:,i)+h/24*(9*f(x(i+1),yp)+19*f(x(i),ymp(:,i))-5*f(x(i-1),ymp(:,i-1))+f(x(i-2),ymp(:,i-2)));

endymp(:,i+1)=yc;

endyam4=ymp;plot(yam4(1,:),yam4(2,:));xgrid