Finite Difference Approximation, Mathematica code

Click For Summary
SUMMARY

The discussion centers on programming a three-component decay chain using finite difference approximation in Mathematica. The user initially faced an issue with the magnitude of the y-axis values being too high, despite the curve appearing correct. The governing differential equation for the second chain is correctly stated as dN2/dt = -λ2N2 + λ1N1. Ultimately, the user discovered that the analytical solution was plotted incorrectly, confirming that the numerical solution was accurate.

PREREQUISITES
  • Understanding of finite difference approximation
  • Familiarity with differential equations, specifically decay chains
  • Proficiency in Mathematica programming
  • Knowledge of numerical methods for solving differential equations
NEXT STEPS
  • Study the implementation of finite difference methods in Mathematica
  • Learn about numerical stability in solving differential equations
  • Explore optimization techniques for Mathematica code
  • Investigate analytical versus numerical solution comparisons in decay problems
USEFUL FOR

Students and researchers in computational physics, mathematicians working with numerical methods, and developers looking to optimize Mathematica code for solving differential equations.

timman_24
Messages
52
Reaction score
0

Homework Statement



I have to program a three component decay chain using finite difference approximation. I understand finite difference and have written my code, but I have an error I can not find which is giving me an erroneous answer. The curve is correct, but the magnitude of the y-axis values is too high. I've gone through it for hours but can not find the mistake.

2. The attempt at a solution

It should be pretty basic stuff, but there is something screwy in there. My problem starts with the second chain (N2a and N2b.)

The governing differential equation for that section is:
\frac{dN_{2}(t)}{dt}=-\lambda _{2}N_{2}(t)+\lambda _{1}N_{1}(t)

My forward approximation is:
N_2b=(\lambda_{1} N_1a - \lambda_{2} N_2a) (tchange) + N_2a

For some reason result of this is far too large, but the curve looks to be correct. Here is an image of the real versus approx graph I get:

Approx:
[PLAIN]http://img19.imageshack.us/img19/9539/n2approxgraph.png
Real (analytical):
[PLAIN]http://img222.imageshack.us/img222/5544/n2realgraph.png

Here is a portion of the code in question:

Code:
lambdaA = .6931;
lambdaB = .13863;
tchange = .05;
ttotal = 75;
A0 = 100;
B0 = 0;
C0 = 0;

timesteps = ttotal/tchange;
list = Table[{i*tchange, 0}, {i, timesteps}];
list[[1, 2]] = A0;
list2 = Table[{i*tchange, 0}, {i, timesteps}];
list2[[1, 2]] = B0;
list3 = Table[{i*tchange, 0}, {i, timesteps}];
list3[[1, 2]] = C0;
For[i = 1, i < timesteps, i++,
 N1a = Evaluate[Extract[list, {i, 2}]];
 N1b = Evaluate[N1a (1 - tchange*lambdaA)];
 list[[i + 1, 2]] = N1b;
 N2a = Extract[list2, {i, 2}];
 N2b = Evaluate[(lambdaA N1a - lambdaB N2a) tchange + N2a];
 list2[[i + 1, 2]] = N2b;
 N3a = Extract[list3, {i, 2}];
 N3b = Evaluate[tchange*lambdaB*N2b + N3a];
 list3[[i + 1, 2]] = N3b;]

Thanks a lot!
 
Last edited by a moderator:
Physics news on Phys.org
Oh my god! Its not wrong! I plotted the analytical solution incorrectly. How about that? The numerical solution being correct but the analytical being wrong. Ugh!

Well, anyway.

Anyone see any ways to streamline the code?
 

Similar threads

  • · Replies 1 ·
Replies
1
Views
3K
Replies
1
Views
2K
  • · Replies 10 ·
Replies
10
Views
2K
  • · Replies 2 ·
Replies
2
Views
2K
  • · Replies 1 ·
Replies
1
Views
2K
  • · Replies 4 ·
Replies
4
Views
3K
  • · Replies 3 ·
Replies
3
Views
3K
  • · Replies 1 ·
Replies
1
Views
2K
  • · Replies 1 ·
Replies
1
Views
2K
  • · Replies 2 ·
Replies
2
Views
2K