# Wprowadzenie do  Maximy

## Wstêp

### Maxima Interfaces

### How to Read this Tutorial

input;

/* clear working space */
 reset()$ kill(all)$

## Podstawy

### Jak zacz¹æ

2+3;

4*5$

2+3; 4*5;

2+3*(4+5)
+6*7-8/4;

/* clear working space */
 reset()$ kill(all)$

### wxMaxima

### Maxima jako Calculator

### Podstawowe dzia³ania 

5+3;

22-7;

121/11;

25*4;

2^8;

(2+3)^4-(5+6+7)/(2+4);

2^100;

/* clear working space */
 reset()$ kill(all)$

### Sta³e

%pi; %e; %inf;

### Funkcje ró¿ne

cos(%pi);

exp(2);

abs(-5);

10^4;

1000^(1/3);

8^(2/3);

log(%e);

ln(%e);

/* clear working space */
 reset()$ kill(all)$

### Zmienne i funkcje u¿zytkownika

x;

x: 2+3;

x;

y: 2*x^2;

sqr(x):= x^2;

sqr(4);

Przyk³ad definicji funkcji zamieniaj¹ce stopnie na radiany 

deg2rad(x):= x*%pi/180;

sin(deg2rad(45));

norm(x,y):= sqrt(x^2+y^2);

norm(3,4);

values;

functions;

kill(x);

values;

functions;

kill(all);

values;

functions;

x=3;

x;

/* clear working space */
 reset()$ kill(all)$

### Wykorzystanie wyjœcia jako wejœcia do dalszych obliczeñ

2+3;

% * 20;

% - 10;

%o1 - 10;

/* clear working space */
 reset()$ kill(all)$

### Help - wykorzystanie

? sqrt;

?sqrt;

apropos("sqrt");

/* clear working space */
 reset()$ kill(all)$

### Zmienne systemowe

### Liczby - manipulacje

2/10 * 11 - 2 - 4/20;

2.0/10.0 * 11.0 - 2.0 - 4.0/20.0;

19 / 3;

2^1000;

sqrt(2);

14 / 6;

sqrt(12);

sqrt(8);

exp(3);

atan(1);

tan(%pi/4);

19.0 / 3;

2.0^1000;

sqrt(2.0);

sqrt(12.0);

exp(3.0);

atan(1.0);

tan(%pi/4.0);

123;

123.0;

ustawienie precyzji wyœwietlania

fpprintprec: 4$

sqrt(2.0);

fpprintprec: 0$

sqrt(2.0);

/* clear working space */
 reset()$ kill(all)$

### Wyniki obliczeñ - 

ev(19/3, numer:true);

19/3, numer;

sin(4), numer;

%pi, numer;

float(19/3);

float(sin(4));

float(%pi);

obliczanie silni

100!;

ustawianie precyzji obliczeñ

fpprec: 500$

bfloat(%pi);

reset()$

set_display(ascii)$

100!;

fpprec: 500$

bfloat(%pi);

reset()$

set_display(xml)$

/* clear working space */
 reset()$ kill(all)$

### Listy - struktury przechowywania danych 

L: [1,4,9,16,25,[36,49]];

emptylist: [];

L[4];

length(L);

L[3]: -99;

L;

Dopisanie elementu do listy

append(L,[1,2,3]);

Usuwanie elementu z listy

delete(1,L);

L: delete(4,L);

L;

/* clear working space */
 reset()$ kill(all)$

### £añcuchy

s1: "przyk³ad ³añcucha.";

/* clear working space */
 reset()$ kill(all)$

### Rozwi¹zywanie równañ 

solve(x^2+3*x-1=0, x);

float(%);

sol: solve(x^2+3*x-1, x);

ev(x, sol[1]);

rhs(sol[2]);

solve(exp(f(x))=10^y, f(x));

eq1: 3*x^2-y^2=2;

eq2: x^2+y^2=2;

solve([eq1,eq2], [x,y]);

s: solve(x^3+2*x^2-5*x-6=0, x);

z: ev(%pi^x, s[2]);

solve([x^(-1/2)*y^(1/2)=1, x^(1/2)*y^(-1/2)=1], [x,y]);

ev([x^(-1/2)*y^(1/2)-1, x^(1/2)*y^(-1/2)-1], [x=1,y=1]);

solve(x^5-x^4+2*x^3+x^2-x+5, x);

allroots(x^5-x^4+2*x^3+x^2-x+5);

algsys([x^5-x^4+2*x^3+x^2-x+5],[x]);

/* clear working space */
 reset()$ kill(all)$

### Liczby zespolne

solve(x^2-4*x+13=0,x);

s: solve(x^3-x+1/3=0,x);

trigsimp(rectform(s));

float(%);

float(rectform(s));

realpart(%);

/* clear working space */
 reset()$ kill(all)$

### Pakiety dodatkowe

pakiet do obliczeñ na wektorach

load(vect);

/* clear working space */
 reset()$ kill(all)$

### LISP

## Obliczenia na symbolach

### Obliczenia symboliczne

solve(a*x^2+b*x+c=0,x);

solve((s-t)*f(y)^2+abs(s+t)*f(y)+(s+t)=0, f(y));

foo(3*bar+2*bar) + bar^2 + log(%e);

x: (a+beta)*(a-beta);

y: x-a^2;

foo(u):= log(u)-x+alpha;

foo(y);

### Rozwijanie, przekszta³canie i upraszczanie wyra¿eñ

(x+y)*(x-y);

expand(%);

factor(%);

(x+y)^4;

expand(%);

factor(%);

factor(-8*y-4*x+z^2*(2*y+x));

factor(20!);

e: (x+y)^2/(x^2-y^2);

expand(e);

ratexpand(e);

(x^2-y^2)/(x+y)^2;

ratsimp(%);

(%e^x-1)/(1+%e^(x/2));

radcan(%);

num(a/b);

denom(a/b);

(3+cos(x)^4+8*sin(x)+4*(cos(x)^2-sin(x)^2)
  -6*cos(x)^2*sin(x)^2+sin(x)^4)/(8*cos(x)^3);

trigsimp(%);

/* clear working space */
 reset()$ kill(all)$

### Sztuczki

y: 2;

'y^2 = y^2;

f(x):= x^2;

f(a);

'f(a);

f(x):= x^2+y^2;

f(x):= x^2+''y^2;

f(x):= exp(-x^2);

fx(x):= ''(diff(f(x),x));

f(x):= exp(-x^2);

fx(x):= ''diff(f(x),x);

y: 2$

f(x):= exp(-x^2)$

f(y);

'f(y);

f('y);

'(f(y));

/* clear working space */
 reset()$ kill(all)$

### Operatory Logiczne

3>4;

is() - sprawdza wartoœæ logiczn¹ wyra¿enia

is(3>4);

not(3>4);

3>4 or 7<8;

is(x>0);

a: (x+1)*(x-1);

b: x^2-1;

is(a=b);

is(equal(a,b));

/* clear working space */
 reset()$ kill(all)$

### Conditional Evaluation

fabs(x):= if (x>=0) then x else -x;

/* clear working space */
 reset()$ kill(all)$

### Za³o¿enia  

is(x>=0);

is(U(x)>=0);

assume(x>0,U(x)>0);

is(x>=0 and U(x)>=0);

facts pokazuje listê uczyninych za³o¿eñ

facts();

assume(y>0);

facts();

forget(U(x)>0);

facts();

assume(a=1);

assume(equal(a,1));

assume(x+y<1);

facts();

is(x*y>0);

is(x<1);

facts();

assume(x>0);

assume(x<0);

facts();

declare(n,integer);

facts();

/* clear working space */
 reset()$ kill(all)$

### Ustawienie precyzji wewn¹trz wyra¿eñ

ev(bfloat(%e), fpprec:100);

bfloat(%e), fpprec:100;

/* clear working space */
 reset()$ kill(all)$

## Kreœlenie 

### Rysunki 2D

wxplot2d(sin(x), [x,-%pi,%pi]);

wxplot2d(sin(x), [x,-%pi,%pi], [y,-2,2]);

wxplot2d([discrete, [10,20,30,40,50], [.6,.9,1.1,1.3,1.4]]);

wxplot2d([discrete, [[10,.6],[20,.9],[30,1.1],[40,1.3],[50,1.4]]],
  [style,points]);

wxplot2d([parametric, t*sin(2*%pi*t), t*cos(2*%pi*t), [t,0,2],
  [nticks, 1000]]);

wxplot2d([x^3, [parametric,cos(t),sin(t),[t,-%pi,%pi]]],  [x,-1,1]);

kreœlenie funkcji uwik³anych

load(implicit_plot)$

implicit_plot(x^2-y^3+3*y-1=0, [x,-4,4], [y,-4,4]);

/* clear working space */
 reset()$ kill(all)$

### Kreœlenie w 3D

wxplot3d(2^(-x^2+y^2),[x,-3,3],[y,-2,2]);

Kreœlenie wstêgi Mobusa  zadanej parametrycznie 

moeb_x: cos(s)*(3+t*cos(s/2))$

moeb_y: sin(s)*(3+t*cos(s/2))$

moeb_z: t*sin(s/2)$

wxplot3d([moeb_x,moeb_y,moeb_z], [s,-4,4], [t,-4,4]);

funkcja kreœlona w postaci konturu

wxcontour_plot(u^3+v^2, [u,-4,4], [v,-4,4], [legend,false],  [gnuplot_preamble,"set cntrparam levels 12"]);

/* clear working space */
 reset()$ kill(all)$

### Opcje polecenia plot2d, wxplot2d itd

wxplot2d([exp(-x^2/2), sqrt(%e)*exp(-x)], [x,0,2], 
  [y,0,1.5], [nticks, 100], 
  [legend,"first function", "second function"], 
  [xlabel, "abscissa"], [ylabel, "ordinate"], 
  [style, [lines,10],[lines,1]], [color, magenta,blue]);

Kreœlenie wg danych zapisanych w liœcie [[x1,y1],[x2,y2],..] 
oraz oblicznych wg funkcji zadanej wzorem

data_list: [[0,-4], [0.5,-1], [1,0.25], [1.2,0.4], [1.5,0.1],
  [1.7,0.3], [2,1], [2.2,0.75], [3,1.1]]$

wxplot2d([log(x), [discrete, data_list]], [x,0,5],
  [legend, "theory", "experiment"],
  [style, [line,2], [points]],
  [color, blue, red], [point_type, triangle]);

Ustawienie palety kolorów

wxplot3d(sin(x*y), [x,-2,2], [y,-2,2],
  [palette,false], [color, red, blue]);

wxplot3d(sin(x*y), [x,-2,2], [y,-2,2], [grid, 50,50],
  [mesh_lines_color,false], [palette, [hue, 0.5,0.7,0.6, 0.2]]);

dok³adnoœæ kreœlenia powierzchni

wxplot3d(sin(x*y), [x,-2,2], [y,-2,2], [grid, 5,5]);

/* clear working space */
 reset()$ kill(all)$

## Pakiety wspomagaj¹ce kreœlenie

load(draw)$

### Draw Scenes

### Graphic Objects in 2D

draw2d(explicit(sin(x), x,-%pi,%pi));

draw2d(color=blue, explicit(sin(x), x,-%pi,%pi),
  xrange=[-4,4], yrange=[-2,2]);

draw2d(implicit(x^2-y^3+3*y-1=0, x,-4,4, y,-4,4));

draw2d(color=blue, nticks=1000, 
  parametric(t*sin(2*%pi*t),t*cos(2*%pi*t), t,0,2));

draw2d(color=blue, point_type=filled_circle, point_size=2, points([[10,.6],[20,.9],[30,1.1],[40,1.3],[50,1.4]]));

draw2d(color=blue, point_type=none, points_joined=true, points([10,20,30,40,50],[.6,.9,1.1,1.3,1.4]));

draw2d(color=blue,key="x^3", explicit(x^3, x,-1,1),
  color=red,key="circle", parametric(cos(t), sin(t), t,-%pi,%pi));

/* clear working space */
 reset()$ kill(all)$

### Graphic Objects in 3D

draw3d(explicit(2^(-x^2+y^2), x,-3,3, y,-2,2));

draw3d(explicit(2^(-x^2+y^2), x,-3,3, y,-2,2), surface_hide=true);

draw3d(enhanced3d=true, explicit(2^(-x^2+y^2), x,-3,3, y,-2,2));

draw3d(implicit(x^2+sin(y)=z^2, x,-2,2, y,-2,2, z,-2,2));

draw3d(nticks=500, parametric(sin(t),cos(t),t, t,0,6*%pi));

moeb_x: cos(s)*(3+t*cos(s/2))$

moeb_y: sin(s)*(3+t*cos(s/2))$

moeb_z: t*sin(s/2)$

draw3d(parametric_surface(moeb_x,moeb_y,moeb_z, s,-%pi,%pi, t,-1,1));

draw3d(explicit(x^3+y^2, x,-4,4, y,-4,4), 
contour=map, contour_levels=12);

draw3d(explicit(x^2+y^2, x,-4,4, y,-4,4), 
contour=both, contour_levels=15);

xc:[1,2,3,4,5]$

yc:[5,4,2,7,4]$

zc:[2,4,3,1,6]$

draw3d(color=red, point_type=filled_circle, points(xc,yc,zc));

/* clear working space */
 reset()$ kill(all)$

### Graphic Options

set_draw_defaults(color=red, grid=true)$

draw2d(explicit(cos(x),x,-2*%pi,2*%pi),
  terminal='png, file_name="cosine");

draw2d(nticks=100, yrange=[0,1.5], 
  xlabel="abscissa", ylabel="ordinate", 
  color=magenta, line_width=10, key="first function", 
  explicit(exp(-x^2/2), x,0,2), 
  color=blue, line_width=1,key="second function", 
  explicit(sqrt(%e)*exp(-x), x,0,2));

datalist: [[0,-4], [0.5,-1], [1,0.25], [1.2,0.4], [1.5,0.1],
   [1.7,0.3], [2,1], [2.2,0.75], [3,1.1]]$

draw2d(color=blue, key="theory", explicit(log(x), x,0,5), 
  color=red, point_type=filled_up_triangle, point_size=3,  key="experiment", points(datalist));

draw3d(key="Gauss", color=blue, 
  explicit(20*exp(-x^2-y^2)-10,x,-3,3,y,-3,3), 
  key="Plane", color=red, yv_grid=10, 
  explicit(x+y,x,-5,5,y,-5,5), 
  key="Circle", color=dark_green, line_width=3, 
  parametric(2*sin(t),2*cos(t),1, t,0,2*%pi), 
  surface_hide=true);

/* clear working space */
 reset()$ kill(all)$

## Algebra liniowa

### Wektory

[1,0,1]+[1,2,3];

[1,0,1]-[1,2,3];

5*[1,2,3];

[1,0,1].[1,2,3];

[1,2,3]*[1,2,3];

x: [1,2,3,4];

x[3];

y[3];

z: [z[1],z[2]];

/* clear working space */
 reset()$ kill(all)$

### Macierze

A: matrix([1,2],[3,4]);

B: matrix([0,1],[2,3]);

A+B;

5*A;

A.B;

A^^2;

A^^(-1);

A*B;

A^2;

A: matrix([1,2],[3,4]);

A[2];

A[2,1];

A[2,1]: 10;

A;

A: matrix([1,2],[3,4]);

x: [1,2];

A.x;

x.A;

/* clear working space */
 reset()$ kill(all)$

### Operacje na macierzach

\*  macierz jednostkowa *\;

ident(2);

/*  macierz diagonalna */;

diagmatrix(2,3);

/*  definiwanie macierzy */;

A: matrix([1,2],[3,4]); 

/* transpozycja */;

transpose(A);  

/* odwracanie macierzy */

invert(A);

/* wyznacznik macierzy*/

determinant(A);

/* clear working space */
 reset()$ kill(all)$

### Rz¹d macierzy i inne 

A: matrix([1,2,3],[4,5,6],[7,8,9]);

rank(A);

/* clear working space */
 reset()$ kill(all)$

### Lista, dodatki

first(f(x,y,z));

rest(f(x,y,z));

last(f(x,y,z));

A: matrix([1,2,3],[4,5,6],[7,8,9]);

CS: columnspace(A);

first(CS);

second(CS);

/* clear working space */
 reset()$ kill(all)$

### Równania liniowe

A: matrix([1,2,2],[3,3,8],[0,7,4]);

b: [-1,-15,2];

A^^(-1).b;

eq1: x[1]+2*x[2]+2*x[3]=-1;

eq2: 3*x[1]+3*x[2]+8*x[3]=-15;

eq3: 7*x[2]+4*x[3]=2;

solve([eq1,eq2,eq3], [x[1],x[2],x[3]]);

solve([eq1,eq2], [x[1],x[2],x[3]]);

%, %r1=2;

/* clear working space */
 reset()$ kill(all)$

### Wartoœci i wektory w³asne macierzy

A: matrix([1,2,1],[1,4,4],[7,3,5]);

charpoly(A,x);

B: matrix([1,0],[6,6]);

eigenvalues(B);

eigvals: eigenvalues(B)[1];

multiplicities: eigenvalues(B)[2];

es: eigenvectors(B);

eigvals: es[1][1];

multiplicities: es[1][2];

vectors: es[2];

C: matrix([3,1],[0,3]);

eigenvalues(C);

eigenvectors(C);

D: matrix([3,0],[0,3]);

eigenvectors(D);

/* clear working space */
 reset()$ kill(all)$

## Obliczenia

### Granice

limit(sin(x)/x,x,0);

limit(atan(x),x,minf);

limit(log(x),x,0);

limit(sin(1/x),x,0);

limit(sin(1/x)/x,x,0);

limit(log(x),x,0,plus);

limit(x/abs(x),x,0,plus);

limit(x/abs(x),x,0,minus);

limit(x/abs(x),x,0);

/* clear working space */
 reset()$ kill(all)$

### Ciagi, sumy,iloczyny

sum(z(i),i,1,5);

sum(z[i],i,1,5);

sum(i^3,i,1,10);

sum(z[i],i,1,n);

sum(z[i],i,1,inf);

sum(1/4^i,i,1,inf);

sum(1/4^i,i,1,inf), simpsum;

sum(1/i^2,i,1,inf);

sum(1/i^2,i,1,inf), simpsum;

'product(i^3,i,1,10);

product(z(i),i,1,minf);

product(x+i*(i+1)/2,i,1,4);

/* clear working space */
 reset()$ kill(all)$

### Ró¿niczkowanie 

diff(x^a,x);

diff(x^a,x,2);

diff(diff(x^3*y^2,x),y);

diff(x^3*y^2,x,1,y,1);

diff(x^3*y^2,x,2,y,1);

diff(sin(y*x));

diff(sin(y*x));

%, del(x)=dx, del(y)=dy;

'diff(exp(-x^2),x)=diff(exp(-x^2),x);

/* clear working space */
 reset()$ kill(all)$

### Ró¿niczkownaie rozszerzenia

diff(f(x,y),x);

diff(f,x);

depends(f,x);

diff(f,x);

derivabbrev: true$

diff(f,x);

derivabbrev: false$

diff(f(x(t)),t);

depends(x,t);

diff(f,t);

dependencies;

kill(f)$

dependencies;

depends(f,[x,t]);

depends(x,t);

diff(f,x);

diff(f,t);

kill(all)$

diff(a*x^2+b*y^2);

declare([a,b],constant);

diff(a*x^2+b*y^2);

/* clear working space */
 reset()$ kill(all)$

### Jakobian i gradient

f(x,y):= [x^2+y^2, x^2-y^2];

jacobian(f(x,y),[x,y]);

f(x,y):= [f[1](x,y),f[2](x,y)];

jacobian(f(x,y),[x,y]);

jacobian([x^3+y^2],[x,y]);

/* clear working space */
 reset()$ kill(all)$

### Szereg Taylora i inne

taylor(exp(x),x,0,3);

T(x):= ''(taylor(exp(x),x,0,3));

T(3);

T(y);

taylor(exp(x^2+y),[x,y],[0,1],2);

powerseries(sin(x),x,0);

niceindices(%);

/* clear working space */
 reset()$ kill(all)$

### Minima and Maxima of Univariate Functions

f(x):= x^4-3*x^2+2;

plot2d(f(x), [x,-2,2], [y,-1,3]);

fx: diff(f(x),x);

crit: solve(fx=0,x);

fxx: diff(f(x),x,2);

fxx, crit[1];

fxx, crit[2];

fxx, crit[3];

/* clear working space */
 reset()$ kill(all)$

### Minima i maxima funkcji wielu zmiennych

hessian(x^3+4*x*y+y^2,[x,y]);

depends(g,[x,y])$

hessian(g,[x,y]);

f(x,y):=1/6*x^3 - x + 1/4*x*y^2;

wxplot3d(f(x,y), [x,-2,2], [y,-2.5,2.5]);

wxcontour_plot(f(x,y), [x,-2,2], [y,-2.5,2.5], [legend,false],
  [gnuplot_preamble, "set cntrparam levels 12"]);

[fx,fy]: [diff(f(x,y),x), diff(f(x,y),y)];

crit: solve([fx=0,fy=0],[x,y]);

hess: hessian(f(x,y),[x,y]);

H1: ev(hess, crit[1]);

eigenvalues(H1);

H3: ev(hess, crit[3]);

eigenvalues(H3);

H2: ev(hess, crit[2]);

[H2[1,1], determinant(H2)];

H4: ev(hess, crit[4]);

[H4[1,1], determinant(H4)];

/* clear working space */
 reset()$ kill(all)$

### Ca³kowanie 

integrate(sin(x)^3, x);

integrate(x/sqrt(b^2-x^2),x);

integrate(exp(-x^2),x);

integrate(tan(log(x)),x);

integrate(x^a,x);

assume(equal(a,-1));

integrate(x^a,x);

integrate(cos(x)^2*exp(x),x,0,%pi);

integrate(x^2*exp(-x^2),x,minf,inf);

metoda prostokatów

quad_qags(x^(1/2)*log(1/x),x,0,1);

/* clear working space */
 reset()$ kill(all)$

## Równania ró¿niczkowe (ODE)

### Rozwiazania analityczne równañ ró¿niczkowych

x^2*'diff(y,x) + 3*y*x = sin(x)/x;

x^2*diff(y(x),x) + 3*y(x)*x = sin(x)/x;

ode2(x^2*'diff(y,x) + 3*y*x = sin(x)/x, y, x);

method;

ode2(x^2*diff(y(x),x) + 3*y(x)*x = sin(x)/x, y(x), x);

sol1: ode2(x^2*'diff(y,x) + 3*y*x = sin(x)/x, y, x);

ic1(sol1, x=%pi, y=0);

sol2: ode2('diff(y,x,2) + y*'diff(y,x)^3 = 0, y, x);

ic2(sol2, x=0, y=0, 'diff(y,x)=2);

/* clear working space */
 reset()$ kill(all)$

### Uk³ady równañ ó¿niczkowych

desolve(diff(y(x),x,2)+y(x)=2*x, y(x));

atvalue(y(x), x=0, 1)$

atvalue(diff(y(x),x), x=0, 2)$

desolve(diff(y(x),x,2)+y(x)=2*x, y(x));

kill(y)$

eqn1: diff(x(t),t)=x(t)-2*y(t);

eqn2: diff(y(t),t)=-x(t)+y(t);

desolve([eqn1,eqn2],[x(t),y(t)]);

atvalue(x(t), t=0, 1)$

atvalue(y(t), t=0, 1)$

desolve([eqn1,eqn2],[x(t),y(t)]);

/* clear working space */
 reset()$ kill(all)$

### Portret fazowy

load(plotdf)$

plotdf(exp(-x)+y)$

plotdf(exp(-x)+y, [trajectory_at,2,-0.1])$

plotdf(exp(-x)+y, [trajectory_at,2,-0.1], [direction,forward], [versus_t,1])$

/* clear working space */
 reset()$ kill(all)$

plotdf([x-2*y, -x+y], [trajectory_at,1,1], [versus_t,1])$

/* clear working space */
 reset()$ kill(all)$

## Æwiczenia z ekonometrii

### Monopoly

p(q):= 124-6*q;

MR(q):= ''(diff(q*p(q),q));

MC(q):= 40;

eq1: MR(q)=MC(q);

q0: ev(q, solve(eq1,q));

p0: p(q0);

pi(p,q):= p*q-40*q;

pi(p0,q0);

plot2d([MR(q),p(q),MC(q)],[q,0,11],
  [color,green,blue,red], [legend, "MR(q)","p(q)","MC(q)"]);

/* clear working space */
 reset()$ kill(all)$

### Profit Maximization

p(q):= 50-1/2*q;

TC(q):= 2+60*q-8*q^2+q^3;

TR(q):= ''(expand(p(q)*q));

TP(q):=''(expand(TR(q)-TC(q)));

cp: solve(diff(TP(q),q)=0,q);

float(cp);

[TP(0),TP(rhs(cp[1])),TP(rhs(cp[2])),limit(TP(q),q,inf)], numer;

q0: rhs(cp[2]), numer;

TP0: TP(q0);

p0: p(q0);

plot2d([TP(q),diff(TP(q),q),TC(q)/8,TR(q)/8], [q,0,6],
  [legend, "TP(q)", "MP(q)", "TC(q)", "TR(q)"]);

/* clear working space */
 reset()$ kill(all)$

### Price Discrimination

q: [q[1], q[2]];

p(q):= [20-2*q[1], 10-q[2]];

C(q):=(q[1]+q[2])^2+5;

TP(q):=p(q).q-C(q);

MP(q):= ''([diff(TP(q),q[1]), diff(TP(q),q[2])]);

cq: solve(MP(q), q);

hessian(TP(q), q);

float(eigenvalues(%));

q0: ev(q, cq);

MR(q):= ''([diff((p(q)*q)[1],q[1]), diff((p(q)*q)[2],q[2])]);

MC(q):= ''([diff(C(q),q[1]), diff(C(q),q[2])]);

MR(q0); MC(q0);

/* clear working space */
 reset()$ kill(all)$

### Cournot Duopoly

q: [q[1],q[2]]$

p(q):= a-b*(q[1]+q[2]);

c:[c[1],c[2]]$  d:[d[1],d[2]]$

C(q):= d+c*q;

pi(q):= p(q)*q - C(q);

R: [solve(diff(pi(q)[1],q[1]),q[1])[1],
  solve(diff(pi(q)[2],q[2]),q[2])[1]];

solve(R,q);

/* clear working space */
 reset()$ kill(all)$

### Cobb-Douglas Production Function

plot3d(3*x^(1/2)*y^(1/3), [x,0,5], [y,0,5]);

u(x,y):= c*x^alpha*y^beta;

hess: hessian(u(x,y), [x,y]);

factor(determinant(hess));

/* clear working space */
 reset()$ kill(all)$

### Stationary Points and Cobb-Douglas Function

f(x,y):= 6*x^(1/2)*y^(1/3)-3*x-2*y;

plot3d(f(x,y), [x,0,3], [y,0,3], [z,0,1]);

fx: diff(f(x,y), x);

fy: diff(f(x,y), y);

solve([fx=0,fy=0],[x,y]);

assume(x>0,y>0);

solve(fx=0,x);

yo: solve(ev(fy,%)=0,y);

xo: solve(ev(fx,yo)=0,x);

kill(all)$

f(x,y):= c*x^a*y^b-r*x-w*y;

fx: diff(f(x,y), x);

fy: diff(f(x,y), y);

assume(x>0,y>0)$

solve([fx=0,fy=0],[x,y]);

solve(fx=0,x);

yo: solve(ev(fy,%)=0,y);

xo: solve(ev(fx,yo)=0,x);

/* clear working space */
 reset()$ kill(all)$

### Constraint Optimization

Y(K,L):= K^2*L;

g(K,l):= K+L;

F(K,L,mu):= ''(Y(K,L) + mu*(3-g(K,L)));

gradF: jacobian([F(K,L,mu)],[K,L,mu])[1];

cp: solve(gradF);

H: matrix([0, diff(g(K,L),K), diff(g(K,L),L)],
  [diff(g(K,L),K), diff(F(K,L,mu),K,2), diff(F(K,L,mu),K,1,L,1)],
  [diff(g(K,L),L), diff(F(K,L,mu),L,1,K,1), diff(F(K,L,mu),L,2)]);

ev(determinant(H), cp[1]);

ev(determinant(H), cp[2]);

cp[1];

/* clear working space */
 reset()$ kill(all)$

kill(all);